Marais. J Kazim S Ali

Aït Tmazirte

N Cosys-Leost

Syed Ali Kazim

Dr Juliette Marais

Interferences in Safety Critical Land Transport Application: Notch Filtering vs Wavelet Transform, an Experimental Analysis

a master's degree in Electronic Engineering and later he completed a specialization degree in

INTRODUCTION

The localization function designed for terrestrial applications inherits high-level performance requirements to guarantee accuracy, availability, continuity and operation safety. These requirements are perceived to be much more stringent for forthcoming generation as it aims to bring more intelligence, autonomy, connectivity and shared mobility services. It necessitates reliable positioning and navigation among other key requirements to accomplish some of these goals. It implies considering and managing all the risks and threats associated with a positioning system to ensure a fail-safe positioning solution, particularly in the case of satellite-based positioning systems. Global Navigation Satellite System (GNSS) is recognized as a Game Changer technology for many applications that can bring economic and ecological incentives, particularly in the railway sector for signaling applications. However, its performance is mainly associated to receiver operating conditions. This makes it very challenging for a satellite-based positioning system to meet performance requirements as the satellite signals are vulnerable to obstacles present around the receiver. It is particularly prone to unintended and malicious radio transmissions from jammers. There exist several approaches for the detection and mitigation of interference signals which are mainly implemented at the Front-end, pre-correlation, post-correlation, and navigation levels in the receiver chain [START_REF] Dovis | GNSS Interference Threats and Countermeasures[END_REF]. It includes signal processing techniques [START_REF] Anyaegbu | An integrated pulsed interference mitigation for GNSS receivers[END_REF][START_REF] Borio | Loop analysis of adaptive notch filters[END_REF]Borio, et al., 2008;[START_REF] Borio | Optimal global navigation satellite system pulse blanking in the presence of signal quantisation[END_REF][START_REF] Dovis | Performance comparison of transformed-domain techniques for pulsed interference mitigation[END_REF] to filter the interference signal that filters IQ samples either in time, frequency, or time-frequency/scale domain; Spatial filtering [START_REF] Fu | Suppression of multipath and jamming signals by digital beamforming for GPS/Galileo applications[END_REF][START_REF] Gupta | Space-frequency adaptive processing (SFAP) for RFI mitigation in spread spectrum receivers[END_REF] which requires processing signals acquired from several antennas mounted in an array. This permits to digitally steer the beam pattern in the direction of the satellite signals while suppressing the unwanted signal from the interference source; Vector tracking [START_REF] Lashley | Performance analysis of vector tracking algorithms for weak GPS signals in high dynamics[END_REF][START_REF] Nunes | Low-complexity VDLL receiver for multi-GNSS constellations[END_REF] allows to track the satellite signal in cooperation with other satellite channels. In this manner, it enhances the tracking capability of largely degraded signals consequently increasing the overall performance. Sensor fusion [START_REF] Groves | Combating GNSS interference with advanced inertial integration[END_REF]) makes use of the information from the inertial sensors to counter the impact of interference.

In this study, we are interested in the detection and mitigation of interference signals at the pre-correlation level. The primary goal is to evaluate performance from the user's perspective. Several countermeasures at the pre-correlation level show high performance in mitigating the effect of an interfering signal. However, their performance is assessed mainly at the signal level (at either acquisition or tracking stage). Recently, some studies presented the impact of interference and mitigation techniques effectiveness at the position level (Borio & Gioia, 2021;Gioia & Borio, 2021). Similarly, we also addressed interference mitigation from a safety perspective [START_REF] Kazim | On the impact of temporal variation on GNSS position error models[END_REF]. We provided KPIs assessment in the presence of interfering signal and after enabling mitigation block implemented within the receiver. This study is the extension of our previous work where we now have implemented two interference mitigation techniques namely adaptive notch filter (ANF) and wavelet packet decomposition (WPD) to counter the interference signal. For each technique, it presents an in-depth discussion and a comparative analysis at multiple levels (frequency estimation, acquisition, tracking, and position). Furthermore, it also provides an assessment at the positioning level discussing previously mentioned KPIs using the Stanford representation. This study provides a performance assessment for two different types of time-varying CW interference signals namely frequency hopping and chirp signal The paper is organized into several sections. Section 2 provides some details on the state of the art interference mitigation techniques namely Adaptive Notch Filter (ANF) and Wavelet Packet Decomposition (WPD). Section 3 provides an overview of the used positioning algorithm based on the weighted least square (WLS) solution. Three well-known weighting schemes based on elevation, carrier to noise ratio (C/N0) and a hybrid model are briefly presented. Moreover, it also presents a classical manner to estimate the protection level using Hslope estimation. Section 4 describes the experimental setup used for data recording and it also provides details on the interference scenarios and acquisition campaigns performed during the day at regular intervals. Finally, a performance assessment based on the experimental results for the two cases: interference without any mitigation and after applying interference countermeasures for the chirp and frequency-hopping scenario is presented before discussion.

INTERFERENCE SUPPRESSION AT THE PRECORRELATION STAGE

In literature, the issue of radio frequency interference is largely investigated at the pre-correlation level. Several techniques have been proposed for the detection, characterization, mitigation and in some cases localization of the interference source. The most common mitigation method includes Pulse Blanking (PB), Adaptive Notch Filter (ANF), Short Time Fourier Transform (STFT), Wigner-Ville Distribution (WVD), Wavelet Transform (WT), and Kerhunen-Loève Transform (KLT). In general, each technique aims to estimate the interference components and subsequently applies a mitigation process for the suppression. In this work, Adaptive Notch Filter (ANF) and Wavelet Packet Decomposition (WPD) have been investigated as interference mitigation techniques. This section will provide a brief overview of the two considered mitigation techniques.

Adaptive Notch Filter (ANF)

An Adaptive Notch Filter (ANF) is a widely used method for the suppression of narrowband interference signals.

It is comprised of a notch filter to reject a very narrow portion of the signal and an adaptive unit to continuously estimate the instantaneous frequency of the jamming signal. Thus, in this manner, it effectively tracks and suppresses the interfering signal without affecting much of the useful component of the signal. The transfer function of the singlepole Infinite Impulse Response (IIR) notch filter is given as:

𝐻(𝑧) = 1 -𝑧 0 𝑧 -1 1 -𝑘 𝛼 𝑧 0 𝑧 -1 (1)
where 𝑘 𝛼 𝜖 [0,1) represents the pole contraction that controls the bandwidth of the notch, 𝑧 0 represents filter zero in the complex plane which in the lock state corresponds to the interference frequency. The transfer function of the filter can be split into the moving average (MA) and autoregressive (AR) blocks. The output from each block can be expressed as

𝑥 𝑟 [𝑛] = 𝑥[𝑛] + 𝑘 𝛼 𝑧 0 𝑥 𝑟 [𝑛 -1] (2) 𝑦 𝑚 [𝑛] = 𝑥[𝑛] -𝑧 0 𝑥 𝑟 [𝑛 -1] (3)
Where 𝑥 𝑟 [𝑛] and 𝑦 𝑚 [𝑛] represents the output from the AR and the MA function respectively. Ideally, a very narrow rejection bandwidth allows the filter to largely suppress interference frequency while preserving mostly useful content of the signal. This condition is obtained when the notch frequency is equivalent to the interference instantaneous frequency. The 3dB bandwidth of the filter can be estimated as:

𝐵 3𝑑𝐵 ≈ (1 -𝑘 𝛼 )𝑓 𝑠 𝜋 10 ⁄ (4)
where a narrow notch is achieved as the value of 𝑘 𝛼 approaches to unity (𝑘 𝛼 → 1). The function of the adaptive block is to move filter zero 𝑧 0 progressively in a complex plane until it converges to interference frequency. It is based on the iterative normalized least mean square (LMS) algorithm, which minimizes the cost function. It could be estimated either by minimizing the expectation of the filter output energy or by minimizing the expectation of MA output energy.

The minimization is achieved from the filter output through an iterative rule to estimate filter zero 𝑧 0 in the runtime.

𝑧 0 [𝑛] = 𝑧 0 [𝑛 -1] -𝜇[𝑛]𝑔( 𝐽[𝑛]) (5)
where An adaptive notch filter serves as a very simple but very effective technique to suppress interference, particularly for chirp signals. However, its performance depends on the appropriate selection of the filter parameters [START_REF] Qin | An Assessment of Impact of Adaptive Notch Filters for Interference Removal on the Signal Processing Stages of a GNSS Receiver[END_REF].

The two important parameters are notch width 𝑘 𝛼 and algorithm step 𝛿 which must be tuned properly to achieve better performance. Furthermore, the main drawback of ANF is that it is not a very effective method for dealing with multiple interfering signals simultaneously present in the GNSS band. [START_REF] Landry | Impact of interference on a generic GPS receiver and assessment of mitigation techniques[END_REF] presented an all-pass filter with a Gauss-Newton algorithm as an adaptive algorithm to enhance the convergence rate to suppress CWI and swept CW interference. (Borio, Camoriano, & Presti, 2008) presented multi-pole adaptive filter by cascading several two-pole IIR notch filters to progressively blocks multiple CW interfering signals.

Discrete Wavelet Transform (DWT)

Wavelet analysis is well-known among the families of transformation techniques used to characterize the signal components. In wavelet transform, a discrete signal is decomposed using a set of orthogonal functions which are derived by scaling and shifting the mother wavelet function. In this manner, it uses a scalable window to transform a time domain signal into a time-scale (frequency) domain. Thus, it overcomes the shortcoming of fixed time window operation as in the case of Short Time Fourier Transform (STFT).

In Discrete Wavelet Transform (DWT), the decomposition stages are mainly realized by bandpass filtering and subsampling process. The signal is passed through a bank of wavelet-based bandpass filters with different cut-off frequencies at different scales. The filtering operation splits the incoming signal into low and high frequency bands at each decomposition stage. The low pass filter, equivalent to averaging operation extracts the coarse approximation whereas the highpass filter similar to differencing operation acquires the detailed information of the signal. Later, a decimation process equivalent to a downsampling operation is performed to reduce sample redundancy induced by the filtering process. It is accomplished by eliminating every alternative sample to remain in compliance with Nyquist's sampling criterion. The same procedure is repeated at each level until the desired decomposition level is reached. Consequently, at every decomposition stage, the time resolution is reduced by the factor of ½ (half of the samples describe the decomposed signal) gaining twice in the frequency resolution. This reduces the uncertainty in the frequency as the signal splits into low and high frequency bands. The decomposition stages of a discrete signal 𝑥(𝑘) for different bands at different scales is shown in Figure 1.

