Zaher El Chami1

Antoine Dinh-

Tuan Pham2

Christine Serviexe3
email: christine.serviere@inpg.fr

Alexandre Guerin1
email: alexandre.guerin@orange-ftgroup.com

A NEW EM ALGORITHM FOR UNDERDETERMINED CONVOLUTIVE BLIND SOURCE SEPARATION

This paper presents a new statistical method for separating more than two sound sources from a two-channel record ing. It is based on a probabilistic model of the Interchannel Level/Phase Difference presented in [1] and the model para meters are estimated using the maximum likelihood criterion and an Expectation-Maximization algorithm. The source separation task is achieved by soft time-frequency masking of the observation. These masks are derived from the es timated source position model. Algorithm performance is evaluated on the real and synthetic convolutive mixtures data of the first audio source evaluation campaign [2] as well as the Signal Separation Evaluation campaign (SiSEC) [10]. Promising results are obtained when comparing to the other methods presented in these two campaigns.

INTRODUCTION

Blind Source Separation (BSS) is a widely used technique that aims at recovering a set of N original sources based only on their M observed mixtures. This task is more difficult when the mixing model is not instantaneous but convolutive and gets even harder in the underdetermined case (N > M). Indeed, when N is larger than M, no algebraic linear solution can be found to separate the sources, even if the mixing matrix is identified. Still, with the source sparseness assumption, researchers have found a way to build non linear masks for the separation task.

Most of the BSS methods that consider the sparseness as sumption deal with the two-channel case from which cues or features like Interchannel Level/Phase Difference (ILD /IPD) are used. In [START_REF] Araki | K-means based underdetermined blind speech separation[END_REF] a KMeans algorithm classifies these features in clusters and binary separating masks are estimated. In [START_REF] Michael | EM localiza tion and separation using interaural level and phase cues[END_REF], rather than clustering the observed (ILD /IPD), the authors propose to model these features as Gaussian variables with the assumption of a dominant path. Then, after estimat ing the model parameters, soft separating masks have been derived. Nevertheless, in real-world situations, the underly ing linear phase assumption deriving from the dominant path hypothesis reveals not applicable due to early reflections and reverberation. Also, the proposed Gaussian model, even if it facilitates the equation computation, does not have any theoretical background.

In this paper a Model Based Underdetermined (blind) Source Separation (MBUSS) is considered. In our previous work [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF] a theoretical distribution for the (log(ILD) /IPD) features is presented, but no separation algorithm based on this distribution was explicitly provided. This paper, pro poses an estimation procedure for this theoretical distrib ution parameters, based on an Expectation-Maximization (EM) algorithm where the Maximization step is speeded up by a Quasi-Newton algorithm. Unlike [START_REF] Araki | K-means based underdetermined blind speech separation[END_REF], probabilistic soft masks are computed instead of binary ones, thus less distortion and artifact are audible in the extracted sources. Unlike [START_REF] Michael | EM localiza tion and separation using interaural level and phase cues[END_REF], separation is performed independently in each fre quency band, therefore no linear phase assumption is made and wide band sources can be extracted even with the pres ence of frequency aliasing. Of course separating indepen dently in each frequency band has the traditional permu tation problem drawback. This permutation alignment is corrected using the ratio envelope of the extracted sources as in [START_REF] Sawada | Measuring dependence of bin-wise separated signals for permutation alignment in frequency-domain BSS[END_REF].

MODEL BASED SOURCE SEPARATION

We describe here briefly the (log(ILD) /IPD) probabilistic model that will be used in the separation task. For more detailed analysis, please refer to [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF]. Consider the convolutive two-channel mixture model:

N N Xj (t) = L xj i ) (t) = L L aji (k) Si (t -k) (1)

i=l i=1 k

where xj i >Ct) is the contribution of the i th source to the j th sensor. si(t), i = 1, .. ,N are the sources and aji(k) is the impulse response of the acoustic channel separating source i from microphone j with j = 1, 2

£ -1 N Xj(t,w) = L w(k)Xj(t + k)e -j w k = LX; i )(t,W) (2) k=O i=l
where w(k) is a window (e.g. Hanning) and X?)(t,w) is the STFT of xj i )(t). Sparseness of the sources in the Time Fre quency (TF) domain is the key assumption in solving the underdetermined separation problem. It means that each given source is non negligible on only a few number of TF slots. It generally implies that the sources have nearly dis joint supports in the TF domain in the sense that for each TF slot there can be at most one dominant source, and it is this last assumption that will be actually assumed [START_REF] Araki | K-means based underdetermined blind speech separation[END_REF]. If q is the index of this dominant source at (t, w), then the j th observa tion in eq. ( 2) can be approximated to Xj (t,w) :::::; x; q ) (t,w) and the ratio between the two observations at (t, w) will then be:

The approximation X?)(t,w) :::::; Aji(w)Si(t,W), with Aji(W) as the Fourier transform of aji(k), can be found in most of the TF source separation methods. With it, the last term in

(3) can be reduced to a constant in each frequency band w and thus: R(t,w) = X1(t,w) � xi q )(t,w) � A1q(w) . X2(t,w) x� q >Ct,w) A2q(W)

