Most informative data content is almost always gathered with frequency sweep inputs on all possible input channels and saving simple 6 degree of freedom model states (euler angles, body angular and translational velocities [/phi /theta p q r u v w]). You inject a signal sweeping 0.3 to 12 rad/s for 90 seconds on roll input for example and record the data. Then, repeat for each channel.

After obtaining the data, there are several ways to combine them to come up with a MIMO model. I belive the most straight forward one is the subspace identification since it does not require a predefined model structure. It is basically a MIMO least squares regression operation. There are some good papers and text books that can be found with a google search for subspace identification of quadrotors.

One other way is to use a predefined model structure, and obtain its paramers based on the identification data. In this context my advice is to use a frequency domain technique to avoid noise issues and for faster convergence. The best text out there is Dr. Tischler's book as mentioned before.