The decomposition level can be mathematically expressed as:

𝑦 𝑙𝑜𝑤 [𝑛] = ∑ 𝑥[𝑘]𝑔[2𝑛 -𝑘] ∞ 𝑘=-∞ (7) 𝑦 ℎ𝑖𝑔ℎ [𝑛] = ∑ 𝑥[𝑘]ℎ[2𝑛 -𝑘] ∞ 𝑘=-∞ (8)
where 𝑦 𝑙𝑜𝑤 and 𝑦 ℎ𝑖𝑔ℎ represent the decomposed signal output from the low pass and high pass filters respectively and after performing downsampling by factor 2. The scaling vector 𝑔[𝑛] and wavelet vector ℎ[𝑛] represent impulse response of low pass and high pass filter respectively. Discrete Wavelet Transform (DWT) can be further extended to form a decomposition tree also known as Wavelet Packet Decomposition (WPD). The two methods differ in the sense that, in WPD, the detailed branch (high frequency band) is also decomposed similarly to the approximation branch (low frequency band) to form a tree-like structure. This process generates several additional wavelet packets after passing a signal through a uniform wavelet filter bank. Each decomposition stage forms 2 𝑁 wavelet packets representing some portions of the input signal where N represents the decomposition level. A three-level wavelet decomposition packet using uniform filter banks is shown in Figure 2.

The main purpose of applying signal decomposition is to break down the incoming signal to get a representation where interference content can be easily identified and isolated from the useful signal. Therefore, the interference mitigation process further requires detection, mitigation and reconstruction steps.

The detection phase, at first, requires calibration of each wavelet packet in the absence of interference. In this way, a threshold is approximated in a statistical manner represented by variations in the standard deviation in time-scale representation. An interference signal is considered to be present whenever the wavelet coefficients surpass the set threshold value previously defined in the calibration phase.

In the mitigation phase, all the wavelet coefficients that exceed the set threshold defined during the calibration process (in the absence of interference) are given a null value also known as coefficient blanking. Thus, the blanking method allows the removal of extra energy which is due to the presence of an interference signal.

In the reconstruction step, the signal is reconstructed using the inverse wavelet transform. This is done using the reconstruction filters and upsampling operation while back tracing all the steps done in the decomposition process.

The performance of wavelet techniques is well assessed in the presence of CW and pulse interference signals. In [START_REF] Anyaegbu | An integrated pulsed interference mitigation for GNSS receivers[END_REF], Wavelet Transform is used to counter the effects of pulsed interference. [START_REF] Seo | Performance of Interference Mitigation with Different Wavelets in Global Positioning Systems[END_REF] provided a comparative study using different wavelet filters for the removal of a narrowband signal. [START_REF] Musumeci | Use of the wavelet transform for interference detection and mitigation in global navigation satellite systems[END_REF], 2013) investigated a wavelet-based interference mitigation method for pulsed and narrowband signals and analyzed the effectiveness of mitigation in the acquisition and tracking stages. [START_REF] Musumeci | A Comparative Analysis of Adaptive Notch Filtering and Wavelet Mitigation against Jammers Interference[END_REF] extended further the analysis and evaluated mitigation performance at position level in the presence of a narrowband interference signal.

INTEGRATION ALGORITHM, WEIGHTING SCHEMES AND PROTECTION LEVEL COMPUTATION

The primary goal of this study is to analyze the performance of previously presented interference mitigation techniques from the user's point of view. For this purpose, we evaluate some important performance indicators such as availability, accuracy and safety at the position level. This section presents some details on the used positioning algorithm, weighting schemes and Horizontal Protection Level (HPL)

Weighted Least Square (WLS) Position Estimator

The most common method to compute the position from a set of pseudorange measurements is based on Least Square (LS) method. In LS, the user position is estimated recursively relatively to the linearization point by using pseudorange errors. The linearized pseudorange equation can be written as:

𝛿𝜌 = 𝐻 × 𝛿𝑋 + 𝜖 ( 9
)
where 𝛿𝜌 is the difference between the actual and the predicted pseudorange estimated from the linearization point, 𝐻 is the observation matrix representing satellite-receiver information, 𝛿𝑋 the error state vector, 𝜖 is the measurement error.

The Weighted Least Square (WLS) estimator can be expressed as:

𝛿𝑋 = [𝐻 𝑇 𝑊𝐻] -1 𝐻 𝑇 𝑊 × 𝛿𝜌 ( 10
)
where W is the observation weighting matrix given by

𝑊 = ( 1 𝜎 1 2 ⁄ ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 1 𝜎 𝑚 2 ⁄ ) (11)
The measurement weights are inversely proportional to the error variance (𝑤 𝑖 = 1 𝜎 𝑖 2 ⁄ ). In WLS, different weights are given to the measurements based on the measurement noise. Therefore, measurements with a small error variance carry more weight hence more confidence to contribute to the position estimation. On the other hand, measurements with large error variance would get lower weight as a result have either less or no contribution at all in the position estimation.

Weighting Schemes

In this discussion, we present some of the most common weighting models that include carrier to noise ratio (C/N0), satellite elevation and a hybrid model (combination of elevation and C/N0 model). In general, weighing models are somehow used to give confidence to the measurements.

The carrier to noise ratio (C/N0) is a measure of the power level of the received signal relative to the noise power level per unit bandwidth. It is the key parameter that defines the quality of the signal. [START_REF] Hartinger | Variances of GPS Phase Observations: The SIGMA-ɛ Model[END_REF] presented a generalized expression to estimate the variance from the C/N0 measurements. The expression for sigma-ε is expressed as:

𝜎 2 = 𝑎 + 𝑏. 10 -0.1 × 𝐶/𝑁 0 (12)
where, a (in m 2 ) and b (in m 2 Hz) are the model parameter that corresponds to receiver characteristics.

Satellite elevation can also be used to evaluate the quality of the signal. The elevation model exploits the fact that low elevation satellites are subjected to more errors than high elevation satellites. The variance for the elevation-based weighting model is given as:

𝜎 𝑖 2 = 1 𝑠𝑖𝑛(𝑒𝑙 𝑖 ) 2 ⁄ ( 12
)
The hybrid model takes into consideration of both elevation and carrier to noise ratio to estimate the variance. A hybrid model with a k factor as the LOS indicator [START_REF] Tay | Weighting models for GPS Pseudorange observations for land transportation in urban canyons[END_REF] increases the effect of deweighting when NLOS signal is detected and is given as:

𝜎 𝑖 2 = 𝑘 10 -𝐶/𝑁 0 𝑖 10 𝑠𝑖𝑛(𝑒𝑙 𝑖 ) 2 (13)

Horizontal Protection Level (HPL) Computation

The Horizontal Protection Level (HPL) is a statistical bound used to guarantee that the true position error is well confined. In general, the protection level is estimated by taking into consideration the impact of measurement errors and bias on the estimated position. Thus, it provides information to the user about confidence in the estimated position. It is mainly used to evaluate the reliability of the system. HPL can be traditionally expressed as the combination of a noise term 𝐻𝑃𝐿 𝑛 and a biased term 𝐻𝑃𝐿 𝑏 :

𝐻𝑃𝐿 = 𝐻𝑃𝐿 𝑛 + 𝐻𝑃𝐿 𝑏 (14)
The noise term is calculated according to the error propagation as for HPL computed with SBAS data (mostly for aeronautic applications):

𝐻𝑃𝐿 𝑛 = 𝐾(𝑃 𝑚𝑑 ) × 𝑑 𝑚𝑎𝑗𝑜𝑟 (15)
where 𝐾 is an inflation factor that permits to meet integrity risk requirements. It is generally obtained conservatively, by selecting the 4 degrees of freedom value in the χ 2 table; 𝑑 𝑚𝑎𝑗𝑜𝑟 is the error uncertainty along the semi-major axis of the error ellipse.

The bias term is the most affected by the change in the weighting strategy and the term HSlope represents the sensitivity of HPE to the bias of 𝑖 𝑡ℎ satellite.

𝐻𝑃𝐿 𝑏 = 𝑚𝑎𝑥 𝑖 (𝐻𝑆𝑙𝑜𝑝𝑒 𝑖 . 𝜎 𝑖 ) × 𝑝 𝑏𝑖𝑎𝑠 (16) where the slope is given as:

𝐻𝑆𝑙𝑜𝑝𝑒 𝑖 = √ (𝐻 𝑁

EXPERIMENTAL PROTOCOL 1. Data Acquisition

For this study, the Stella record and playback tool from M3Systems have been used to record raw IQ signals that use Universal Software Radio Peripheral (USRP 2954R) from National Instruments as the front end to perform initial filtering and downconversion. The recording system is connected to an antenna (Adsotech ASH11661) mounted on the rooftop of our building for relatively open-sky conditions and optimal satellite visibility. This also allows minimizing the signal disturbances induced by the presence of obstacles in the close vicinity of the antenna. As far as the jammer is concerned, we have used a customizable interference source to transmit the interference signal via USRP 2910.

The interference and GNSS signals are then combined through a radio frequency combiner and the interfered signal is fed to the signal recorder. The complete experimental setup is shown in Figure 3. We presented two interference scenarios: 1) frequency-hopping and 2) chirp signal. This allows us to analyze the robustness of each mitigation method against two interference signals considered. The spectrogram representing the time-frequency relationship of the two interference signals is shown in Figure 4. Table 1 provides some details of the interference used cases and the signal recording system. In total, we have recorded 18 files (9 for each interference scenario), each of 90 seconds length. The signal acquisition has been performed during 3 different periods of the day to get different satellite configurations and visibility. The interference signal of 30-second duration is injected after 55 seconds in each recording. The recorded files are later processed individually using a Matlab-based software receiver to conduct performance analysis.

The main purpose of this study is to investigate the impact of interference and the importance of mitigation from a safety perspective. Therefore, it provides performance assessment for some key positioning performance indicators such as accuracy, availability, continuity, and safety. A multi-objective comparative discussion has been made on the impact of interference and also on the performance of mitigation techniques at different stages of the receiver.

Results and Analysis