The above ratio, being (approximately) time independent but frequency and source dependent, has been widely used as dominant source indicator in the TF binary masks ap proach. In [START_REF] Araki | K-means based underdetermined blind speech separation[END_REF], the modulus and argument of R(t,w), which are no more than the well-known Interchannel Level Differ ence (ILD) and Interchannel Phase Difference (IPD), have been clustered by a KMeans algorithm. The cluster members define the binary separating mask and the cluster centers give an estimation of Alq(W)/A2q(W). But does the above approximation (4) hold when dealing with a long tap im pulse response, i.e with a reverberating environment?

1 One source Ratio STFT distribution model

To simplify the notations in this section, we will omit the w parameter. From now, we implicitly work in given frequency band w. It has been demonstrated in [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF] that, even if one source is observed, R(t) is not constant in time. Fig. 1 plots the real and imaginary part of R(t) for a 10-second speech source placed in a moderate reverberant conditions (T60 = 250ms). This one source observed ratio should be considered as a (complex) random variable, the distribution of which has been derived in [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF]. More precisely this paper provides the theoretical joint density of the real part and imaginary parts

x and y of 10gR( q )(t) = 10g[ Xi q )(t,w)/X� q )(t,w)] :

x = 10gIR( q )(t)1 = 10gILDq, y = arg R( q )(t) = IPDq
where ILDq/IPDq are the ILD/IPD of Sq. Note that we did not considered the traditional feature R( q )(t) because it admits an infinite variance (see [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF]). To estimate the joint density of (x, y), we assume that L is large enough so that, by the Central Limit Theorem, the pair xi q )(t) and x� q )(t) can be considered as Gaussian (complex circular) [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF] with variance O"i ,q (t) and O"� ,q (t) and complex cross-correlation f3q. It has been demonstrated in [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF] that the variance ratio 0"1,q/0"2,q and f3q does not vary in time and may be considered as specific features for each source S q , hence of the variable log R( q ) (t). Based on the complex gaussian circular assump tion, the couple of variables (x, y) admits the following joint density [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF]: p (x, yl rq, pq) = Ppq (x -log Irql, yarg rq) [START_REF] Sawada | Measuring dependence of bin-wise separated signals for permutation alignment in frequency-domain BSS[END_REF] where rq = (0"1,q/0"2,q)e i arg(3Q , pq = lf3ql and 1

1 _ p 2 pp (x,y) = 47r ( cosh x -p cos x) 2
The parameters r q and pq are specific to the source Sq posi tion in space. Note that it is the source position in space which is modeled and not the source itself. Thus, this "'space'" position model can be applied independently from the source model (gaussian, laplacian, ... ). More in details, r q corresponds to the mean of log R( q ) (t) and thus depends only on the source position in space where Irql is equal to the ILDq and arg ( r q) is equal to the IPDq. As for Pq, it will stand for the reverberation degree of the acoustic path separating Sq from the set of microphones. This can be viewed from its definition as the modulus of the cross-correlation between the two observations where high reverberation causes low cross-correlation between microphones and vice-versa. For example, in free field conditions, i.e anechoic environment, the cross-correlation is maximum and pq = 1. ------------'1 ---------:-:---,-:--,-, ,-, -, arg(R(t)) --arg( A,JA, q ) = .

� o ���� � �� , . �� ,

(6) i=l

This model is given at the frequency w for the set of consid ered time points t E T with p = [ Pl, .. ,PN], r = [ rl, .. , rN] and J.l = [ J-ll,", J-lN] where J-li is an added parameter that denotes the a priori probability of the i th source in the con sidered frequency band. This parameter reveals necessary since all sources are not equiprobable in a given frequency band, depending for instance on the mean pitch of a per son. Simulations with equiprobable hypothesis showed that the estimation of the model parameters p and r is biased, hence producing degraded performance in terms of source separation.

Soft mask separation

The set of parameters p, r and J.l, which depend on the frequency band w, are the parameters of the mixture model and need to be estimated in order to separate the sources.

Once the parameters of the probabilistic model given in [START_REF] Nesta | A novel robust solution to the permutation problem based on a joint multiple TDOA estimation[END_REF] are estimated, the a posteriori probability that i th source is dominant at the TF point (t, w) can be obtained directly as follows:

(7)
where Ri(t) = R(t)/ri. Then, source separation can be eas ily performed in each frequency band w by applying these above probabilities to the observations. Note that sources are independently separated in each frequency band, thus permutation ambiguity remains and needs to be solved. In this paper, the correlation between the ratio envelope is used as in [START_REF] Sawada | Measuring dependence of bin-wise separated signals for permutation alignment in frequency-domain BSS[END_REF]. Note that other methods, based on linear phase as sumption [START_REF] Nesta | A novel robust solution to the permutation problem based on a joint multiple TDOA estimation[END_REF], can not be used due to reverberation. Fig. 2 shows the flow chart of the Model Based Source Separation approach where Mi (t, w) refer to the masks that extract the i th source. Mi(t,W) is constructed from the a posteriori probabilities (7) after correcting the permutation ambigui ties. Thus, Mi(t,w) = 7rrrw ( i )

with IIw = [IIw(1), ... , IIw(N)] t: "I (f, m) M I (f,m) J\(f, m) lJ (t) ,2 c 'iii
,- As said before, to build the soft separating mask, the Mi, ri, and Pi parameters need to be estimated for each source at each frequency band w. To do so, we propose to use as criterion the maximum of the log-likelihood of the data {log IR(t) I, arg R(t)} , t E T. Under the independence assumption between the time set of observations, it is given by: L(p,r,ll) = � log { � '71WpJ l O g IRi(t)1 ,arg Ri (t)] } . (8) The above log-likelihood is hard to maximize. However, it may be recasted as the log-likelihood for a model with miss ing observations or hidden variables. These variables are the indexes that indicate which source is dominant at each (t, w) point (here w is fixed and hence not displayed). In this con text the log-likelihood can be maximized by the well known Expectation Maximization (EM) algorithm [START_REF] Dempster | Maxi mum likelihood from incomplete data via the EM algo rithm[END_REF],

:J xl (f) . § § .2 «: t: en ,2 U ...l x2 (t) '" c 0 � '" � ::J f- o(! "- ::J :; 0 en xi (f) f- 'ol Q. U c f- r/l '" § c "'" "- "- ,2 en "- co XM(t xM (t,W) v u '" � '" '" "'" � en en S "- '" co '" 5 f- > CO � .s " "- "N (f, m) cc. MN(r,m Y,v (fm) YN(r) v f- 0 �

THE EM ALGORITHM

This algorithm operates in two steps as described below.

1 The E-step

This step computes the conditional expectation of the full log-likelihood given the data {log IR(t)l,arg R(t)},t E T. The expected log-likelihood will be computed at generic new parameters M�, r� and p� and the conditional expectation is computed relatively to the model specified by the current parameters Mi, ri and Pi. The result can be shown to be: (

"'� {' [I R(t) 1 R(t) ]}
) 9
where 7ri(t) is the a posteriori probability given in ( 7) and computed at the current parameter Mi, ri and Pi.

The M-step

This step maximizes the above conditional expectation of the full log-likelihood with respect to the generic parameters M�, r� and p�. The maximum point is then taken as the new pa rameter. It is easily seen that the maximization of (9) with respect to M� (under the constraint �:l M� = 1) and with re spect to the set (r�, p�) can be performed independently. The first maximization yields the new Mi:

Mi = � tET [7ri(t) l /ITI
where ITI denotes the number of points in T. The second one is reduced to the maximization of:

I I '" [I R(t) 1 R(t) ]
C (ri ' Pi) = L 7ri (t) log pp ; log ---:;;-,arg ---:;;-tET "

for each i = 1, '" N with respect to r� and p�. Maximizing the above expression (10) is not easy and can not be done without using an iterative algorithm. In [START_REF] Ei-Chami | A new-model based underdetermined source separa tion[END_REF], we overcome this problem by replacing the theoretical den sity ( 5) in (10) with another similar and easier one to handle. Thus, we were able to maximize it analytically. In this work, as the theoretical density shall be used, we propose not to maximize the expected full log-likelihood (10) but to just make it increase. In fact, from the EM algorithm theory, in creasing the expected likelihood (10) would be sufficient to increase, also, the marginal likelihood [START_REF] Ei-Chami | A new-model based underdetermined source separa tion[END_REF]. Therefore, given the current parameters ri and Pi, we limit ourselves to find ing the new parameters r; and P; that increase the objective function (10). Such new parameters can indeed be found analytically as described in the following basic step: and ( = �u 7ri(U) , �i(t) = log IR(t)l, Ai = log Ir;jril, <Mt) = arg Ri(t), (;li = arg(r;jri) and R(z) is the real part of the complex number z. The maximum of the above equa tion (11), with respect to (r�, p�), is strictly larger than zero (because it is equal to zero when (r�,p�) = (ri,Pi)). Thus, the new couple (r;, pD that maximizes (11) will increase, at the same time, the objective function (10). Therefore, by a simple variable substitution, the problem of increasing the likelihood is now reduced to maximize (11) with respect to where sign {z} = z/Izl = e i argz. Finally, to increase (10),

we are led to assume the following one iteration basic step:

where ai = I:t xi( t) IRi( tW, bi = I:t Xi(t), and Ci I:t Xi(t)Ri(t) with Xi(t) = �wi(t)/IRi(t)l•

Combining with the Quasi Newton step

The basic step could lead to a slower convergence of the algo rithm as the obtained increase of the expected log-likelihood (10) can be much less than the maximum achievable. To overcome this problem we will consider a Quasi-Newton (Q N) algorithm. In fact, when the old parameter ri and Pi are close to the maximum solution of (10), the Quasi-Newton algorithm would have a quadratic convergence to the point maximizing (10). Therefore, the algorithm should converge within one single Q-N iteration. However, in contrast to the basic step, the Q-N step does not guarantee the increase of the expected log-likelihood. It does not even guarantee the new P; to be in the interval (0,1). Thus, the following strat egy is adopted in the M-step:

• Compute the new parameter P; and r; based on the basic step as in section 3.3.1

• Compute the other new estimation of p; of the Q-N step and test if it belongs to the interval (0,1)

• If not, adopt the new parameter P; and r; of the basic step, otherwise compute the other new parameter f; of the Q-N step and test if it and p; of the Q-N step lead to a larger increase of the expected log-likelihood than the basic step; if so, adopt these parameter f; and p; otherwise adopt those of the basic step P; and r;. Computations (not detailed here) show that the Q-N step is given by:

Pi -I:t x2( t) {[IRi( t) 12 + IF -4R[Ri(t)J2}/11"i(t) ri exp[(3/2) (ai + bi)/(ai -bi)]sign[3ci -R(Ci)]
Note that these formula involve the already computed quan tities ai, bi, Ci.

EM Initialization and Stop Criteria

A simple way to initialize our algorithm is to choose ran domly N times points tl, .. , tN and initialize ri by ri = R(ti), i = LN. The a posteriori probability 11"i(t) can be initial ized by formula (7) using the density [START_REF] Sawada | Measuring dependence of bin-wise separated signals for permutation alignment in frequency-domain BSS[END_REF] with Pi set to one, leading to:

11"i(t) = :[R(t),ri] I:i=l d[R (t),ri] where d [R (t) ,ri] = cosh log [R (t) /ri] -cos arg [R (t) /ri]'
Since the log-likelihood increases monotonically at each EM iteration, the EM iteration process is stopped when the in crease of log-likelihood becomes insignificant(e.g 10-6 ).

EXPERIMENTS AND RESULTS

In order to evaluate our algorithm, we simulated speech and music mixtures in a reverberant noise-free environment by convolving speech and music samples with filter impulse re sponses coming from the first audio source separation cam paign [START_REF] Vincent | First stereo audio source separation evaluation campaign: data, algorithms and results[END_REF]. Also, it was recently tested in the Signal Separa tion Evaluation Campaign (SiSEC) [10] where results showed that our algorithm compares favorably with the others pre sented in this campaign in most of the real and synthetic convolutive situation. As these results can be accessed at [10], only the results for the first audio source separation campaign will be presented. Signals are 16 kHz-sampled and have a lOs duration. Four types of mixtures were generated: 4 female speakers, 4 male speakers, 3 musics with one of them is drums and 3 other musics (no drums). Two different sets of source positions are used, one for speech and another for music mixtures. Angles and distances of these positions are given in Table 2. Our algorithm uses Hanning windowed 2048 sample frames, and reconstruction is achieved using the overlap and add method with 75% overlap.

Results

Detailed performance are given in Table 1 whereas the mean of each case and the overall performance OVAP (mean per formance on all mixture types) are plotted in Figure 3. When looking at the SDR which computes the global separation performance, OVAP column, the proposed MBUSS gives the best results as compared to the others. It shows a slight ad vantage for the new MBUSS 3.61dB compared to the old one 3.5dB and a more efficient as compared to KMeans 3.42dB. Furthermore, the intermediate errors show that the compro mise operated by the algorithms is different.

Comparing the new MBUSS to the old one, results show that the proposed algorithm gives better results in terms of SAR (SAR=8.05 and 7.19 respectively) and, almost, the same in terms of SIR (SIR=7.04 and 7.27 respectively), which justifies the use of the theoretical distribution and not the approximated one as in [START_REF] Ei-Chami | A new-model based underdetermined source separa tion[END_REF]. Comparing the new MBUSS with the KMeans, Fig. 3 (a),(b) and (c) show that the MBUSS method performs better for speech mixture, whereas KMeans performs better on music ones. One ex planation could lie in the nature of the signal themselves: music sources are very resonant, hence the energy is highly concentrated on some (usually harmonic) frequencies, which ensures a quasi-disjoint spectro-temporal supports between instruments. Thus, hard masking would be more suitable fil ter for this type of mixtures than a soft filter. The OVAP on Fig. 3 (d) also shows that the KMeans algorithm favors weak interferences (SIR=8.34) to the cost of more degraded sep arated speech (SAR=7.22), in contrast to our method that presents stronger interference (SIR=7.04) but with a clearer separated speech signal (SAR=8.05).

These results are confirmed by informal listening tests: less audible artifact are noticed in the proposed MBUSS, at the expense of more interference, which reveals nicer to lis ten to since our auditory system is very sensitive to artifact (gurgling noise) and less to interference. One explanation (partially) lies in the nature of the mask: as compared to the binary mask used in the KMeans method, the soft mask slightly smooths the output, limiting the isolated errors (the artifact like musical noise), but favoring the presence of inter ference.Note that the mask is used " as is": we are inclined to think that some "wise" smoothing (linear or nonlinear) could give better objective and subjective results, reducing even more the artifact.

CONCLUSION

In this paper, a new method to solve the underdetermined blind separation of audio mixtures problem has been pre sented. It is based on the sparseness assumption and on a theoretical model for the interchannel cues (log (ILD) ,IPD) given in [START_REF] Pham | Modeling the short time fourier transform ratio and ap plication to underdetermined audio source separation[END_REF]. By an EM algorithm, we were able to estimate the parameters of this model and then build the time fre quency separating masks. The algorithm demonstrated its

I 1 '

 1 -, -, IR(t)1 --IA ,JA,ql � 0 .• 1 .

Figure 1 :

 1 Figure 1: Ideal and observed (IPD,log(ILD)) for a single source mixture at the 100Hz center frequency

Figure 2 :

 2 Figure 2: Basic scheme of Model Based Source Separation Approach

L

  L 7r i (t) log MiP p ; log ----:;:

3. 3

 3 Relaxing the M-step(10)

3. 3 . 1

 31 Basic stepUsing Jensen ' s inequality, it can be easily proved that C (r�, p�) -C (ri ' Pi) has the following lower bound:-2( log L Wi(t) {cosh[�i(t) -Ai] -P; cos [<Mt) -(;li l } (11) IRi (t) I] -P cos [arg Ri(t)] 2IRi(t)l7ri(t) IRi(t)12 + 1 -2pR[Ri(t)]

5 Table 2 :

 52 Source positions from the set of microphone 4. 1 Comparison Algorithm and performance mea surementThe proposed algorithm (referred to as new MBUS) is com pared with two other algorithms: the first one (referred to as old MBUS) is our previous method based on the same model-based EM approach, but with a simplifying probabil ity distribution of the couple of variable (IPD ,log(ILD))[START_REF] Ei-Chami | A new-model based underdetermined source separa tion[END_REF]. The second one (referred to as Kmeans) is the KMeans un derdetermined source separation presented in[START_REF] Araki | K-means based underdetermined blind speech separation[END_REF], in which the ILD /IPD observations are grouped into N clusters using a KMeans algorithm: each cluster center gives an estimation of the mixing matrix and the cluster point sets give the bi nary TF separating masks. The separation performance was evaluated for each estimated source i by the same criteria used in SiSEC: Signal to Interference Ratio (SIRi), Image to Signal Ratio (ISRi), Signal to Distortion Ratio (SDRi) and Signal to Artifact Ratio (SA�). For a detail description of these criteria and of their computation, the reader may refer to [10].

6Figure 3 :

 3 Figure 3: (a),(b),(c) presents respectively the SDR, SAR and SIR performance for each type of mixture; (d) presents the OVAP

Table 1 :

 1 MeEIl OVAP SOR 3,91 2,68 3,71 4,17 3,61 4,46 3,40 4,37 5,27 4,37 2,42 5,80 4,81 4,34 1,65 3,49 0,35 1,83 3,61 neN ISR6,17 5,99 6,52 6,87 6,39 6,58 7,28 7,07 7,52 7,11 4,27 8,61 8,61 7,16 1,76 6,38 0,40 2,85 6,00 MBUSS SIR 5,75 2,62 5,37 6,89 5,16 6,75 3,45 6,77 10,13 6,78 4,07 8,80 6,57 6,48 13,38 6,15 11,88 10,47 7,04 SA.R 6,21 5,52 6,19 6,08 6,00 7,78 7,44 6,82 7,72 7,44 5,50 11,24 8,96 8,56 15,90 9,73 7,56 11,06 8,05old SOR 3,25 1,97 3,69 4,11 3,25 4,43 3,84 4,44 5,50 4,55 1,69 5,87 5,19 4,25 1,64 3,02 0,34 1,67 3,50 Results for synthetic recording with a 5 cm microphone spacing and two different types of sources: speech and music. the overall performance OVAP of the source separation is presented in the last column. All figures are given in dBs

	4 Men Speedl 52 53 54 Mean S1 53 MBUSS 4 Women Speedl 3 Music, No drums 3 Music, with drums S1 52 53 54 Mean S1 52 53 MeEIl S1 52 ISR 6,14 4,17 6,87 7,21 6,10 6,34 8,03 7,89 8,65 7,73 4,31 8,30 10,50 7,70 1,80 7,16 0,42 3,12 [81 SIR 3,73 2,21 5,20 6,27 4,35 8,45 5,35 7,24 9,53 7,64 3,15 9,78 7,15 6,69 14,85 6,94 11,94 11,24 7,27 6,27 SA.R 5,99 5,04 6,25 6,17 5,86 6,78 6,95 6,68 7,29 6,93 5,60 11,26 9,00 8,62 11,59 6,09 6,01 7,89 7,19 SOR 3,50 2,24 3,16 3,57 3,12 4,03 3,10 3,77 5,22 4,03 1,82 5,58 5,15 4,19 1,72 4,62 0,36 2,23 3,42 KmeEIls ISR 7,64 5,90 7,38 8,65 7,39 6,16 10,55 9,20 10,06 8,99 5,29 8,00 10,28 7,86 1,82 8,75 0,42 3,66 7,15 [31 SIR 7,01 4,62 6,36 7,08 6,27 12,51 5,52 8,04 9,81 8,97 4,13 8,49 7,05 6,55 14,73 7,77 13,68 12,06 8,34 SA.R 4,58 2,40 4,24 4,54 3,94 4,69 5,65 5,53 6,31 5,54 7,34 11,27 9,68 9,43 16,26 10,92 7,75 11,64 7,22

ai, Ai, p;. By developing the numerator terms cosh[�i(t) -Ai] and COS[¢i(t) -ai], the maximum is obtained when ' 2 Pi -1 I:t Wi(t) sinh �i(t) tanh I:t Wi (t) cosh �i(t) arg L wi(t)sign {Ri(t)} [I:t Wi (t) COS¢i(tW + [I:twi(t)sin¢i(tW [I:t Wi (t) cosh �i (t) J2 -[I:t Wi (t) sinh �i (t) J2


domain: spi phys

text:
Michael Niestroy

Robert Luckner

Carsten Döll

Flight Control Systems for Aircraft Engaged in Air Wake Surfing for Efficiency

This paper gives an overview of a flight control architecture developed during studies for a trail aircraft flying in formation with a leading aircraft in order to save fuel, commonly referred to as air wake surfing for efficiency (AWSE). The paper starts with an overview of the motivation for formation flight followed by a generic flight control architecture and accompanying algorithms that achieve formation flight. Finally, a summary of some of the results from simulation of a formation of turboprop transport aircraft and formations of a tanker or high-speed transport aircraft with fourth generation fighters is presented.

I. Nomenclature

II. Introduction

In order to reduce the fuel consumption of an aircraft, fuel-saving formation flight has been studied for several decades. An overview about already conducted investigations including flight tests is for example given in [START_REF] Vechtel | Flight dynamics simulation of formation flight for energy saving using LES-generated wake flow fields[END_REF]. The first manual flight tests have been conducted in 1984 by the Technical University of Braunschweig, Germany [START_REF] Beukenberg | Flugversuche zur Messung der Leistungsersparnis im Verbandsflug[END_REF], the first ones with a modified autopilot in 1987/88 also by the Technical University of Braunschweig [START_REF] Beukenberg | Aerodynamics, performance and control of airplanes in formation flight[END_REF]. The basic principle is that a trailing aircraft flies within the outer section of the wake vortex flow field behind a leading aircraft where an updraft can be encountered like shown on Figure 1. This principle is commonly called "wake surfing". The trailing aircraft must be flown within a very narrow section of the flow field where the upwind velocity leads to the best fuel savings. This position within the flow field is commonly called "sweet spot" and is typically slightly outside of the wake vortex core requiring some wingtip overlap of the two aircraft (in case that both aircraft are of the same type). Applying classical trim computations yields the sweet-spot location as shown on Figure 2 for an A320 following another A320. The corresponding rolling moment induced by the flow field is plotted on Figure 3. Fuel savings of more than 10% can be achieved in a practical area (deep blue) around the sweet spot. It is noticeable that the rolling moment is close to zero at the sweet spot but with a large gradient along the lateral direction.

The white area on the left side of these figures represents the region of the flow field where the induced roll moment and the downdraft cannot be countered by the ailerons and the elevators of the trailing aircraft. Therefore, a flight into this hazardous section of the wake vortex should be avoided.

The flight tests have illustrated that it is very tedious for pilots to track the sweet spot manually for a long time without encountering the hazardous area of the wake vortex. Therefore, an automatic flight control system on the trailing aircraft is needed for safe and effective operation. The next sections describe a candidate control architecture and present some simulation results.

½ -1 nautical mile separation

III. Wake Surfing Flight Control System Architecture

The wake surfing flight control system development began by researching existing formation control architectures and verifying/judging whether or not they made sense for complex nonlinear simulations that have been created for wake surfing, see for example [START_REF] Yadlin | A method and applications for tracking airplane trailing wakes[END_REF], [START_REF] Kaden | Modeling wake vortex roll-up and vortex-induced forces and moments for tight formation flight[END_REF], [START_REF] Misaka | Large-Eddy Simulation of aircraft wake evolution from roll-up until vortex decay[END_REF].

Much of the current work in the field of wake surfing flight control systems is for "ideal" situations and does not take into account the nonlinear disturbances within the flow field. Numerous papers also ignored noisy position signals, signal lag, and atmospheric turbulence. Very simple aircraft models were used in some papers, in other cases the control algorithms were only demonstrated for two-dimensional tracking. It was necessary to test many different control algorithms to determine which control scheme or variation of control scheme offered the best tracking and robustness in the presence of the multiple realistic characteristics that were modeled in beforementioned simulations considering position sensor noise, signal lag (especially concerning ADS-B signal), and light turbulence.

The flight control system in total is either developed for the wake surfing (as in this paper) or consists of some add-on functions coupled to the outer loop autopilots that were previously designed around a classical inner loop control system, as shown in references [START_REF] Hanson | An overview of flight test results for a formation flight autopilot[END_REF], [START_REF] Kienitz | Modeling and design considerations for robust control of aircraft in tight cruise formation[END_REF], [START_REF] Koloschin | Towards automatic control of fuel-saving formation flight[END_REF], and [START_REF] Beukenberg | Aerodynamics, performance and control of airplanes in formation flight[END_REF].

A successful wake surfing flight control architecture is discussed in the following sections. It consists of an inner loop control that uses dynamic inversion to achieve the desired dynamics of the aircraft. The outer loop formation flight control is a blend of various algorithms to achieve both zero steady-state error tracking and a reasonably rapid response to disturbances.

The control architecture and algorithms that follow were created for the generic turboprop transport 6DoF model GTTM [START_REF]Revolutionary Configurations for Energy Efficiency (RCEE) Final Report Contract[END_REF] but can be applied to other aircraft by changing the desired dynamics, the outer loop gains, and the thrust controller.

A. Dynamic Inversion Inner Loop Control

An inner loop controller was designed using dynamic inversion with desired dynamics is shown in Figure 4. Desired roll rate and pitch rate commands from the outer loop controller are used to calculate desired roll rate, pitch rate, and yaw rate accelerations for a given set of desired dynamics (where the actuator dynamics Act were also included). The effector blender calculates the aileron, elevator, and rudder commands from the difference of the desired angular rates and the current angular rates.

The desired roll rate acceleration ̇ was tuned by adjusting the roll mode time constant, Roll, as shown in Equation (1) to meet the handling qualities of the class of aircraft.

̇ = + ( -) (1)
The desired pitch rate acceleration ̇, as seen in Equation ( 2), was tuned by adjusting the short period natural frequency sp, its damping ratio sp, and the desired pitch time constant Tθ2 to meet the longitudinal handling qualities of the class of aircraft. The desired yaw rate acceleration ̇, Equation (3), was tuned by adjusting the Dutch roll mode natural frequency dr, its damping ratio dr, and the desired yaw rate time constant TR1 to meet the appropriate handling qualities metrics.

̇ = + ( -) - (2)
̇ = + 2 ( + ) + - 2 - 1 1 (3)
As shown in Equations ( 4) and ( 5), the thrust inner loop control calculates a desired thrust value Thrust based on the outer loop commanded body x-axis acceleration ̇and the current flight condition. This desired thrust value is used in a reverse lookup on the engine table to find the required throttle input TLACmd for the desired thrust and flight condition. The vertical position outer loop control, as seen on Figure 5 and defined in Equation ( 6), is maintained by generating pitch rate commands QCmd that are fed to the inner loop controller. The pitch rate command is generated by summing the vertical position error Zerr, the derivative of the vertical position error ̇ , and the double derivative of the vertical position error ̈ with gains Kp, Kd, Kdd applied to each value. An iterative approach was used to tune these gains. It was found that the double derivative term is the largest contributor to maintaining stability. The same control algorithms were used on both the leading and trailing aircraft. However, the gains on the trail aircraft were modified to allow faster response and better tracking. A slower response was desired for the leading aircraft since high frequency motion on the part of the lead can drive the trail aircraft to become unstable.

ℎ = cos 2 57.3 ̇ -( - -sin ) - ( cos -sin ) (4) = ( ℎ , ℎ, ) - (5)
The vertical outer loop controller drives the vertical tracking error to zero by commanding a pitch rate QCmd.

= -

+ ̇ + ̈ (6)

C. Forward Outer Loop Control

Figure 6 Forward Outer Loop Controller.

The forward position outer loop control, as seen on Figure 6 and Equation ( 7), was done in a manner similar to the vertical controller. A desired body X acceleration command, ̇ , was generated by summing the forward position error Xerr and the derivative of the forward position ̇ with appropriate gains Kp and Kd applied to those values. Since there is not a forward motion control mode for desired dynamics, the dynamics of the engine with a thrust response transfer function was included in the acceleration command ̇

. For this aircraft, the same control algorithm was used for both the leader and followers; however, for other aircraft a velocity hold algorithm was used. Since either algorithm produces adequate performance, it is left to the designer to determine what works best for them given their aircraft dynamics.

The outer loop forward control maintains the forward position relative to the lead aircraft by commanding a forward acceleration ̇ . The lateral position tracking control, as seen on Figure 7 and Equations ( 8), ( 9), and ( 10) was the most difficult task during the control design since there are strong lateral-directional non-linearities in the flow field right at and near the sweet spot (see Figure 3) that the controller is trying to maintain. For this reason, extra control augmentation was necessary to achieve and maintain stable tracking in formation. The roll control command was first generated using the same method of summing the lateral position error Yerr, with the lateral position error derivative ̇ , and the lateral position error double derivative ̈ . Lag filters were added to the derivative signals to help limit the noise propagation and amplification. A lead filter was added to the controller to increase the error signal when the aircraft moved rapidly in the lateral direction, which helped to boost the control signal that drove the aircraft back to a stable position.

Open loop control that was dependent on the lateral movement in the flow field was added to the control command .

Basically, the open loop control command determined the incremental change in roll rate due to lateral motion in the flow field, CRdy. This incremental roll rate value was then subtracted from the roll command that was generated. The removal of the roll command allowed the controller to use the rolling moment of the flow field to its advantage to maintain lateral control of the aircraft.

The outer loop lateral position control included a controller similar to the vertical controller , but also included extra roll command augmentation to account for the strong roll nonlinearities in the flow field. The outer loop lateral control maintains the lateral position relative to the sweet spot by commanding a roll rate .

= 1 + . + 1 + . + 1 + . (8)
=

= - [START_REF]Revolutionary Configurations for Energy Efficiency (RCEE) Final Report Contract[END_REF] ., An important part of the formation tracking task is applying adequate filters, as seen on Figure 8, to the GPS position signals since there will be some amount of noise in the signal. Since the derivatives of these signals are taken, it is especially important to filter the GPS position inputs. Second order transfer functions were used to filter the noisy position signals. The appropriate frequency and damping were found through iteration. A more advanced signal filtering method may be desired for future work.

The ADS-B position information from the leading aircraft can be received with significant delays or and may contain data dropouts. This also must be taken into account in the input signal filtering block of the overall control scheme.

IV. Results

Formation position controllers were developed for the generic turboprop transport model similar to a C130, generic high-speed mobility aircraft model similar to a C5, and a generic fourth generation fighter model similar to an F16. This section summarizes some of the results from simulations of the controllers implemented in MATLAB/Simulink for the trail aircraft along with lead aircraft dynamics and the flowfield disturbances caused by the lead aircraft wingtip vortices.

Table 1 shows simulation results of the turboprop transport aircraft alone and in various formations. The reference to system fuel saved or system lift/drag indicates those are the respective values of the aircraft taken as a group, so the fuel burn and savings includes that expended by the lead aircraft and all the trail aircraft. The ideal numbers are produced by looking at the theoretical savings when positioning the trail aircraft at the sweet spot. The simulated columns are results of the Simulink simulation over a period of time, including the effects of light turbulence on all the trail aircraft. Table 2 shows the results of simulations of the generic fourth generation fighter in formation with a high-speed mobility aircraft and generic tanker aircraft. Note the high-speed transport produces very strong vortices that can result in over 35% fuel savings on the fighter when flying at the theoretical sweet spot. With the vortex flow model for the respective lead aircraft and the type of controller described earlier but with the appropriate gains, simulation results show a fuel savings of just over 30% when in formation with the transport and approximately 12% when in formation with the tanker. Since the fighter is significantly smaller than the leaders and consume much less fuel, the fuel savings in this table are for the fighters themselves, not including the transport or tanker leads.

V. Conclusions and Perspectives

The benefits of wake surfing for fuel efficient operations was introduced and reasons for the need of an automatic flight control system on the trailing aircraft were highlighted. A candidate control architecture with an inner loop and complementary outer control loops was illustrated. Simulation results for various combinations of turboprop transports, high speed mobility aircraft, tankers, and fighters were presented. The results show that the candidate controllers are capable of accurate relative positioning, including sensor noise and delay, atmospheric turbulence, vortex wake, and aircraft dynamics.

=

  Aircraft angle of attack, rad = Aircraft sideslip angle, rad = Aircraft pitch angle, rad Roll = Roll mode time constant, sec Act = Actuator natural frequency, rad/sec sp = Short period natural frequency, rad/sec dr = Dutch roll natural frequency, rad/sec dr = Dutch roll damping ratio sp = Short period damping ratio Alt = Aircraft altitude, ft = Body y-axis acceleration, ft/sec 2 = Reference wing span, ft = Coefficient of lift = Coefficient of drag CRdy = Incremental change in roll rate due to lateral motion in the flow field, rad/sec g = Acceleration due to gravity, ft/sec 2 = Body axis roll moment of inertia, slug-ft 2 Kp = Controller gain on the error signal, proportional gain Kd = Controller gain on the first derivative of the error signal, derivative gain Kdd = Controller gain on the second derivative of the error signal, double derivative gain 1 LM Fellow, Control Sciences, Advanced Development Programs, Lockheed Martin Aero, AIAA Associate Fellow 2 Professor for Flight Mechanics and Flight Control, Aerospace Department -Flight Mechanics, Flight Control and Aeroelasticity. 3 Researcher in Flight Control and Handling Qualities, Information Processing and Systems Department. Ki = Controller gain on the integral of the error signal, integral gain Kff = Controller gain on the flow field = Aircraft mass, slug ̇= Angular roll acceleration, rad/sec 2 ̇= Angular pitch acceleration, rad/sec 2 ̇= Angular yaw acceleration, rad/sec 2 = Commanded roll rate, rad/sec = Roll rate command due to lateral offset error, rad/sec = Roll rate command augmentation, rad/sec ̇ = Desired roll acceleration, rad/sec 2 = Commanded pitch rate, rad/sec ̇ = Desired pitch acceleration, rad/sec 2 ̇ = Desired yaw acceleration, rad/sec 2 = Dynamic pressure, lb/ft 2 p = Roll rate, rad/sec q = Pitch rate, rad/sec r = Yaw rate, rad/sec = Reference wing area, ft 2Tθ2 = Pitch attitude numerator zero, sec TR1 = Yaw rate time constant, sec = Throttle lever angle command ℎ = Thrust desired to maintain longitudinal separation, lbs = Throttle lever angle at trimmed flight u = Body x-axis velocity component, ft/sec v = Body y-axis velocity component, ft/sec w = Body z-axis velocity component, ft/sec ̇= X-axis linear acceleration, ft/sec 2 ̇ = Commanded x-axis acceleration, ft/sec 2 = Aircraft true velocity, ft/sec = Error between actual follower location and desired distance from lead aircraft, ft ̇ = Rate of change error between actual follower location and desired distance from lead aircraft, ft/sec Yerr = Lateral position error, ft ̇ = Lateral position error first derivative, ft/sec ̈ = Lateral position error second derivative, ft/sec 2 Zerr = Vertical position error, ft ̇ = Vertical position error first derivative, ft/sec ̈ = Vertical position error second derivative, ft/sec 2

Figure 1

 1 Figure 1 Fuel saving formation flight is achieved by the trail aircraft flying within the updraft section of the lead aircraft wingtip vortex, referred to as wake surfing [1].

Figure 2

 2 Figure 2 Possible fuel savings within the flow field for an A320 following another A320 [4].

Figure 3

 3 Figure 3 Aileron deflection is necessary to compensate for the rolling moment induced by the flow field [4].

Figure 4

 4 Figure 4 Dynamic Inversion Inner Loop Controller.

,Figure 5

 5 Figure 5 Vertical Outer Loop Controller.

Figure 7

 7 Figure 7 Lateral Outer Loop Controller.

Figure 8 2

 8 Figure 8 2 nd Order Input Signal Filters.

Table 1 Summary of simulation results for varying numbers of C130 aircraft in fuel efficient formation flight.

 1

Table 2 Summary of simulation results for varying numbers of fourth generation fighter aircraft in fuel efficient formation flight with transport and tanker aircraft.

 2

Acknowledgments

This article has been prepared in the framework of the NATO CSO group Applied Vehicle Technologies AVT-279 "Formation Flying for Efficient Operations". The work has been financially supported by the affiliation organizations of the authors. Complete results and further documentation can be found in the AVT-279 final report available from the NATO CSO organization.