This section provides experimental results of a comparative performance analysis between a classical Adaptive Notch Filter (ANF) and a Wavelet Packet Decomposition (WPD) for the two considered mitigation techniques. It evaluates mitigation performance at multiple levels including interference identification and suppression and their impact on tracking and positioning level. Finally, it provides an analysis of positioning performance from the user perspectives which is not only limited to positioning accuracy but also discusses system availability and safety.

Identification and Mitigation of Interference

To suppress the effects of jamming, it is very important to understand the characteristics of the interfering signal. This allows selecting a countermeasure strategy in time, frequency, or time-frequency domains, where interference components can be easily identified to possibly isolate them from the useful signal. For this purpose, we applied ANF and WPD techniques to estimate the time-varying frequency of the interference signal. Figure 5. illustrates the estimated frequency by the adaptive notch filter in the presence of frequency hopping (left) and chirp (right) signal. As previously discussed, the effectiveness of ANF-based mitigation depends on the filter parameters that require appropriate tuning depending on the characteristics of the interference signal. For simplicity, we used the same configuration (𝑘 𝛼 = 0.8 𝑎𝑛𝑑 𝛿 = 0.05) of ANF for each of the two interference scenarios. Here, ANF appears to track very closely the interference frequencies in comparison to the scenarios illustrated in Figure 4. However, in the frequency hopping scenario, the filter has shown significant oscillations around the interference frequency resulting in higher estimation noise. It is important to note that no extensive analysis is conducted on filter convergence or frequency estimation at this level.

Figure 6 illustrates a decomposed form of the interfered signal represented in the time-scale domain after applying wavelet packet decomposition (WPD). Here, a 5-level decomposition is realized using the "Symlet" function providing 32 wavelet packets, each representing some portion of signal content. The detection floor (threshold mask) is represented by the black plane which separates the clean and the interfered signal component. It is estimated by individually calibrating each wavelet packet in the nominal conditions. From the time-scale representation, it is evident that interference signal components have occupied approximately 50% of the wavelet scales which shows a high signal activity compared to higher scales. After the detection process, interference mitigation is performed in which all the coefficient values exceeding the threshold mask are set to a null value, also known as coefficient blanking. The modified wavelet packets are then fed to the synthesis filter bank to achieve signal reconstruction.

Tracking

This section aims to assess the tracking capabilities of the receiver before and after applying the mitigation techniques. Figure 7 clearly shows the impact of interference on satellite tracking in presence of frequency hopping (left) and chirp signal (right). The effects of interference can be seen approximately after 58 seconds with the emergence of the interference signal. At this point, without mitigation, the receiver completely lost track of the satellite signal (red curve). Instead, it started to track the interference signal which overpowered the satellite signal. After the removal of interference, in each mitigation scenario, the receiver started to track back the satellite signal. In the case of the chirp signal (right), the two mitigation techniques show similar performance. However, in the frequency hopping case (left), WPD-based mitigation appears to perform better with lower tracking noise compared to ANF filtering (blue).

A similar trend appears in the carrier-to-noise ratio (C/N0) estimation representing the tracking quality of the signal. Figure 8 shows that the C/N0 value abnormally fluctuates in presence of an interference signal which after applying interference mitigation dropped a few dB from the nominal trend (before interference injection) due to the filtering effects and some residual interference. The two techniques show relatively similar performance in the removal of a chirp signal. However, in the frequency hopping scenario (left), WPD-based mitigation (green) appears to better recover the GPS signal resulting in higher C/N0 than ANF (blue). Therefore, WPD appears to be the most suitable method for the removal of frequency hopping signals.

Positioning

The final contribution of this work is to provide an assessment of the impact of interference and its mitigation at the position level particularly focusing on the safety concerns. For this purpose, elevation, carrier to noise ratio (C/N0) and hybrid (elevation and C/N0) weighting models are applied for position estimation and also to estimate the protection level. The system performance is analyzed using the Stanford representation providing a quick and clear view of the KPIs including accuracy, availability and safety.

The positioning points in the Stanford diagram are distributed into several zones, each describing the state of the localization function. The significance of each zone is attributed to the capability to bound the true positioning error while respecting the application requirements defined by the alarm limit (AL). The points are classified into normal operation (white), misleading information -MI (pink and orange), system unavailable (yellow), and hazardous misleading information -HMI (red). Note, an Alarm limit (AL) of 30m is purposely chosen for this analysis which could vary depending on the application requirements.

Firstly, we present a reference case to show the positioning performance in the nominal conditions. This includes all the instances (epochs) before interference injection. In this manner, it generalizes the achievable performance with a specific configuration of the software receiver.

In the absence of interference, no significant difference in the positioning accuracy is observed for each weighting model as seen in Figure 9. On average, the mean and the standard deviation of position error appear to be approximately 4.5m and 2.5m respectively. In a safety context, each model provides nearly 99% of normal operations ensuring the positioning integrity, however, with few misleading events. When C/N0 based model is applied, some points are declared unavailable by the integrity monitoring criterion. As expected, the hybrid model tends to be the most conservative with higher PL values.

Figure 10 represents positioning performance in the presence of interference with no mitigation applied for the chirp and the frequency hopping scenario. We can clearly see that the nominal operation reduced significantly no matter the weighting model implemented. This highlights the importance of a mitigation process for the removal of the interference signal which when not mitigated raises safety concerns for the localization system and consequently of the entire system. Interference clearly shows an adverse impact on positioning accuracy represented by a higher concentration of points in HMI (> 50 m error). Moreover, the C/N0 and hybrid model compared to the elevation-based weighting model in some way played a role in managing the effects of interference by increasing the protection level. This resulted in compensating some HMI with unavailability for system protection which is nevertheless enough to completely isolate positioning points affected by interference. additional undue unavailability (yellow) which is more the case in the hybrid model. The drop in C/N0 (as seen in Figure 8) appears to be a probable cause of such overprotection since each measurement experiences similar behavior thus increasing the overall protection level.

Figure 12 presents the positioning performance after applying ANF (upper panel) for suppressing the frequency hopping signal. It can be seen that the effectiveness of ANF has reduced significantly in the frequency case compared to the previously presented chirp signal. Interference mitigation has resulted in a significant decrease in HMIs events however there remain some points with more than 30 m of positioning error. The elevation-based model still appears to be the most unsafe with 865 MIs (pink) and 49 HMIs (red) whereas C/N0 and the hybrid model manage to remove the HMI events however resulting in higher system unavailability. It can also be seen that WPD-based mitigation outperformed ANF in the removal of the frequency hopping signal. It retrieved a similar performance as the reference scenario (as seen in Figure 9). It has significantly increased the nominal operations (white) and more importantly, it completely ensures system safety.

CONCLUSION

This paper reports an experimental analysis of interference mitigation with either ANF or WPD techniques. With this study, we contributed to address the impact of interference and the effectiveness of mitigation strategy from a safety (integrity) and user perspective. We presented three cases: 1) referencewithout interference, 2) in the presence of interference but without any mitigation technique, and 3) after applying mitigation technique. For comparison analysis, we prepared two interference scenarios; linear chirp and frequency hopping. Furthermore, classical weighting schemes based on elevation, carrier to noise ratio (C/N0) and hybrid (combination of elevation and C/N0) model are used to estimate the position and the protection level. The main conclusions can be summarized as follows. A detection and mitigation layer in the localization function is very essential to effectively remove the disturbances induced by the jamming signal. Indeed, as shown, an inadequate strategy can lead to hazardous consequences endangering the safety of the system eventually resulting in a fatal accident. However, interference mitigation at the signal level cannot guarantee to bring back the localization function to normal operation meeting the performance requirements. In most of the presented scenarios, mitigation significantly reduced HMIs but with increasing unavailability, as a tradeoff to compromise availability for the system safety. This could be improved by introducing a complementary source of positioning to keep the system in nominal conditions. The primary concern of an adaptive notch filter (ANF) is the appropriate selection of the filter parameters. It requires proper tuning that could vary depending on the characteristics of the jamming signal. Certainly, interference characterization could help in refining the choice of these parameters to improve mitigation performance and to retrieve the nominal situation. Additionally, the undue unavailability could also be reduced by the use of specific weighting calibration for the interference case.

WPD performed exceptionally well for the frequency hopping case as it nearly retrieved the same performance as in the nominal conditions. However, it could not provide similar performance for the removal of the chirp signal. Conceivably, changing wavelet function, filter length and (or) decomposition level could help in improving the performance against the chirp signal. This prospect needs to be investigated in future work.

Figure 1

 1 Figure 1 An Illustration of 3-Level Discrete Wavelet Decomposition (DWT); Approximation (A) And Detail (D) Block Represent the Low Pass and High Pass Filtering and The Subsampling Operation

Figure 2

 2 Figure 2 An Illustration of 3-Level Discrete Wavelet Packet Decomposition (WPD); Approximation (A) And Detail (D) Block Represent Low Pass and High Pass Filtering Respectively, and Subsampling Operation.

Figure 3

 3 Figure 3 GNSS Signal Recording System Connected to an Interference source (Left) and Antenna The Roof Mounted on The Roof of The Building to Capture Satellite Signals (Right).

Figure 4

 4 Figure 4 Interference Scenario with Frequency Hopping (Upper Panel) And Chirp Signal (Lower Panel) At Baseband.

Figure 5

 5 Figure 5 Interference Frequency Tracked by Adaptive Notch Filter (Adaptation Step =0.05 and Pole Contraction Factor =0.8) In Case of Frequency Hopping (Left) and Chirp Signal (Right).

Figure 6

 6 Figure 6 Time-Scale Representation with 5 Level Decomposition of GPS L1 Signal Interfered with Frequency Hopping (Left) and Chirp Signal (Right).

Figure 7

 7 Figure 7 Satellite (PRN 02) Tracking State in The Presence of Frequency Hopping (Left) and Chirp (Right) Interference; No Mitigation (Red), ANF (Blue) And WPD (Green).

Figure 8

 8 Figure 8 Carrier to Noise Ratio (C/N0) Estimated After Satellite Tracking Stage for PRN 02 In The Presence of Frequency Hopping (Left) and Chirp (Right) Interference; With No Mitigation (Red), After Applying ANF (Blue) And WPD (Green).

Figure 9

 9 Figure 9 Stanford Diagram Representing Performance in Nominal Condition Using Different Weighting Schemes; Elevation (Left), C/N0 (Middle) And Hybrid (Right).

Figure 11 Figure 11

 1111 Figure11represents the positioning performance achieved after applying ANF (upper panel) and WPD (lower panel) techniques for the removal of the chirp signal. It can be seen that each mitigation technique significantly improves normal operations, however, they showed similar performance in combating the chirp signal. The elevation model appears to maximize the normal operation (white) to more than 96% but some points remain unbounded resulting in more MIs (pink). On the contrary, the C/N0 and Hybrid model appears to be over-protective thus it brings some

Figure 12

 12 Figure 12 Stanford Diagram Representing Performance After Mitigation of Frequency Hopping Signal by ANF (Upper Panel) And WPD (Lower Panel), Using Weighting Schemes; Elevation (Left), C/N0 (Middle) And Hybrid (Right)

  with 𝐻 + = [𝐻 𝑇 𝑊𝐻] -1 𝐻 𝑇 𝑊 with indexes E and N representing the components alond the East and North axes in an ENU frame; 𝑆 = 𝐼 -𝐻𝐻 + ; 𝜎 𝑖 is the standard deviation of the 𝑖 𝑡ℎ measurement error; 𝑝 𝑏𝑖𝑎𝑠 = √𝑁𝑆𝑆𝐸 (Normalized Sum of Squared Error), 𝑝 𝑏𝑖𝑎𝑠 represents the bias in the space of test statistics. Finally, HPL is expressed as[START_REF] Walter | Weighted RAIM for Precision Approach[END_REF]:

	,𝑖 + ) 2 + (𝐻 𝐸,𝑖 + ) 2 𝑆 𝑖𝑖	(17)
	𝐻𝑃𝐿 = 𝐾(𝑃 𝑚𝑑 ) × 𝑑 𝑚𝑎𝑗𝑜𝑟 + 𝑚𝑎𝑥 𝑖 (𝐻𝑆𝑙𝑜𝑝𝑒 𝑖 . 𝜎 𝑖 ) × 𝑝 𝑏𝑖𝑎𝑠	(18)

Table 1

 1 Receiver Configuration Used to Record GPS L1 signal with Interference signal characteristics.

ACKNOWLEDGEMENTS

The present work benefits from the financial support of the ANR French Research Agency within the project LOCSP No 2019-CE22-0011

received a Ph.D. degree in electronics and the "


Beate Orberger

Serge Miska

Alina Tudryn Geops

Christiane Wagner

Michel Fialin

Omar Boudouma

Colette Derré

Ghasem Nabatian

Maryam Honarmand

Iman Monsef

Abdolreza Ghods

Iron-oxide mineralogy of banded iron formations in the Takab region, North Western Iran

Keywords: Banded iron formation, Precambrian, Iran, magnetite, maghemite, hematite, goethite

Siliceous banded and nodular iron ores intercalated in mica schists and marbles of Precambrian age, were studied from the NW part of Iran (Takab region). The iron ore is composed of magnetite, which was partly transformed into maghemite and hematite along fractures and cleavages during a major oxidizing deformation event. Quartz shows undulose extinction. Both quartz and iron oxides partly form 120° triple junctions, but essentially grain boundary migrations. In the nodular ore, grain boundary migration is more pronounced leading also to grain isolation quartz is partially annealed. Iron oxides are agglomerated, elongated and disrupted. These finding indicate a longer recrystallization episode in the nodular ore. The iron oxide grains contain inclusions of euhedral, zoned and corroded zircons and euhedral monazite. Magnetite contains traces of Al. Hematite and goethite show higher Al, Si and Ca contents. Later hydrothermal solutions precipitated goethite surrounding the magnetite-hematite-maghemite grains and replacing hematite. Barite and uraninite occur in fractures of iron oxides, Mn-Ba-Pb oxi-hydroxides and scheelite occur interstitial to iron oxides.

Introduction

Iran is the 10 th largest iron ore producer and the 9 th iron ore exporter in the world. According to forecasts, iron production will triple until 2025 up to 66.2 Mt of iron ore [START_REF] Hastorun | Recent trends in the nonfuel minerals industry of Iran[END_REF]. Most of the iron deposits are located in the northeast and central part of Iran, while iron ore deposits of the northwest part are still little studied. We carried out a preliminary mineralogical study on different iron ore samples from the Takab region in north western Iran (Fig. 1) in order to define iron oxide minerals in their structural and fluid context.

Geological setting

The Takab iron ore deposits belong structurally to the Sanandaj-Sirjan Zone (Fig. 1). The study area comprises Precambrian-Early Paleozoic metamorphic rocks, Miocene sedimentary rocks and Quaternary deposits. The Precambrian rocks have been thrust over the Cenozoic sedimentary units. The iron ore deposits of this area are classified being of volcano sedimentary origin. Their formation might be related to the back arc extension do to the Prototethys subduction beneath the Cadomian magmatic arc [START_REF] Ghorbani | The Economic Geology of Iran: Mineral Deposits and natural Ressources[END_REF][START_REF] Nabatian | Iron and Fe-Mn mineralisation in Iran: implications for Tethyan metallogeny[END_REF]. The most important iron deposits, known at present, occur in schists and marbles. The mineralization forms thin layers and lenses interlayered with quartzite [START_REF] Nabatian | Iron and Fe-Mn mineralisation in Iran: implications for Tethyan metallogeny[END_REF]; Fig. 2A).

The iron oxide mineralogy is essentially described as hematite, goethite, with barite as accessory mineral [START_REF] Ghorbani | The Economic Geology of Iran: Mineral Deposits and natural Ressources[END_REF]. This is the first detailed petrological mineralogical study on banded and nodular iron ores from the Takab region reflecting its complex geological history.

Analytical Methods

Optical microscopy on polished thin sections was performed in transmitted and reflected light mode, SEM and electron microprobe analyses at the Centre Camparis, Sorbonne Université, UPMC Univ Paris 06, Paris, France. XRD and Curie Balance analyses were carried out at the GEOPS, Université Paris Saclay, Université Paris Sud, Orsay, France.

Results

The studied samples are banded and nodular iron ores.

Both massive banded and nodular ore occur as decimetric layers within deformed (folded) mica schist (Fig. 2A). banded iron ore shows coarse and discontinuous banding with interstitial quartz (Fig. 3A). The nodular ore is composed of mm-sized iron oxide agglomerates, partly elongated and disrupted (Fig. 3B).

The banded iron ore is characterized by micro quartz showing 120° triple junctions or slight grain boundary migration of lobate crystal surfaces. Some quartz grains show also 120° triple junctions with iron oxide grains. Quartz grains show undulose extinctions. In iron oxide rich bands, poikilotopic textures of iron oxides or quartz grain are characteristic. Iron oxide veins formed perpendicular to the iron oxide bands (Fig. 4A). Pyrite relicts are observed in goethite. Barite, interstitial to quartz, is rarely present.

The nodular iron ore shows bimodal iron oxide and quartz grain sizes: 1) some tenths of micron and 2) some hundreds of microns. Both quartz types show undulose extinction (Fig. 4B). The large quartz grains occur mainly in the vicinity of iron oxide agglomerates of a few hundreds of microns. They are composed of individual small oriented quartz crystals. Small iron oxide grains show sometimes 120° triple junctions with quartz grains.

Grain boundary migration also led to grain isolation (Fig. 4B). X-ray diffraction on 3 banded and 1 nodular ore samples indicate that all samples contain similar major mineralogy quartz, magnetite, hematite, maghemite, goethite in different proportions. Small peaks are attributed to witherite (BaCO3).

The iron oxi-hydroxide mineralogy is complex. Euhedral magnetite is transformed into maghemite and hematite (Fig. 4C, D, 5A, B). Trellis structures, which are characteristic for martite, are also observed (Fig. 5E). Chemical zoned goethite occurs interstitial to the magnetite-hematite-maghemite grains (Fig. 5C). Goethite may also host neoformed hematite inclusions (Fig. 5B). Euhedral, zoned and corroded zircon, as well as euhedral monazite are observed as inclusions in magnetite (Fig. 5D,E). In the matrix, a few hundreds of micron large oxi-hydroxides composed of Mn-Ba-Pb occur (Fig. 5F). Scheelite (CaWO4) was also observed. Secondary barite and uraninite occur in fractures of the iron oxides together with goethite.

Curie Balance analyses show that the material is highly magnetic and magnetite is the dominate phase. In two banded iron ore samples maghemite and /or goethite are indicated.

Electron microprobe analyses were performed on iron oxihydroxides from banded and nodular ores. A special program was designed to analyse oxygen in order to trace the iron oxide and hydroxide variabilities. The magnetite and hematite contain traces of Al (< 0.5 wt. %). Goethite hosts variable amounts of Al (up to 0.7 wt. %), Si (up to 1.7 wt. %) and Ca (up to 0.2 wt. %). Figure 6 shows the iron versus the oxygen contents of representative iron oxhydroxides and the iron contents versus the total (wt. %).

The conversion of Fe 2+ into Fe 3+ is accompanied by a decrease of the total iron content and a slight increase of oxygen. Goethite has lowest totals.

Discussion and Conclusions

This study highlights the complex history of this iron ore and gives useful hints for future research, which will contribute to elucidate the origin of the ore forming elements iron and silica, but also the accessory elements (Mn, Ba, Pb, U, W…). Based on the present results, it can be suggested that: (1) The primary iron ore is banded and crystallized statically from an iron and silica rich solution.

Zircon and monazite were included during this step in the initial iron oxide, magnetite; (2) Regional deformation led to recrystallization of the iron oxides through grain boundary migration and quartz grain isolation. The micro quartz grains form neoblasts and annealed to the larger quartz grains. The nodular ore represents most likely a highly deformed banded iron ore, where iron oxide bands were disrupted and recrystallized. Sulfides were present at that state. Oxidizing fluids accompanied this step and transformed the magnetite into maghemite and hematite.

(3) Aqueous solutions carrying Mn, Ba, Pb, Ca, W, U, S, CO2 led to the formation of goethite and Mn-Ba-Pb oxihydroxides and associated scheelite, barite and uraninite. This episode is most likely contemporaneous to the W mineralisation in this region [START_REF] Ghorbani | The Economic Geology of Iran: Mineral Deposits and natural Ressources[END_REF].

Figure 1 .

 1 Figure 1. Tectonic and structural map of Iran (modified after Alavi 1991; Aghanabati 1998, 2005). The location of the study area is shown by the black star. Abbreviations are as following: Zagros, Zagros ranges: KRSZ, Kermanshah Radiolarites subzone; SSZ, Sanandaj-Sirjan magmatic-metamorphic zone; UD: Urumieh-Dokhtar magmatic arc; CIM, Central Iranian microcontinent (includes the Yazd, Posht-e-Badam (PB), Tabas, and Lut blocks); Alborz, Alborz ranges, western Alborz-Azarbayjan, KTZ Khazar-Talesh-Ziveh structural zone; CIZ, Central Iranian zone; Sistan, East Iran ranges; Makran, Makran zone; KD, Kopeh-Dagh ranges; Zabol, Zabol area; and TMR, Cenozoic magmatic rocks.

Figure 2 .

 2 Figure 2. A. Iron ore layer (mainly magnetite) in mica schist. B. fine banded iron ore; C. nodular iron ore (Takab region, NW of Zanjan).

Figure 3 .

 3 Figure 3. Polished thin section scans. A. Banded iron ore composed of alternating quartz and iron oxide bands; B. Nodular iron ore in quartz matrix.

Figure 4 .

 4 Figure 4. A+B: Optical Microscopy (transmitted light). A. Micro quartz hosting iron oxide bands and an iron oxide vein perpendicular to banding. Quartz grains show grain boundary migrations and rarely 120° triple junctions; B. Iron oxide nodules and large quartz grains showing undulose extinction at large grain scale. The large quartz grains are composed of individual small elongated quartz grains. C+D: Reflected light microscopy. C. Iron oxide grain in quartz matrix showing hematite (he, bright) in magnetite (mag). Goethite occurs between the iron oxide grain and the quartz grains; D. Close up of an iron oxide grain showing magnetite (Mag, grey), hematite (He, bright) and probably maghemite (Magh, light grey).

Fig. 5 .

 5 Fig. 5. Scanning electron microprobe images in back scattered mode. A. Euhedral magnetite partly transformed into hematite and surrounded by goethite; B. Goethite interstitial to magnetite with euhedral hematite inclusions. C. Goethite showing oscillatory zoning; D. Euhedral magnetite with corroded zircon (zr) inclusion; E. Magnetite showing trellis hosting euhedral monazite; F. Mn-Ba-Pb oxi-hydroxide interstitial to magnetite. Magnetite hosts uraninite inclusion.

Figure 6 .

 6 Figure 6. Electron microprobe analyses of iron oxides and iron hydroxides of some iron ores from the Takab region. Upper diagram: Fe versus O (wt. %), lower diagram: Fe versus Total (wt. %). Centre Camparis, Sorbonne Université, UPMC, Univ Paris 06, France).

Acknowledgements

We thank the Institute of Advanced Sciences of Basic Science and University of Zanjan, Zanjan for logistic support during fieldwork and the stay of BO in Zanjan. Campus France and the French embassy are thanked for financing the stay of BO at IASBS, Zanjan. The UMR7193 ISTeP, UPMC is thanked for financial support of the analyses.


Sébastien Kolb
email: sebastien.kolb@ecole-air.fr

Olivier Montagnier
email: olivier.montagnier@ecole-air.fr

Laurent Hetru
email: laurent.het@free.fr

Thierry Faure
email: thierry.faure@ecole-air.fr

Real-Time Detection of an Aircraft Deep Stall and Recovery Procedure

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Real-time detection of an aircraft deep stall and recovery procedure

Sébastien Kolb * , Olivier Montagnier † , Laurent Hétru ‡ and Thierry M.

II. Introduction

D stall is a kind of stall affecting aircraft longitudinal dynamics. It corresponds to a stable equilibrium at high angle-of-attack (AoA) for which the pitch control surface is almost or totally ineffective [START_REF] Torenbeek | Flight Physics: Essentials of Aeronautical Disciplines and Technology, with Historical Notes[END_REF][START_REF] Abzug | Airplane Stability and Control : A History of the Technologies that Made Aviation Possible[END_REF]. Deep stall is a hazardous aircraft loss-of-control which implies an established stall flight with high descent velocities and no easy recovery procedure. In the history of aeronautics, several aircraft have encountered this type of stall leading to crashes, more especially during flight tests (e.g. Javelin Gloster, BAC 1-11, Canadair Regional Jet 600 [START_REF] Roberts | Safety recommandation A-94-101[END_REF], etc.).

The existence of this equilibrium depends on aircraft design factors and on some peculiar flight conditions. Taylor and Ray [START_REF] Taylor | A systematic study of the factor contributing to post-stall longitudinal stability of T-Tail Transport configurations[END_REF][START_REF] Taylor | Deep Stall aerodynamic characteristics of T-Tail aircraft[END_REF] have studied several aircraft designs in a wind tunnel to identify geometrical factors promoting deep stall risk. They noticed that T-tail or aft mounted engines are critical factors. The flight of the T-tail inside the separated wake of the stalled main wing is generally considered as an explanation of the phenomenon [START_REF] Taylor | A systematic study of the factor contributing to post-stall longitudinal stability of T-Tail Transport configurations[END_REF][START_REF] Taylor | Deep Stall aerodynamic characteristics of T-Tail aircraft[END_REF][START_REF] White | Simulator studies of the deep-stall[END_REF][START_REF] Faure | Aerodynamic features of a two-airfoil arrangement[END_REF]. Byrnes et al. [START_REF] Byrnes | Effect of horizontal stabilizer vertical location on the design of large transport Aircraft[END_REF] have studied the impact of T-tail position and surface, and swept wing angle on the aircraft aerodynamics and flight dynamics. In these experiments, deep stall equilibrium occurs for T-tail and classical small tail configurations but only in the case of swept main wings. Swept wings induce pitch-up effect after stall leading to deep stall equilibrium [START_REF] Kenneth | Design Guide for Pitch-Up Evaluation and Investigation at High Subsonic Speeds of Possible Limitations Due to Wing-Aspect-Ratio Variations[END_REF][START_REF] Harper | A Review of the Stall Characteristics of Swept Wings[END_REF].

Nevertheless, in the last decade, several aircraft (e.g. Airbus A400M, Comac ARJ21, Beechcraft 1900, Bombardier CRJ1000, etc.) and also small unmanned aerial vehicle (UAV) have been designed with a T-tail and sometimes with aft mounted engines. In general, the normal laws of aircraft are protected from high AoA (AoA protection) as in the case of the C-17A [START_REF] Iloputaife | Design of a Deep Stall Protection for the C-17A[END_REF]. Other airplanes have a stick pusher or a stick shaker to prevent such dangerous situations. In spite of everything, some flights involve high AoA like flight tests for a new aircraft which necessarily reach the stall or post stall regions. For an aircraft potentially presenting a deep stall risk, a "deep stall recovery system" (rocket, parachute, etc.) is generally implemented to get out from deep stall equilibrium if necessary [START_REF] Wall | A400M To Shed Strain Gauges That Slow Flight Testing[END_REF]. But this kind of device can be fatal if it is activated at the wrong time instant. Another case concerns peculiar system failures involving alternate law (or secondary law) where AoA protection is switched off as in the flight crash AF447 Rio-Paris. The final report on the accident [START_REF]Final report on the accident on 1st June 2009 to the Airbus A330-203 registered F-GZCP operated by Air France flight AF 447 Rio de Janeiro[END_REF] states that the AoA had reached 40°. Lastly, in normal flight conditions, some inertial and aerodynamic aircraft configurations, possibly combined, can also lead to deep stall entry such as bad aft center balance, overload, unexpected high-lift devices use.

Therefore, in accordance with Goman [START_REF] Goman | Application of bifurcation methods to nonlinear flight dynamics problems[END_REF], it is interesting to know as soon as possible the precursory signs of deep stall entry in order to implement a warning for this kind of untypical situation follow-up and an automatic or manual procedure to escape from this hazardous phenomenon. Recently, several works have concerned the analysis of deep stall dynamics applying bifurcation theory [START_REF] Wang | Study of deep-stall characteristics and longitudinal special phenomena of T-tail aircraft[END_REF][START_REF] Shi | Bifurcation analysis of polynomial models for longitudinal motion at high angle of attack[END_REF][START_REF] Kolb | Nonlinear analysis of the longitudinal flight dynamics of a transport airplane[END_REF][START_REF] Hétru | Tandem airfoil vortex interaction in deep-stall[END_REF][START_REF] Xin | Bifurcation analysis and stability design for aircraft longitudinal motion with high angle of attack[END_REF][START_REF] Kolb | Nonlinear analysis and control of an aircraft in the neighbourhood of deep stall[END_REF], others have concerned the analysis of deep stall landing for UAV [START_REF] Taniguchi | Analysis of Deepstall Landing for UAV[END_REF][START_REF] Kubo | High angle of attack flight control of a tail-sitter unmanned aircraft[END_REF][START_REF] Mathisen | Non-linear Model Predictive Control for Guidance of a Fixed-Wing UAV in Precision Deep Stall Landing[END_REF] but none concern the real-time detection of deep stall and recovery procedure. On the other hand, due to the development of more and more autonomous aircraft, many recent works concern the loss-of-control prevention, detection and recovery [START_REF] Richards | Vehicle Upset Detection and Recovery for Onboard Guidance and Control[END_REF][START_REF] Lombaerts | Piloted Simulator Evaluation of Safe Flight Envelope Display Indicators for Loss of Control Avoidance[END_REF][START_REF] Lombaerts | Design and piloted simulator evaluation of adaptive safe flight envelope protection algorithm[END_REF]. Deep stall studies are generally carried out with the assumption of longitudinal flight and such hypothesis is adopted hereafter. However, an asymmetric phenomenon like the roll off instability could exist. In fact, pilots try first to find a solution keeping the longitudinal flight and hesitate to escape from a crash thanks to a lateral procedure with the risk of spin.

The aim of the present study is to propose a real-time identification of deep stall entry as soon as possible and find the ad-hoc way out without knowing the aircraft aerodynamic above the stall angle. The proposed algorithm must be simple and robust because it cannot obviously be tested on a real airplane (unless on a small UAV or an aircraft designed for that purpose [START_REF] Sim | Flight characteristics of a manned, low-speed, controlled deep stall vehicle[END_REF]). To avoid a false detection in normal flight condition, it would be only triggered when the aircraft would be in alternate law or in flight test phase below the stall angle. Therefore, this study concerns an actual existing aircraft but not the conception parameters involved during the design phase. In order to reach this goal, the first step is to analyze the path towards the deep stall equilibrium and the second is to identify the dynamical features of the phenomenon. The first part of the paper is about the implementation of a realistic model of a deep stall prone aircraft. The second part concerns the analysis of the aircraft model behavior applying dynamical system theory in order to have time evolution of the physical variables and modes properties near the equilibrium points. The following part proposes the real-time deep stall detection algorithm. Finally, the real-time detection is successfully applied and some recovery procedures are evaluated.

III. Aircraft model

A. Governing equations

As it is considered classically, it is assumed here that aircraft motion is longitudinal during deep stall phenomenon. In the longitudinal plane, the earth, aerodynamic and aircraft frames of reference are denoted (x 0 , z 0 ), (x a , z a ) and (x b , z b ), respectively. In consequence, the nonlinear governing equations of the aircraft are in the aerodynamic frame for the forces and the aircraft frame for the moment (see the nomenclature for the list of the symbol definitions):

V = - ρV 2 S 2m C D (α, q, δ e , δ f ) + T(ρ, δ T ) m cos α -g sin γ (1) γ = ρV S 2m C L (α, q, δ e , δ f ) + T(ρ, δ T ) mV sin α - g V cos γ (2) α = q - ρV S 2m C L (α, q, δ e , δ f ) - T(ρ, δ T ) mV sin α + g V cos γ (3)
q = ρV 2 Sc 2I y C m (α, q, δ e , δ f ) (4)
h = V sin γ ( 5
)
where it is assumed that thrust does not create any moment and is aligned with the x b direction. This mathematical model of longitudinal flight dynamics can be rewritten in the state space form:

X = f (X, U) (6)
where the state vector X and the control vector U are:

X = (V, γ, α, q, h) ; U = δ e , δ T , δ f (7)
Hereafter, computation of Eqs. 1-5 is called full model and computation of only Eqs. 3-4 is called reduced model.

B. Aerodynamic coefficients

For the deep stall study, it is necessary to find realistic aerodynamic coefficients at very high AoA of a deep stall prone aircraft. These data are easier to find for a fighter aircraft in literature [START_REF] Goman | Application of bifurcation methods to nonlinear flight dynamics problems[END_REF][START_REF] Nguyen | Simulator Study of Stall/Post-Stall Characteristics of a Fighter Airplane With Relaxed Longitudinal Static Stability[END_REF] than for a transport aircraft. Nevertheless, the Learjet 23 D is a good candidate, because it was a well known deep stall prone aircraft and because Soderman and Aiken have obtained the static aerodynamic coefficients during full scale wind-tunnel tests [START_REF] Soderman | Full scale wind tunnel tests of a small unpowered jet aircraft with a t-tail[END_REF]. In consequence, all static coefficients C L , C D and C m s t at i c (for q = 0) were digitized hereafter. These experimental measurements have been interpolated with spline functions in the experimental data range and extrapolated linearly outside the bounds. Figures 1a and1b give the lift C L and drag C D coefficients for a large range of AoA, respectively without flaps (δ f = 0°) and with full flaps (δ f = 40°) configurations and extreme elevator deflection angle values. Figure 1c gives C m s t at i c for the same configurations at wing quarter chord.

The static part of the pitching moment allows to determine the AoA values at equilibrium (Fig. 1c). From the flight dynamics point of view, deep stall occurs when there exists a stable equilibrium at high AoA for a given elevator deflection angle, center of gravity and configuration (flaps up/down, etc). A longitudinal equilibrium requires, among others, a zero pitching moment and is statically stable when the α-slope of the pitching moment coefficient is negative (i.e. ∂C m /∂α < 0) [START_REF] Etkin | Dynamics of Flight : Stability and Control[END_REF].

Moreover, the loss of tail efficiency in deep stall also affects the pitch rate derivative of the pitching moment coefficient (∂C m /∂q ). Indeed the lift of the horizontal tail creates normally an opposite moment to the aircraft pitch motion, followed by a modification of the local AoA of the tail versus that of the wing. But this effect is reduced with the degraded aerodynamics of the horizontal tail, when it is inside the separated wake of the stalled main wing. Here, a simple continuous piecewise linear function is proposed for the pitch rate derivative inspired from the model of Montgomery and Moul [START_REF] Montgomery | Analysis of deep-stall characteristics of T-Tailed Aircraft Configurations and Some Recovery Procedures[END_REF] and whose absolute value is reduced in the deep stall conditions:

∂C m ∂q = C m q c V where C m q =          -9.4 α [20 • , 50 • ] 0.61 × (α -20) -9.4 α ∈ [20 • , 35 • ] -0.61 × (α -35) -0.25 α ∈ [35 • , 50 • ] (8)
It should be noticed that recent experimental studies on SUPRA aircraft project showed that the C m q increase could be higher than the one proposed in Montgomery and Moul [START_REF] Abramov | Aerodynamic model of transport airplane in extended envelope for simulation of upset recovery[END_REF][START_REF] Khrabrov | Aerodynamic model development and simulation of airliner spin for upset recovery[END_REF]. The overall aircraft pitching moment takes the following form with its diverse contributions, that is to say the quarter chord pitching moment C m static , the static margin δ sm and the pitch damping C m q : with

C m α, q, δ e , δ f = C m static α, δ e , δ f + C L (α, δ e , δ f )δ sm + c V C m q (α)q (9)
δ sm = x cg -x ac c (10)
where x cg and x ac are the algebraic position of the center of gravity and the wing aerodynamic center as in [START_REF] Soderman | Full scale wind tunnel tests of a small unpowered jet aircraft with a t-tail[END_REF]. With this definition, the static margin is negative for an aft center of gravity. The aircraft pitching moment is given for q = 0 and δ sm = -10% in Fig. 1d. Thrust in Newtons is defined using the following expression: [START_REF] Soderman | Full scale wind tunnel tests of a small unpowered jet aircraft with a t-tail[END_REF][START_REF] Stengel | Flight Dynamics[END_REF]. The thrust contribution on the aircraft pitching moment is neglected hereafter. It can be noted that for a T-tail aircraft with engines on the rear fuselage, the moment is significantly smaller than in a conventional setup where the engines are placed under the wings. Finally, it is assumed as usual that ∂C L /∂q = ∂C D /∂q = 0.

T = 26243.2( ρ 1.225 ) 0.

IV. Analysis of aircraft model behavior

In this part, characteristics properties of the aircraft dynamics near deep stall are identified in order to predict the convergence of its trajectory in a phase portrait toward a low or a high AoA equilibrium.

A. General aircraft dynamics

An aircraft may meet deep stall for several ranges of elevator deflection angles in a given configuration (static margin, flaps, etc.). For such critical cases, there are three equilibria for one fixed elevator deflection angle δ e . They correspond to the AoA for which the overall pitching moment is equal to zero (C m (α, δ e ) = 0 on Fig. 1c for δ sm = 0% and on Fig. 1d for δ sm = -10%). These equilibria are: the classical low AoA equilibrium, an unstable equilibrium at medium AoA (called saddle point) and the stable equilibrium at high AoA associated to deep stall. These points are marked by color dots for δ f = 40°and δ e = 0°in Fig. 1d. Figure 2a and 2b show the phase portrait of the Learjet 23D, typical of a deep stall prone aircraft, computed with the reduced model and the full model, respectively. A stable manifold can be computed with the reduced model (thick black curve in Fig. 2a), it starts from the saddle point and delimits the basin of attraction of the equilibrium at high AoA. Note that the stable manifold is computed using a backward time integration starting from a point in the neighborhood of the saddle point. Two types of behaviors can be distinguished in Fig 2a . First, when the airplane converges towards a stable equilibrium point, the distance between two angles-of-attack for which α = 0 becomes more and more reduced. Secondly, in the neighborhood of the saddle point, all the trajectories are attracted by this equilibrium and then repelled from it. The trajectories remaining at the right of the saddle point equilibrium stay in the basin of attraction as shown by several authors [START_REF] Goman | Application of bifurcation methods to nonlinear flight dynamics problems[END_REF][START_REF] Montgomery | Analysis of deep-stall characteristics of T-Tailed Aircraft Configurations and Some Recovery Procedures[END_REF][START_REF] Thomas | Longitudinal motions of aircraft involving high angles of attack[END_REF]. Other trajectories (outside the basin of attraction) converge towards the low AoA equilibrium. Another simple remark never pointed out in previous studies, to the best of our knowledge, is the way of convergence toward a stable equilibrium. Figure 2a shows a fast convergence of the short-period mode to the low AoA in only one cycle while the convergence to the high AoA is slow and takes many cycles. Slow convergence can be considered like an opportunity allowing sufficient time to initiate the recovery procedure. This difference is explained by the large difference in damping for the two cases that will be explained and quantified in part IV.C.

Another point is that the reduced model computations are carried out with the hypotheses of constant slope and velocity. For the low AoA equilibrium, the slope is necessarily null, corresponding to the light green basin of attraction. On the contrary, for the high AoA equilibrium, the aircraft gets a very negative slope. For the studied configuration, the airplane is on a -33°slope and its velocity is decreased to 63 m/s. These conditions give a new stable line (dotted black in Fig. 2a) and a new basin of attraction. It can be noticed that both basins of attraction are very close. Finally the slope has little influence on the airplane behavior. Hereafter, only the basin of attraction for a slope equal to zero is considered since this is the one indicating whether or not the airplane is entering in deep stall. A second point is that the basin of attraction can be plotted only for the reduced model, because only the planar case insures that the trajectories are not crossing each others [START_REF] Montgomery | Analysis of deep-stall characteristics of T-Tailed Aircraft Configurations and Some Recovery Procedures[END_REF][START_REF] Guckenheimer | Nonlinear Oscillations Dynamical Systems and Bifurcations of Vector Fields[END_REF]. This is what Fig. 2b shows, where the full model is considered. Then, it can be noticed that the basins of attraction previously computed provide a blurred criterion for the convergence or divergence toward the high AoA equilibrium.

Fig. 3 presents an airplane at low AoA equilibrium whose pilot applies a pitch-up command during 1s followed by a return to δ e = 0. In this realistic airplane flight with a change on pilot command, the basin of attraction computed with the reduced model provides a good estimate of the final equilibrium for both the reduced and full calculations. In addition, the basin of attraction for a full pitch-down command is plotted in the figures (in light red), showing that there is no possible recovery of the airplane after one cycle. Figure 4 presents the bifurcation diagram corresponding to the AoA at equilibrium as a function of the elevator position δ e . It shows two limit points (LP), also called saddle-node bifurcations or fold bifurcation, at δ e = -10°and δ e = 5°delimiting three zones [START_REF] Guckenheimer | Nonlinear Oscillations Dynamical Systems and Bifurcations of Vector Fields[END_REF]. For an elevator deflection angle below -10°and above 5°, there is only one stable equilibrium. For an elevator deflection angle between -10°and 5°, there are three equilibria, two stable ones and an unstable one. Near an elevator deflection angle of δ e ¬ ¬ -10°in a stabilized situation at low AoA, a pitch-up command can produce a jump to the stable branch at high AoA which will not be so easy to recover from since a reverse pitch-down command may change the equilibrium state value, but the aircraft may remain at high AoA. There is a phenomenon which looks like an hysteresis. It is a hazardous situation and as long as the aircraft is on the stable branch at low AoA, the pilot may expect and foresee the happening of change of AoA leading to deep stall.

B. Effect of altitude, velocity, flaps and pitch damper

The basins of attraction for diverse configurations of the Learjet 23 D are now compared.

Figure 5a shows the basin of attraction with flaps up for various altitudes between 3000 m and 10000 m. The deep stall basin of attraction is rather limited. The altitude has little influence on its extension around the high AoA equilibrium point. However, the branch at the left top corner is influenced by the altitude. It means that the AoA for which the aircraft can enter into the basin of attraction (Fig. 3) is lower at high altitude (h = 10000 m). Figure 5b shows the basin of attraction with flaps down for the same altitudes but for a slightly lower speed corresponding to the maximum speed the flaps down configuration can bear. Larger basins of attraction are observed, showing that this configuration is more deep stall prone. As a consequence, it is advisable to retract the flaps when a situation prone to deep stall is detected for this type of T-tail aircraft.

Figure 5c shows the basin of attraction with flaps up for various velocities between 80 m/s and 160 m/s at h = 6000 m. The velocity has little influence on the basin of attraction except for 80 m/s which corresponds to the region of the reversed command of the power curve.

Effects of pitch damping on the aircraft behavior have been widely studied in the past by Thomas and Collingbourne [START_REF] Thomas | Longitudinal motions of aircraft involving high angles of attack[END_REF]. Figure 5d shows the basin of attraction with pitch damper for various altitudes between 3000 m and 10000 m and a velocity of 80 m/s. The basin of attraction is larger when the pitch damper is engaged (comparison with Fig. 5b). Thus, it can be advised to switch off the pitch damper when a situation prone to deep stall is detected. This effect can be explained by the design of a pitch damper: a component of elevator deflection angle δ e which is added to the pilot action and proportional to the pitch rate q. It produces more pitch damping which is translated mathematically into a more negative derivative of the pitching moment coefficient with respect to the dimensionless pitch rate C m q . Moreover, since the oscillations are less damped out and the divergence of the vector field is lower when the pitch damper is switched off, the aircraft is less likely to converge towards the deep stall equilibrium and since the damping of the movement is lower, the pilot has a little more time available to react.

C. Short-period mode properties at low and high AoA

The relatively fast motion observed in Fig. 2a corresponds to short-period modes. These modes are classically assumed decoupled with phugoid modes at low AoA [START_REF] Etkin | Dynamics of Flight : Stability and Control[END_REF][START_REF] Stevens | Aircraft Control and Simulation[END_REF] but also here at high AoA [START_REF] Hetru | Étude expérimentale et numérique de l'interaction aérodynamique entre deux profils: application au risque aéronautique du décrochage profond[END_REF]. Note that the simulations of the reduced model seem quite close to the ones of the full model (Fig. 3).

By linearizing the dynamics of the lift and pitch equations near a low AoA equilibrium for a fixed airspeed V e , altitude and angle-of-attack α e , analytical expressions of the pulsation ω spm and of the damping ζ spm can be obtained classically for the short-period mode [START_REF] Etkin | Dynamics of Flight : Stability and Control[END_REF][START_REF] Stevens | Aircraft Control and Simulation[END_REF]:

ω 2 spm = - V 2 e ρcS 2I y ρSV e 2m C L α + T mV e c V e C m q + C m α ( 11
)
2ζ spm ω spm = ρSV e 2m C L α - c 2 ρSV e 2I y C m q + T mV e (12)
where C L α , C m α and C m q = (V e /c) × ∂C m /∂q are computed for the values: α = α e , q = 0, δ e and δ f . It can be shown that a quasi-similar expression can be obtained in the case of high AoA [START_REF] Hetru | Étude expérimentale et numérique de l'interaction aérodynamique entre deux profils: application au risque aéronautique du décrochage profond[END_REF]:

ω 2 spm = - V 2 e ρcS 2I y ρSV e 2m C L α + T mV e cos α e c V e C m q + C m α ( 13
)
2ζ spm ω spm = ρSV e 2m C L α - c 2 ρSV e 2I y C m q + T mV e cos α e (14)
but the values of derivatives C L α , C m q are very different because they are measured near α e which is very large in this case. Table 1 shows an example of short-period mode computation for the Learjet 23 D with the Eqs. 11-14. This numerical application shows clearly that the damping of the short-period mode is far smaller at high AoA than at low AoA (approximately 10 times smaller). This phenomena can be explained by the fact that in stall, the lift curve slope C L α is lower (see for example the slope of C L between 15°and 40°vs. the slope at origin on Fig. 1a which is largely reduced or even negative). Moreover, the absolute value of the pitching moment coefficient derivative ( C m q ) is also clearly reduced (see C m q definition in the part III.B). On the other hand, the pulsation is quite the same at low or high AoA and cannot be used as a discriminant condition to find the type of equilibrium. Note also that the period of the short-period mode at high AoA is approximately equal to 4 s allowing to confirm that the convergence to the deep stall equilibrium is relatively slow and the order of magnitude is around one minute (see Fig. 2a).

V. Real-time deep stall detection

Previous main remarks are now used to propose a fast detection algorithm of deep stall without knowing the aircraft aerodynamics at medium and high AoA. The fast detection permits, among others, to have a sufficient level of energy to go out from deep stall and to react as early as possible, as far as possible from the deep stall equilibrium. The main idea of this algorithm is to predict the evolution of the aircraft trajectory in the (α, α) phase portrait using the relative positions of the extrema. Some assumptions should be made before describing the method. Firstly, it is assumed that α and α (or q) are measured in real-time with sufficient confidence up to very high AoA. Even if the data are noisy, it is also assumed that a standard filter allows to obtain sufficiently smooth estimation of these data. Secondly, the aircraft is assumed free and in a quasi longitudinal flight (quasi-symmetric aircraft motions in the vertical plane). Free means that the pilot makes no supplementary action to the pitch command during detection.

To avoid a false detection in normal flight condition, the algorithm is started only when α exceeds α st all but a different choice can be made. Starting from the stall angle, Fig. 6 shows two typical trajectories. The first one converges to a low AoA and the second one to a high AoA. As it was said before, in the two cases, the aircraft is repelled from the saddle point equilibrium. On the figures, α is incrementally denoted α i (i ∈ N * ) when α = 0 and α d j ( j ∈ N * ) when d α/dα = 0. From a mathematical point of view, the main difference between the portraits of Figs. 6a and6b is the number of α d j between two successive α i (α 1 and α 2 in the figures). Consequently, the case of the convergence to a low AoA (Fig 6a) can be easily eliminated with the condition max( j) = 1. In the second case (Fig 6b) and if α d1 ∈]α 2 , α 1 [, aircraft begins to describe a spiral trajectory typical of deep stall. Another step is then necessary to confirm this spiral trajectory and the convergence into a high AoA. The analysis of what happens between α 2 and α 3 can confirm this fact. If there is again only one α d j (α d2 in Fig. 6b) i.e. max( j) = 2 and α d2 ∈]α 2 , α 3 [ and also α 3 ∈]α 2 , α 1 [, the detection of a high AoA equilibrium is more and more probable. Then, an approximation of the AoA at equilibrium and the damping value of the short-period mode can be computed with the mean and the half-cycle turning points approximation formulae:

α e = α d1 + α d2 2 (15)
ζ spm = 1 2π ln α 1 -α 2 α 3 -α 2 (16)
Finally, the detection is considered true if α e α st all and ζ spm ζ L Ao A spm . The identification algorithm is summarized in Fig. 7. Note that the algorithm distinguishes different extrema order if α 1 > α 2 or α 1 < α 2 . Another remark is that the algorithm is written with the (α, α) state values but it could also have been done with the (α, q) system. The last comment is that other kinds of trajectories in the (α, α) portrait can be found in the literature [START_REF] Thomas | Longitudinal motions of aircraft involving high angles of attack[END_REF][START_REF] Wang | Algebraic analysis of stability and bifurcation for nonlinear flight dynamics[END_REF]. However, it has been checked that any other trajectory than spiral is rejected by the present simple algorithm. (reduced model, flaps down δ f = 40°, δ e = 0°, δ T = 1 and δ sm = -10% with initial conditions V 0 = 72 m/s, γ 0 = 0°, α 0 = 44°, 46°, q 0 = 0 rad/s and h 0 = 10 000 m).

VI. Application of the real-time detection algorithm and recovery procedure

In this section, the prediction algorithm is applied on the aircraft model and recovery is performed with a full pitch-down command applied manually or automatically.

In the past, Thomas and Collingbourne [START_REF] Thomas | Longitudinal motions of aircraft involving high angles of attack[END_REF] carried out a large numerical study to identify the entry and exit maneuvers for a deep stall prone aircraft and noticed a strong influence of initial conditions. Figure 8 shows two kinds of pilot's reaction for an airplane in the deep stall basin of attraction. In Fig. 8a, the pilot is supposed to understand the airplane deep stall after 40 s and applies a full pitch down command of δ e = 15°. This is corresponding to an estimated reaction time necessary to understand this unusual flight situation. He does not succeed in going out from the deep stall equilibrium because, after 40s, the position in the phase portrait is completely inside the rose basin of attraction corresponding to δ e = 15°. The aircraft reaches a new deep stall equilibrium with a little lower AoA. In Fig. 8b, the pilot is warned of a fore-coming deep stall by the algorithm (starting from the initial point of the simulation) after approximately 5 s and reacts 5 s later by applying a large pitch-down command of δ e = 15°during 5 more seconds and then release the command. This pilot succeeds in going out from deep stall and reaches a low AoA. This example shows the usefulness of being aware of the dangerous situation as early as possible but it is not sufficient. The next two examples in Fig. 9 show that the time instant when the pitch-down command is applied is as well as important. In the first case (Fig. 9a), the pilot applies the command just after that the AoA reaches the first turning point α 4 , that is to say when both α and α decrease. In this case, the pilot succeeds in going out from deep stall. In the second case (Fig. 9b), the pilot applies the command just after that the AoA reaches the second turning point α 5 , that is when both α and α increase. The pilot does not succeed in coming back to a low AoA even if the command is maintained during 10 s.

The last example shows that, in order to avoid the basin of attraction (δ e = 15°) where deep stall recovery is impossible, it is better to apply a pitch-down command when both α and α decrease. In order to justify this statement, a simplified model for the short-period mode is written in the following form (ignoring the term with C L α ):

α = q (17) q = ρV 2 Sc 2I y C m s t at i c (α, δ e ) + c 2V C m q (α) q (18)
with C m q (α) < 0 and dC m st at i c /dδ e < 0. When applying a pitch-down command (δ e > 0), the variation of the pitching moment is negative ∆C m st at i c (α, δ e ) < 0. Thus the pitch rate derivative q is more negative and the pitch rate decreases even more. In order to be in a situation where the AoA reduces much more and such that it flies easier towards lower AoA, it is necessary that both derivatives of AoA and of pitch rate are negative, which means that both variables α and q decrease. This case is shown in Fig. 9a. This part demonstrates the interest of the detection algorithm but also that the choice of the time instant when the pitch-down command is applied is fundamental to recover the aircraft. Because it is not possible to expect the pilot to apply the command at the right time instant, it can be concluded that the success of the recovery necessitates to automate this last step. It consists simply in applying the larger pitch-down command when α 3 or α 4 is reached, depending on the order of α i and until α st all threshold is crossed. The basic recovery algorithm is summarized in Fig. 10. In order to validate the detection and recovery algorithm for different initial conditions, simulations have been done for three aircraft configurations (aft CG and flap up, neutral CG and flap down, aft CG and flap down) with the reduced model or the full model in Fig. 11. The procedure is the same as in Fig. 3 but the detection and recovery algorithms are now applied. Both models (reduced and full) provide similar behaviors. Simulations show that there exist three kind of situation for the detection algorithms:

• Figs. 11a-11b. With flap up and a static margin of -10%, C m curve has only one equilibrium at low AoA for full pitch down (Fig. 1d, black dotted line, δ e = 15 • ). In this case, even if the aircraft goes into deep stall, the tail command is sufficiently effective to recover at any time. This is shown with the basin of attraction which disappear for higher command than δ e = 5 • (a small basin of attraction is plotted for δ e = 5 • ). Then, the algorithm is always efficient and the time instant for applying recovery procedure has no influence on the success of the procedure. • Figs. 11c-11d. With flap down and a static margin of 0%, C m curve has a deep stall equilibrium for full pitch down (Fig. 1c, red dotted line, δ e = 15 • ). In this case, when the aircraft goes into deep stall, the tail command can be ineffective to recover. The basin of attraction exists for δ e = 15 • but is much smaller than the case where δ e = 0 • . The algorithm is always efficient. The time instant for applying recovery procedure is important because the trajectory of the aircraft in the phase portrait can cross the basin of attraction for δ e = 15 • . In this case, it seems to be adapted to apply the full pitch-down command far from this new basin. • Figs. 11e-11f. With flap down and a static margin of -10%, C m curve has a deep stall equilibrium for full pitch-down (Fig. 1d, red dotted line, δ e = 15 • ) and a large part of this curve is above the abscissa axis. In this case, when the aircraft goes into deep stall, the tail command can be highly ineffective to recover. The basin of attraction for δ e = 15 • is very large. After approximately one cycle the aircraft is in this new basin of attraction and a full pitch-down command is ineffective. The algorithm is efficient for only two cases (orange and gray curves). It can be noticed that the time instant for applying recovery procedure is necessarily at the right of the basin of attraction for δ e = 15 • . In the last two cases, it is always possible to recover from deep stall by retracting the flaps to return to the first situation (Figs. 11a-11b). But retracting flaps is not always possible in particular during hazardous situations.

VII. Conclusion

This study about the loss-of-control of an airplane in deep stall shows that a recovery could be possible in real-time using a simple algorithm for detection and recovery. To reach this goal, the first part is concerned with the analysis of deep stall in a flight dynamics point of view in order to obtain its main properties. This was done with a Learjet 23D aircraft model and non linear dynamic analysis. This analysis pointed out that the manner to reach a low or a high angle-of-attack (AoA) equilibrium is largely different. In the two cases, it is observed a short-period mode with quasi the same pulsation but with a ten times smaller damping for the high AoA equilibrium. It results a characteristically spiral trajectory in the phase portrait. Using that fact, an algorithm is proposed to reject the trajectories which do not have this form. Two criteria are then used to confirm the convergence toward a deep stall equilibrium. The first one is a comparison between α st all and the estimated value of the AoA for this equilibrium, which should be greater. The second one is a comparison between the usual short-period mode damping at low AoA and the estimated value at this moment which should be largely smaller. Using this algorithm, several cases of recovery are studied. Recovery is very effective when the pitch command is applied at the correct time instant that is when both α and α decrease. It is then proposed to fully automate the recovery procedure by applying the larger pitch-down command at the right time instant with a second very basic algorithm.

  Faure § Air Force Research Center (CReA), French Air Force Academy, Salon-de-Provence, France I. Nomenclature c = main wing chord m C L = lift aerodynamic coefficient -C D = drag aerodynamic coefficient -C m = pitching moment aerodynamic coefficient -I y = moment of inertia about the y-axis kg m

Figure 1

 1 Figure 1 Interpolated static aerodynamic model from the NASA full scale Learjet 23D wind tunnel tests for various elevator commands and two flap configurations [29] : a) C L , b) C D , c) C m s t at i c at wing quarter chord ("o" markers corresponds to NASA experimental measurements) d) C m for q = 0 and δ sm = -10%.

  7 δ T where δ T ∈ [0, 1]. The other aircraft data are: m = 4 536 kg, I y = 33940.7 kg m 2 , c = 2.14 m and S = 21.48 m 2

Figure 2

 2 Figure 2 Phase portrait comparison of the two models for flight dynamics computation of the Learjet 23D model: a) reduced model, b) full model and basin of attraction of the reduced model (flaps down δ f = 40°, δ e = 0°, δ T = 0.4 and δ sm = -10% with initial conditions V 0 = 71 m/s, γ 0 = 0°, α 0 = 0°, 5°, 10°, 13°, 20°, 23°, 30°, 45°, 48°, q 0 = 0 rad/s and h 0 = 6000 m).

Figure 3

 3 Figure 3 Phase portrait comparison of two various model for flight dynamics computation of the Learjet 23D model starting from equilibrium and applying various pitch nose-up command during 1s: a) reduced model with basin of attraction, b) full model with reduced model basin of attraction (flaps down δ f = 40°, δ T = 0.4and δ sm = -10% with initial conditions δ e i ni t = {-3, -5, -7, -9, -11, -13, -15}°, V 0 = 71 m/s, γ 0 = 0°, α 0 = 7.45°, q 0 = 0 rad/s and h 0 = 6000 m).

Figure 4

 4 Figure 4 Example of a bifurcation diagram (α, δ e ) of the Learjet 23D model (flaps up δ f = 0°, δ sm = -10%).

Figure 5

 5 Figure 5 Comparison of the basins of attraction for different configurations of the Learjet 23D model (δ sm = -10%, γ 0 = 0°, q 0 = 0 rad/s): a) flaps up for various altitude h for V = 100 m/s b) flaps down for various altitude h for V = 80 m/s, c) flaps up for various velocities V for h = 6000 m, d) flaps down with pitch damper for various altitude h for V = 80 m/s.

Figure 6

 6 Figure 6 Definition of α i and α d j for typical paths in the phase portrait at high AoA of the Learjet 23D model: a) convergence towards an equilibrium at low AoA, b) convergence towards an equilibrium at high AoA.(reduced model, flaps down δ f = 40°, δ e = 0°, δ T = 1 and δ sm = -10% with initial conditions V 0 = 72 m/s, γ 0 = 0°, α 0 = 44°, 46°, q 0 = 0 rad/s and h 0 = 10 000 m).

Figure 7

 7 Figure 7 Real-time algorithm for deep stall prediction.

Figure 8 Figure 9

 89 Figure 8 Comparison of pilot's correcting action near deep stall equilibrium: a) without algorithm detection leading to a failed recovery, b) with algorithm detection leading to a successful recovery (full model, flaps down δ f = 40°, δ e = 0°, δ T = 0.4 and δ sm = 0% with initial conditions V 0 = 100 m/s, γ 0 = 0°, α 0 = 37.5°, q 0 = 0 rad/s and h 0 = 6000 m).

Figure 10

 10 Figure 10 Real-time algorithm for recovery from deep stall.

  model Basin of attraction for � e =15°B asin of attraction for � e =0°S tep 4: neutral pitch command (� e =0°) Step 2: neutral pitch command (� e = 0°) Step 1: pitch-up command (� e =� e,init ) Step 3: full pitch-down command (� e = 15°) Basin of attraction for � e =5°F igure 11 Automatic detection and recovery of the Learjet 23D model for various configurations and static margins (γ 0 = 0°, q 0 = 0 rad.s -1 , h 0 = 6000 m): a-b) flaps up, δ sm = -10%, δ T = 0.4, V 0 = 83 m.s -1 , α 0 = 8.1°, δ e init = {-6, -6.5, . . . -10.5}°, c-d) flaps down, δ sm = 0%, δ T = 0.48, V 0 = 78.3 m.s -1 , α 0 = 4.8°, δ e init = {-12, -12.3, . . . -14.7}°, e-f) flaps down, δ sm = -10%, δ T = 0.42, V 0 = 71 m.s -1 , α 0 = 7.45°, δ e init = {-3, -5, . . . -15}°.

Table 1 Short-period mode approximation for the Learjet 23D at low and high AoA (flaps down

 1 δ f = 40°, δ sm = 0%).

	Parameters Units	Low AoA	High AoA
	α e	deg	4.4	34.2
	δ e	deg	-3	-3
	δ T	-	0.48	1
	V e	m/s	78.3	64.9
	h	km	4	4
	λ	-	-0.789 ± 1.39i -0.0762 ± 1.69i
	ω spm	rad/s	1.60	1.69
	ζ spm	-	0.49	0.045