Guillaume Novella

Christophe Macabiau

Axel Garcia-Pena

Anaïs Martineau

Pierre Ladoux

Impact of Terrestrial Emitters on Civil Aviation GNSS Receivers

Since 1994, he has been working on the application of satellite navigation techniques to civil aviation. He received his Ph.D in 1997 and has been in charge of the signal processing lab of ENAC since 2000, where he also started dealing with navigation techniques for terrestrial navigation.

INTRODUCTION

In the civil aviation field, GNSS is characterized as a Safety of Life (SoL) service. A SoL service is defined by ITU (International Telecommunication Union) radio regulations (ITU, 2020) as a radiocommunication service used permanently or temporarily for the safeguarding of human life and property. As a SoL service, GNSS must be protected from non-SoL systems that may cause degradation of GNSS performance or deny of usage. However, due to the low power of the GNSS received signal at the Earth surface, GNSS is very sensitive to radio frequency interference (RFI). Therefore, there is a major interest in analyzing GNSS RFI in order to verify that required GNSS SoL performance defined in ICAO SARPS are met.

Impact of RFI on GNSS receiver are usually analyzed in terms of carrier to noise power spectral density (𝐶/𝑁 0 ) degradation. Indeed, 𝐶/𝑁 0 reflects GNSS performance in terms of pseudo-range accuracy and service availability. In addition, a wide-band RFI is modelled as an increase of the noise floor. Therefore, effective noise 𝑁 0,𝑒𝑓𝑓 , which is the effective noise power spectral density in presence of RFI, is well adapted to characterize the impact of RFI on GNSS receiver performance.

GNSS RFI are usually divided into two categories. First, aeronautical RFI gathers all RFI sources from systems benefiting of an Aeronautical Radio Navigation System (ARNS) frequency allocation. Aeronautical RFI are often caused by aeronautical systems providing other communication services for aircraft navigation. Thus, aeronautical RFI sources are also part of SoL systems, and compatibility with GNSS must be analyzed to guarantee that both systems meet their required SoL performance. Because of this compatibility analysis and the systematic certification process when it comes to aeronautical systems, aeronautical RFI sources are completely known and monitored in terms of RF characteristics. Therefore, their impact on GNSS acquisition, tracking and demodulation performance (hereinafter referred to GNSS performance) are widely analyzed and detailed, especially in certification documents elaborated by Radio Telecommunications Committee for Aeronautics (RTCA) such as DO-235C (RTCA, 2022) and DO-292 (RTCA, 2004). Second, GNSS is also inherently impacted by spurious emission transmitted by non-aeronautical RFI sources. Non-aeronautical RFI sources gather all systems which do not have an ARNS frequency allocation, but which transmitted power in the GNSS band is not null, because of imperfect frequency isolation and imperfect filtering. Characterization of the impact of non-aeronautical RFI sources on GNSS receiver, despite its difficulty due to the high degree of uncertainty in the RFI characteristics, has a major significance when assessing GNSS receiver capability to meet minimum signal processing requirement defined in standards.

In particular, electronic devices involuntarily transmit spurious emission in the GNSS band. Therefore, these electronic devices behave as a RFI source from a GNSS receiver perspective. Code of Federal Regulations CFR Part 15 classifies electronic devices which may involuntarily radiate in other bands within two categories. Class A gathers digital devices marketed for professional usage in a commercial, industrial and business environment. Class B gathers digital devices marketed for a usage in a residential environment. Because of the high number of electronic devices involuntarily transmitting in the GNSS band, the impact analysis of these electronic devices on the GNSS receiver capacity to meet minimum signal processing requirement is of the utmost importance.

Another classification of these electronic devices can be proposed from an aircraft operation point of view. First, on-board portable electronic devices (PEDs) refer to Class B electronic devices carried on-board the aircraft by passengers. The impact of on-board PEDs on GNSS L1/E1 receivers has been analyzed in DO-235C (RTCA, 2022), and should be included in Minimum Operational Performance Standard (MOPS) DO-292A for GPS L5 and Galileo E5a receivers, adapting the L1/E1 propagation model to the L5/E5a band. Second, terrestrial emitters refer to Class A and Class B equipment located on the ground. This article focuses on the impact of terrestrial emitters on GNSS receiver capacity to meet minimum signal processing requirement. [START_REF] Nguyen | Portable Wireless LAN Device and Two-Way Radio Threat Assessment for Aircraft Navigation Radios[END_REF] analyzed radiated power of several Class A and Class B equipment. It appears that the peak radiated power is constant within the GNSS band. Therefore, terrestrial emitters impact on GNSS receiver can be characterized as an equivalent noise 𝐼 0,𝑡𝑒𝑟𝑟 , where 𝐼 0,𝑡𝑒𝑟𝑟 is the power level in a 1 MHz bandwidth at the receiver antenna port of the aggregate of all terrestrial sources RFI. DO-235B (RTCA, 2008) made a first attempt in estimating 𝐼 0,𝑡𝑒𝑟𝑟 for GPS L1C/A legacy receiver. This estimation was performed near Los Angeles airport. However, very unprecise assumptions on the density of emitters were considered, leading to pessimistic estimation of 𝐼 0,𝑡𝑒𝑟𝑟 . In addition, free-space losses model was considered in (RTCA, 2008) for all phases of flight. However, free-space loss model may not be the most adequate propagation model to characterize radio propagation in the ground vicinity during approach for examples, as this model does not take into account multipath and the consequent potential fading for example. [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] proposed an alternative propagation model for phases of flight below 2500 ft above mean sea level (AMSL). This alternative propagation model was adopted in DO-235C and is based on attenuation empirical data in urban environment presented in [START_REF] Erceg | An empirically based path loss model for wireless channels in suburban environments[END_REF] and in [START_REF] Hata | Empirical formula for propagation loss in land mobile radio services[END_REF]. ITU also developed its own propagation model for ultra-high frequency (UHF) communications. This ITU model is called P-528 and is presented in (ITU, 2021). This article has three main objectives. First, it proposes a refinement of the estimation of 𝐼 0,𝑡𝑒𝑟𝑟 in the L1/E1 frequency band. As the density of terrestrial emitters is linked to the density of population, a more precise estimation of the density of emitters is proposed from European density of population dataset [START_REF] Gallego | A population density grid of the European Union[END_REF]. Second, it proposes a comparison of 𝐼 0,𝑡𝑒𝑟𝑟 in the L1/E1 frequency band obtained with three propagation channel models: free space loss, DO-235C model and ITU P-528. Third, the estimation of 𝐼 0,𝑡𝑒𝑟𝑟 in the L5/E5a band is computed.

To achieve these objectives, this article is divided in three main sections. First, the mathematical approach to estimate 𝐼 0,𝑡𝑒𝑟𝑟 is provided. Second, the main inputs to compute 𝐼 0,𝑡𝑒𝑟𝑟 are described. One of these inputs is the propagation channel model. The existing propagation channel models used for analyzing ground emitters impact on GNSS receiver noise floor are reviewed and ITU-P528 model is introduced. Moreover, inputs include the density of terrestrial emitters distribution, and a more precise estimation of the density of emitters compared to the estimation done in (RTCA, 2008) is proposed. Third, 𝐼 0,𝑡𝑒𝑟𝑟 is computed for different phases of flight in the L1/E1 and L5/E5a frequency bands.

TERRESTRIAL EMITTERS EQUIVALENT NOISE

This section proposes a mathematical model to derive terrestrial emitters equivalent noise. The proposed mathematical formula to derive 𝐼 0,𝑡𝑒𝑟𝑟 is based on the calculation model of [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF].

The power 𝑃 𝑡𝑒𝑟𝑟 from terrestrial emitters received at the GNSS receiver antenna port is modeled as a random variable. Let us define 𝑁 the number of active terrestrial emitters visible by the aircraft GNSS antenna. In [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF], 𝑁 follows a Poisson distribution with mean value 𝑁 given by Equation (1).

𝑁 = ∫ ∫ 𝑑(𝑟, 𝜑)𝑑𝑟𝑑𝜑 2𝜋 0 𝑅 𝐿𝑂𝑆 0 (1)
𝑑(𝑟, 𝜑) is the density of terrestrial emitters at the position (𝑟, 𝜑) and 𝑅 𝐿𝑂𝑆 is the radio line of sight horizon. (𝑟, 𝜑) are the polar coordinates on a horizontal frame centered on the aircraft.

Let us define (𝑋 𝑘 ) 𝑘∈⟦1;𝑁⟧ the contribution of each individual terrestrial emitters to 𝑃 𝑡𝑒𝑟𝑟 . The position of each terrestrial emitter, defined in polar coordinates (𝑟 𝑘 , 𝜑 𝑘 ), is also assumed to be random. The family of random variables (𝑟 𝑘 , 𝜑 𝑘 ) 𝑘∈⟦1,𝑁⟧ is independent and identically distributed. The random position (𝑟 𝑘 , 𝜑 𝑘 ) of an emitter thus only depends on the density of terrestrial emitters and do not depend on the random number of emitters in view by the GNSS receiver, 𝑁. The probability density function of the position of a given terrestrial emitter, 𝑓, (the probability of finding a terrestrial emitter at a given position) depends on the density of terrestrial sources at the corresponding position 𝑑(𝑟, 𝜑). The probability density function of the position of a given terrestrial emitter is given by Equation (2).

𝑓(𝑟

, 𝜑) = 𝑑(𝑟, 𝜑) ∫ ∫ 𝑑(𝑟, 𝜑)𝑑𝑟𝑑𝜑 2𝜋 0 𝑅 𝐿𝑂𝑆 0 = 𝑑(𝑟, 𝜑) 𝑁 (2)
Let us define 𝑋 𝑘 (𝑟 𝑘 , 𝜑 𝑘 ) as the received power in a 1 MHz bandwidth at the receiver antenna port from one terrestrial emitter source located at position (𝑟 𝑘 , 𝜑 𝑘 ). The contribution 𝑋 𝑘 (𝑟 𝑘 , 𝜑 𝑘 ) of each terrestrial emitter to 𝑃 𝑡𝑒𝑟𝑟 depends on three parameters: -The equivalent isotropic radiated power 𝑃 0 in a 1 MHz bandwidth by the RFI source. As discussed in the next section, this term is assumed to be a representative value for all terrestrial emitters. -The receiver antenna gain 𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 ).

-Propagation losses 𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 ).

The expression of 𝑋 𝑘 (𝑟 𝑘 , 𝜑 𝑘 ) as a function of 𝑃 0 , 𝐺 𝑅𝑋 and 𝐿 𝑝𝑟𝑜𝑝 is given in Equation (3).

𝑋 𝑘 (𝑟 𝑘 , 𝜑 𝑘 ) = 𝑃 0 𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 )𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 )

(3)

The received power in a 1 MHz bandwidth at the antenna port from terrestrial emitters is then given by Equation [START_REF] Rtca | DO 292-Assessment of Radio Frequency Interference Relevant to the GNSS L5/E5A Frequency Band[END_REF].

𝑃 𝑡𝑒𝑟𝑟 = ∑ 𝑋 𝑘 (𝑟 𝑘 , 𝜑 𝑘 ) 𝑁 𝑘=1 (4)
The equivalent noise induced by terrestrial emitters on GNSS receiver is then defined by the mean value of 𝑃 𝑡𝑒𝑟𝑟 . The calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is performed in Equation ( 5). (5)

Note that considering the average received power to define 𝐼 0,𝑡𝑒𝑟𝑟 may sound unsafe, since the received power would regularly exceed the mean value. In order to avoid an underestimation of the of the impact of terrestrial emitters on equivalent noise, a 6 dB safety margin is usually added on the estimated non-aeronautical RFI equivalent noise. The goal of this safety margin is to cover for uncertainties in the calculation of non-aeronautical RFI equivalent noise, and to cope with the potential exceedance of the nonaeronautical RFI received power compared to the average received power. [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] showed that the probability that the terrestrial emitter received power 𝑃 𝑡𝑒𝑟𝑟 exceeds 𝐼 0,𝑡𝑒𝑟𝑟 plus the 6 dB safety margin is very low (around 3 10 -4 ).

In order to precisely estimate 𝐼 0,𝑡𝑒𝑟𝑟 , an accurate knowledge of 𝑃 0 , 𝐿 𝑝𝑟𝑜𝑝𝑎 and 𝑑 is needed. The modeling of these parameters is discussed in the following sections. This modeling is a key point to obtain a realistic estimation of 𝐼 0,𝑡𝑒𝑟𝑟 .

INPUT PARAMETERS FOR 𝑰 𝟎,𝒕𝒆𝒓𝒓 CALCULATION

Each of the main input parameters allowing to derive 𝐼 0,𝑡𝑒𝑟𝑟 are now discussed from a civil aviation point of view in a dedicated section. As highlighted in Equation ( 5), main contributors to 𝐼 0,𝑡𝑒𝑟𝑟 include the radiated EIRP by one ground source 𝑃 0 , the receiver antenna gain 𝐺 𝑅𝑋 , propagation losses 𝐿 𝑝𝑟𝑜𝑝 (determined by the propagation channel model) and the density of terrestrial emitters 𝑑. This section discusses these parameters from a civil aviation point of view.

Transmitted EIRP 𝑷 𝟎 by each terrestrial source

For simplicity purposes, it is proposed to select a representative value for 𝑃 0 which would be applied to any individual terrestrial emitter. Transmitted emissions of electronic devices in the GNSS band are constrained by 47 CFR Part 15 (Code of Federal Regulations, 1989). In particular, the field strength of unintentional radiators such as terrestrial emitters, measured at a distance of 3 m from the radiator, shall not exceed 500 µV/m at frequencies higher than 960 MHz. The transmitted EIRP 𝑃 0 is linked to the field strength through the commonly used relation between the field strength and the power on a sphere centered on the radiating source given by Equation ( 6).

𝑃 0 4𝜋𝑑 𝑟 2 = 𝐸 2 120𝜋 (6)
𝐸 is the electric field measured at a distance 𝑑 𝑟 from the radiating source. From 47 CFR Part 15 constraint on the radiated electric field, the maximum EIRP in 1 MHz bandwidth value transmitted by a given terrestrial emitter is thus equal to 𝑃 0 = -71.25 𝑑𝐵𝑊/𝑀𝐻𝑧.

NASA (National Aeronautics and Space Administration) performed measurement on the involuntarily transmitted power by a wide diversity of electronic devices: laptops, WLAN (wireless local area network) devices, Bluetooth devices, FRS (family radio service) and GMRS (general mobile radio service) radios, etc. The results of this measurement campaign are presented in [START_REF] Nguyen | Portable Wireless LAN Device and Two-Way Radio Threat Assessment for Aircraft Navigation Radios[END_REF]. From these measurements, it appears that the transmitted power in the GNSS band is much lower than what is allowed by 47 CFR Part 15. Therefore, [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] decided to use the value 𝑃 0 = -81.1 𝑑𝐵𝑊/𝑀𝐻𝑧 as a representative value of EIRP radiated by terrestrial emitters. For consistency with the previous analysis, this value is adopted in this article.

Receiver antenna gain 𝑮 𝑹𝑿

Civil aviation GNSS receiver antennas are standardized in DO-301 (RTCA, 2006) for GPS L1C/A legacy receiver and DO-373 (RTCA, 2018) for Dual Frequency Multi Constellation (DFMC) GNSS receiver. In particular, the GNSS receiver antenna gain shall not exceed a maximum value represented in Figure 1 as a function of the elevation angle of the incoming signal in the lower hemisphere of the antenna. As an important remark, RFI from terrestrial emitters usually arrive with a negative elevation angle. On the one hand, and since the GNSS antenna is installed on the top of the aircraft, the aircraft fuselage may shadow the signal. On the other hand, the fuselage may also act as a waveguide for the terrestrial emitter RFI signal. In any case, the impact of the fuselage for signals arriving with a negative elevation angle is already taken into account in the antenna gain plots of Figure 1, since the maximum antenna gain at lower hemisphere has been defined from measurements performed on manufactured antennas installed on an aircraft model (see Appendix G of (RTCA, 2008)). A worst case with no polarization losses is considered here.

The constraint on the maximum antenna gain at the lower hemisphere helps in limiting the received RFI power, as most RFI come from below the aircraft. DO-373 antenna is more constraining with a maximum antenna gain of -6 dB for incoming signals with a 0° elevation, against a maximum antenna gain of -5 dB for the DO-301 antenna. In addition, a higher constraint is imposed on antennas installed on aircraft certified for Cat II and Cat III approaches. In this case, the maximum antenna gain is lowered to -13 dB for signals with elevation angles between -90° and -45°. 𝐸(𝑟, 𝜑) = -atan ( 𝐻 𝑎 -ℎ 𝐵 𝑟 ) (7)

Propagation channel model

Propagation channel model is a key parameter impacting the RFI received power. In this application, the targeted propagation channel models are narrowband propagation channels which as are thus completely defined, from a RFI impact point of view, from the propagation losses term. Moreover, since the worst case is searched for, only the statistical propagation losses are inspected without considering the time-evolution.

Initial estimation of 𝐼 0,𝑡𝑒𝑟𝑟 in L1/E1 band performed in (RTCA, 2008) was based on the free space loss model. A refinement of the propagation channel model was later adopted in DO-235C. This new propagation channel model, that is referred as DO-235C propagation channel model in this article, is in fact composed by several empirical models depending on the horizontal distance, 𝑟, between the RFI source and the aircraft. Eventually, ITU developed its own model for UHF propagation channel which is called ITU P-528. This model is particularly adapted for computing propagation losses between a ground source and the aircraft.

The objective of this section is thus to present these three propagation channel models and to compare the propagation losses introduced by each of them as a function of the horizontal distance 𝑟.

Free Space Loss Model

Free space loss propagation channel model was used to estimate 𝐼 0,𝑡𝑒𝑟𝑟 in DO-235B. The propagation losses between the terrestrial RFI source and the aircraft antenna are simply deduced from the distance between the RFI source and the aircraft GNSS receiver antenna. Free space losses mathematical model is given in Equation ( 8).

Height of terrestrial emitter ℎ 𝐵

Terrestrial emitter 𝐸(𝑟, 𝜑)

Height of aircraft antenna 𝐻 𝑎 r 𝐿 𝑝𝑟𝑜𝑝𝑎 = ( 𝜆 4𝜋√𝑟 2 + (𝐻 𝑎 -ℎ 𝐵 ) 2 ) 2 (8)
𝜆 is the wavelength of the GNSS signal center frequency.

Free space loss propagation channel model has the drawback of not considering propagation effects near the ground, such as reflection and diffraction which could potentially lead to fading. Thus, this model is adapted to compute 𝐼 0,𝑡𝑒𝑟𝑟 when the aircraft is at high altitude, but may be not suited for low altitude operations.

DO-235C propagation model

DO-235C propagation model considers that the visible surface area below the aircraft is split in three zones as represented in Figure 3.

Figure 3 Illustration of DO-235C situation for DO-235C propagation model

The selected propagation channel and the propagation losses of each one of the three zones are explained next. First, the selected propagation channel and propagation losses will be described without considering fading. Afterwards, the propagation losses determined for each zone will be modified to consider the fading phenomenon; where the fading phenomenon will be characterized with the same model for each zone but with different numerical values to express the different impact on the different zones. Finally, the process to determine the different regions, or equivalently to determine 𝑅 1 and 𝑅 2 , is described.

Zone 1:

The first zone corresponds to locations where the horizontal distance between the emitter and the aircraft, 𝑟, is lower than a radius R1 (𝑅 1 is addressed later on this article). In that zone, potential reflection of the interfering signal is taken into account using a tworays model (combination of one direct ray plus its multipath) which details are presented next.

The direct and reflected ray path lengths are given by Equation ( 9).

𝑅 𝑑𝑖𝑟 (𝑟) = √(𝐻 𝑎 -ℎ 𝐵 ) 2 + 𝑟 2 𝑅 𝑟𝑒𝑓𝑙 (𝑟) = √(𝐻 𝑎 + ℎ 𝐵 ) 2 + 𝑟 2 (9)
The phase offset between the direct on reflective rays is given by Equation (10).

𝜙(𝑟) = 2𝜋 𝜆 (𝑅 𝑑𝑖𝑟 (𝑟) -𝑅 𝑟𝑒𝑓𝑙 (𝑟)) (10) Zone 1 Zone 2 Zone 3 𝑅 𝐿𝑂𝑆 𝑅 2 𝑅 1
The path loss under two rays model also depends on the amplitude ratio of the reflected signal compared to the direct ray. This amplitude is linked to the ground characteristics described by the relative dielectric constant 𝜀 𝑟 and conductivity 𝜎 𝑐𝑐 . The complex reflection coefficient parameter 𝜌 𝑣 for vertical polarized signal is given by Equation ( 11).

𝜌 𝑣 (𝑟) = (𝜀 𝑟 -𝑖𝑥) sin(𝜃(𝑟)) -√(𝜀 𝑟 -𝑖𝑥) -cos 2 (𝜃(𝑟)) (𝜀 𝑟 -𝑖𝑥) sin(𝜃(𝑟)) + √(𝜀 𝑟 -𝑖𝑥) -cos 2 (𝜃(𝑟)) 𝜃(𝑟) = asin ( 𝐻 𝑎 + ℎ 𝐵 𝑅 𝑟𝑒𝑓𝑙 (𝑟) ) 𝑥 = 𝜎 𝑐𝑐 2𝜋𝑓 𝑐 𝜀 0 (11)
𝑓 𝑐 is the GNSS carrier frequency and 𝜀 0 is free space permittivity. 𝜃(𝑟) is called grazing angle. 𝑖 = √-1.

Eventually, two-rays propagation loss, not considering fast fading, is expressed as in Equation ( 12).

𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 = ( 𝜆 4𝜋 |1 + 𝑅 𝑑𝑖𝑟 (𝑟) 𝑅 𝑟𝑒𝑓𝑙 (𝑟) 𝜌 𝑣 (𝑟)𝑒 -𝑖𝜙(𝑟) | 𝑅 𝑑𝑖𝑟 (𝑟) ) 2 (12)
Representative values chosen in DO-235C and in [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] for ground dielectric constant and conductivity are respectively 𝜀 𝑟 = 7 and 𝜎 𝑐𝑐 = 0.15 𝑆/𝑚. These values are also used in this article.

Zone 2:

The second zone represents the locations where the horizontal distance between the emitter and the aircraft, 𝑟, is between R1 and R2 (𝑅 2 is addressed later on this article). In this zone, the applied propagation channel model depends on the altitude of the aircraft, 𝐻 𝑎 . If the aircraft is below 80m (𝐻 𝑎 < 80), the Erceg-Greenstein model is adopted [START_REF] Erceg | An empirically based path loss model for wireless channels in suburban environments[END_REF]. This model provides inputs to estimate propagation losses if the transmitter is in a flat surface with a light tree density. An airport environment corresponds to that description. If the aircraft is above 80m (𝐻 𝑎 > 80) or if the horizontal distance between the aircraft and the terrestrial emitter is above 8 km, a simple log-log constant slope model is adopted since the Erceg model is no longer applicable.

Erceg propagation losses, if fast fading is not considered, can be expressed as in Equation ( 13).

𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 = ( 𝜆 4𝜋𝑟 0 ) 2 ( 𝑟 0 𝑟 ) 𝑎-𝑏𝐻 𝑎 + 𝑐 𝐻 𝑎 (13)
Where 𝑟 0 = 100 and coefficients 𝑎, 𝑏 and 𝑐 have been empirically determined from the power received by 95 cellular base stations in a wide diversity of terrains. In particular, for a flat area with a light tree density such as an airport, these coefficients are summarized in below in TABLE 1. The third zone corresponds to locations where the horizontal distance between the emitter and the aircraft, 𝑟, is beyond radius 𝑅 2 until the Radio Line-of-Sight radius (RLOS). In this zone, the Hata-Okumura propagation model described in [START_REF] Hata | Empirical formula for propagation loss in land mobile radio services[END_REF]

Where 𝑎(ℎ 𝐵 ) is a correction term to take into account the height of the terrestrial RFI source. In urban environment and for frequencies above 400 MHz, 𝑎(ℎ 𝐵 ) is given by Equation ( 16).

𝑎(ℎ 𝐵 ) = 3.2(log(11.75ℎ 𝐵 )) 2 -4.97 (16)
Injection Equations ( 15) and ( 16) into Equation ( 14), the Hata-Okumura propagation losses can be re-expressed by Equation ( 17). In this case, 𝑅 1 and 𝑅 2 are computed as follows. 𝑅 1 is determined to target a grazing angle verifying tan(𝜃(𝑅 1 )) ≈ 0.5. Some adjustments can be done to avoid strong propagation loss slope breakout. 𝑅 2 is determined such that the elevation angle toward the antenna is around 4°: 𝐸(𝑅 2 ) ≈ -4°, where 𝐸(𝑅 2 ) is computed using Equation (7).

𝐿 𝑝𝑟𝑜𝑝𝑎

ITU P-528 propagation model

ITU P528 model is the method recommended by the International Telecommunications Union to evaluate propagation losses for aeronautical systems transmitting in the 100 MHz -30 GHz frequency range. It can be applied for transmitter and receiver antennas height between 1.5 m and 20,000 m. The objective of ITU P528 model is to be more accurate than free space loss model. Indeed, ITU P528 model includes additional propagation considerations. First, attenuation due to atmospheric absorption is added to free space losses. Second, two-rays model is implemented to account for reflection on the ground of the interfering signal. Third, smooth Earth diffraction losses are also considered. Fourth, ITU P528 model also includes transhorizon propagation, considering reflection on the troposphere. Note that this troposcatter component can be neglected for low power interference sources such as terrestrial emitters but may be not negligeable when it comes to high power interference such as DME/TACAN or jammer. Fifth, ITU P528 ray tracing uses ITU R P676-12 model (ITU, 2019) instead of the traditional "4/3 Earth" method. The main difference is the consideration of the variation of the atmosphere index with altitude. However, ITU R P676-12 model is expected to be very close to the "4/3 Earth" method when the aircraft receiver is at low altitude.

The method to compute propagation losses is widely detailed in (ITU, 2021), and thus it is not developed in this article. In this article, the p parameter of the P-528 model, which refers to the time probability that propagation losses exceed the returned value, is set to 50%. In addition, a C software providing propagation losses under P-528 model is given by ITU. This software has been used to obtain P-528 propagation loss results in this article.

𝐻 𝑎 -ℎ 𝐵 𝑅 1

Density of emitters function

Density of emitters function, 𝑑, has a key role in the received power in a 1 MHz bandwidth, 𝑃 𝑡𝑒𝑟𝑟 , where the number of transmitting terrestrial sources, 𝑁, as well as their position, (𝑟, 𝜑), are linked to the density of emitters function. Several assumptions were considered in the past to estimate the density of emitters function. However, these hypotheses seemed quite conservative leading to an over estimation of 𝐼 0,𝑡𝑒𝑟𝑟 . In this section, first past assumptions and models concerning the derivation of the density of terrestrial emitters function are reminded, mainly the DO-235B models for low altitude and high-altitude phases of flight. Second, new density of emitters function models are proposed, one for intermediate/high-altitude phase of flight (for example en-route) based on population density data and another one for low altitude phases of flight (mainly take-off, approach or taxiing) specifically customized for an airport environment. Note that in this section the mathematical model defining the density of emitters as a function of (𝑟, 𝜑) is referred to as density of emitters function, 𝑑(𝑟, 𝜑), whereas specific numerical values associated to a given position (𝑟 = 𝑟 𝑝𝑜𝑠 , 𝜑 = 𝜑 𝑝𝑜𝑠 ) or to a given zone (𝑟 𝑧𝑜𝑛𝑒_𝑚𝑖𝑛 ≤ 𝑟 ≤ 𝑟 𝑧𝑜𝑛𝑒_𝑚𝑎𝑥 , 𝜑 𝑧𝑜𝑛𝑒_𝑚𝑖𝑛 ≤ 𝜑 ≤ 𝜑 𝑧𝑜𝑛𝑒_𝑚𝑎𝑥 ) are referred to as density of emitters.

Previous density of emitters models -DO-235B

DO-235B en-route (high-altitude) density of emitters function In DO-235B, the density of emitters function, 𝑑(𝑟, 𝜑), during en-route operation is defined as follows. The aircraft is supposed to fly above a metropolitan area, since it is expected that the highest 𝐼 0,𝑡𝑒𝑟𝑟 is obtained in this situation. Density of emitters function is defined identifying three circular zones centered at the aircraft position as illustrated in Figure 5. . , 2008) and [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] consider that a reasonable estimate for the density of emitters in a urban area is 10 -4 𝑚 -2 . Indeed, [START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] mentions that cellphones activity in Los Angeles does not exceed 0.78 active emitter per hectare. Therefore, 10 -4 active terrestrial emitter per square meter seems to be a reasonable upper bound. Zone A radius is set to 30 km in DO-235B. Zones B and C respectively refers to suburban and rural areas. Density of emitters in suburban areas is set to The main drawback of this configuration is the lack of flexibility to the environment configuration. Indeed, this configuration does not take into account the number and size of the cities in view from the aircraft. Therefore, the method presented in this section provides an indicative value of 𝐼 0,𝑡𝑒𝑟𝑟 for en-route operations but cannot be translated to other situations where the environment is different. The 𝐼 0,𝑡𝑒𝑟𝑟 estimation using this model of terrestrial density of emitters is thus a rough estimation as this configuration does not match the actual analyzed scenario. Radius 𝑅 1 and 𝑅 2 values can be modified to make the estimation more or less conservative; however, 𝐼 0,𝑡𝑒𝑟𝑟 estimation is still loose.

DO-235B low and intermediate altitude density of emitters function

In DO-235B, for low altitude operations analysis, the density of emitters function is derived differently from the en-route case. Low altitude operations usually take place in the vicinity of an airport, as these operations refers to take-off or approach phases of flight.

The method used in DO-235B to estimate the density of emitters function to compute 𝐼 0,𝑡𝑒𝑟𝑟 for low altitude operations consists in identifying zones having a different density of emitters depending on their environment, as shown in Figure 6.

Figure 6 Terrestrial emitters environment in DO-235B low-altitude 𝐼 0,𝑡𝑒𝑟𝑟 calculation

The environment is divided into four categories of zones.

𝑅 𝐿𝑂𝑆

-Zone A corresponds to the nominal environment. Zone A gathers all rural and suburban areas. The density of terrestrial emitters is set to 1 3 10 -4 𝑚 -2 in Zone A.

-Zone B corresponds to urban areas. It gathers all cities visible by the aircraft GNSS antenna. The density of terrestrial emitters is set to 10 -4 𝑚 -2 in Zone B. -Zone C corresponds to a zone in which a high number of terrestrial emitters are used simultaneously. For example, DO-235B identified airport terminals as a terrestrial emitters hot-spot. In this zone, the density of terrestrial emitters is chosen to be ten times higher than the density of terrestrial emitters in Zone B (urban environment). Therefore, the density of terrestrial emitters in Zone C is set to 10 -3 𝑚 -2 . -Zone D is a no-emission zone. It represents areas in which no terrestrial emitter is active. For example, oceans are part of Zone D. Moreover, airport movement areas are part of Zone D as no emitters are expected to be active in this area. Density of terrestrial emitter is 0 in zone D.

The density of terrestrial emitter on each zone is given in TABLE 3. As illustrated in Figure 6, Zones B, C and D are defined by the inner and the outer radius, ((𝑟 𝑧𝑜𝑛𝑒_𝑚𝑖𝑛 ≤ 𝑟 ≤ 𝑟 𝑧𝑜𝑛𝑒_𝑚𝑎𝑥 ), and an aperture angle (𝜑 𝑧𝑜𝑛𝑒_𝑚𝑖𝑛 ≤ 𝜑 ≤ 𝜑 𝑧𝑜𝑛𝑒_𝑚𝑎𝑥 ). This density of emitters function, 𝑑(𝑟, 𝜑) has the advantage of being more precise than DO-235B en-route model and to be adaptable to each particular situation. The drawback of this model is that the definition of the border between the different zones remains unprecise. Moreover, another weakness of this model is the assumption that all cities have the same properties in terms of concentration of active emitters (or equivalently the same density of users numerical value).

New proposed density of emitters models

Intermediate/high-altitude density of emitters function defined from population density In order to make the density of terrestrial emitters function more accurate, it is proposed to link the density of active terrestrial emitters (numerical value) to the population density. The European density of population which is used in this article is taken from [START_REF] Gallego | A population density grid of the European Union[END_REF]. Figure 7 represents the population density in inhabitants per square kilometer across western Europe, with a resolution of 1 km*1 km (translated into latitude and longitude coordinates). Number of inhabitants per square kilometer is clipped to 1000 in Figure 7 for illustration purpose, even though the maximum population density from the population density dataset 40,000 inh/km².

Figure 8 Density of mobile phone users as a function of population density for eight Spanish cities

Figure 8 shows that the number of simultaneously active terrestrial emitters (or density of emitters), 𝑑, can be linearly linked to the population density, 𝑑 𝑝𝑜𝑝 . In addition, [START_REF] Nguyen | Portable Wireless LAN Device and Two-Way Radio Threat Assessment for Aircraft Navigation Radios[END_REF] shows that electronic devices such as mobile phone are part of the most powerful terrestrial emitters category. Moreover, density of mobile phone is expected to be higher than the density of other terrestrial emitters. Consequently, it is proposed to define the density of terrestrial emitter function as shown in Equation ( 21). The advantage of this method with respect to previous DO-235B en-route method is its capacity to precisely reflect the terrestrial RFI environment, with a good resolution (data set allow a resolution up to 100m * 100m, in this article a resolution of 1 km * 1 km is sufficient to estimate 𝐼 0,𝑡𝑒𝑟𝑟 at intermediate and high altitude). For example, it takes advantage of the real density of population in rural non-inhabited areas to compute 𝐼 0,𝑡𝑒𝑟𝑟 whereas DO-235B clipped the density of emitter to a non-null value. Therefore, the distribution of emitters function described in this paragraph is expected to give a lower and more precise 𝐼 0,𝑡𝑒𝑟𝑟 value in comparison to terrestrial emitters distribution of DO-235B.

𝑑 (𝑟

Low altitude (airport environment) density of emitters function

The precision on the calculation of 𝐼 0,𝑡𝑒𝑟𝑟 on low altitude and ground operations depends on the precision of the terrestrial emitters distribution in the airport environment. However, the usage of population density map is not relevant in an airport since the dataset does not have a sufficient resolution and thus to use the population density dataset may lead to under-estimation of 𝐼 0,𝑡𝑒𝑟𝑟 .

Airport is a particular environment in terms of terrestrial emitters distribution. On one hand, just few terrestrial emitters are active on movement area which includes runway and taxiway. Indeed, electronic devices carried on-board by passengers are switched off during taxiing, take-off and landing phases. On the other hand, there is a high density of emitters inside the public facilities of the airport. Therefore, the airport environment and the associated density of terrestrial emitters function, 𝑑 𝑎𝑖𝑟𝑝𝑜𝑟𝑡 (𝑟, 𝜑), is proposed to be modelled as follows: 1) Airport movement area is a no-emission zone 2) Terrestrial emitters are only located inside buildings within the airport area; density of terrestrial emitters is set to 10 -4 𝑚 -2 inside airport buildings The density of terrestrial emitter is also set to 10 -4 𝑚 -2 outside the airport area since the airport is located in urban area. . Figure 10 illustrates the density of emitters near two airports: Frankfurt EDDF and Los Angeles KLAX. As the aircraft is inside the airport no-emission zone when it is taxiing or in its final approach phases, the calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is expected to be more precise using this terrestrial emitters' distribution model in comparison to a uniform density of terrestrial emitters assumption which is considered in DO-235C.

RESULTS: COMPUTATION OF 𝑰 𝟎,𝒕𝒆𝒓𝒓

This section derives 𝐼 0,𝑡𝑒𝑟𝑟 for L1/E1 and L5/E5a receivers, using the different propagation channel models and density of terrestrial emitters function models presented respectively in sections 3.3 and 3.4. This section is divided in three parts. First, locations and altitude of interest pairs, at which 𝐼 0,𝑡𝑒𝑟𝑟 is derived, are identified. Moreover, for each identified case, the most relevant propagation channel model and density of terrestrial emitters function model are chosen. Second, a comparison between the three propagation channel models is performed. Third, impact of terrestrial emitters on civil aviation GNSS receiver is deduced computing 𝐼 0,𝑡𝑒𝑟𝑟 for each identified points of interest.

Location and altitude of interest

Impact of terrestrial emitters on GNSS receivers is analyzed at some key points of interest as part of the RFI environment analysis standardization process (DO-235B/C, DO-292A). The objective of RFI environment analysis standardization process is to ensure the capability of the receiver to meet GNSS minimum performance objectives defined in ICAO Annex 10 (ICAO, 2006). As part of the demonstration of GNSS receiver capability to fulfill GNSS minimum performance objectives, 𝐼 0,𝑡𝑒𝑟𝑟 is computed at several altitudes corresponding to different phases of flight:

-Taxiing: The aircraft is located on the ground at an airport gate. The height of the GNSS receiver antenna, which is located on the top of the aircraft, is 𝐻 𝑎 = 5.64 𝑚 (stabilizer height of a Boeing B-737). 𝐼 0,𝑡𝑒𝑟𝑟 is computed at two representative international airports: Frankfurt EDDF in Europe and Los Angeles KLAX in the United States (US). Propagation channel models analyzed here are the free space loss, P-528 and DO-235C models. In addition, the density of terrestrial emitters function of the airport environment, 𝑑 𝑎𝑖𝑟𝑝𝑜𝑟𝑡 (𝑟, 𝜑), is selected for precision purposes. -Cat II decision height (DH): The minimum Cat II precision approach decision height is 100 ft. As in (RTCA, 2022), it is here considered that the aircraft is 15 ft below the nominal approach slope, so 𝐻 𝑎 = 25.94 𝑚. Free space, P528 and DO-235C propagation channel models are compared. The calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is done considering the density of terrestrial emitters function of the airport environment, 𝑑 𝑎𝑖𝑟𝑝𝑜𝑟𝑡 (𝑟, 𝜑), of EDDF and KLAX. For this case, the Cat II/III antenna model is used. -Cat I DH: The minimum Cat I precision approach DH is 200 ft. As in (RTCA, 2022), it is considered that the aircraft is 25 ft below the nominal approach slope, so 𝐻 𝑎 = 53.34 𝑚. Free space, P528 and DO-235C propagation channel models are compared. The calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is done considering the density of terrestrial emitters function of the airport environment, 𝑑 𝑎𝑖𝑟𝑝𝑜𝑟𝑡 (𝑟, 𝜑), of EDDF and KLAX. -Intermediate altitude: Intermediate altitude scenario corresponds to the situation where the aircraft is at Final Approach Fix height point. As population density data is only available in Europe, 𝐼 0,𝑡𝑒𝑟𝑟 is only computed at EDDF FAF. From EDDF Instrument Approach Chart (IAC), TIXAK is the FAF with the lowest altitude among EDDF published RNAV (Radio navigation) procedures. The altitude of TIXAK is 2000 ft AMSL, corresponding to 520 m above ground level (AGL). , Free space and propagation channels models are used whereas DO-235C propagation model is no applicable at this altitude. Calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is done considering terrestrial emitters density derived from using the intermediate/high altitude density of emitters function defined from population density, 𝑑 𝑖𝑛𝑡/ℎ𝑖𝑔ℎ (𝑟, 𝜑). -En-route: The identified European aeronautical RFI hot-spot for L5/E5a GNSS receivers is at high altitude (flight level FL 400), because of a high number of RF in-view DME (Distance Measurement Equipment) and TACAN (Tactical Air Navigation) ground stations, which transmit high power pulsed aeronautical RFI. European aeronautical RFI hot spot is located in 50.5°N, 5.5°E. Calculation of 𝐼 0,𝑡𝑒𝑟𝑟 is done using free-space and P-528 propagation channel models as well as intermediate/high altitude density of emitters function defined from population density, 𝑑 𝑖𝑛𝑡/ℎ𝑖𝑔ℎ (𝑟, 𝜑)..

Figure 11

Propagation losses for DO-235C, P-528 and free-space models in L5/E5 band for taxiing, Cat II and Cat I operations Some conclusions can be deduced from Figure 11. 1) Free-space and P-528 losses are equal as long as 𝑟 remain below a particular breakpoint, 𝑟 𝑏 . Above this breakpoint, 𝑟 𝑏 , atmospheric loss, Earth diffraction and troposcatter effects become non-negligeable, and thus, P-528 losses become more important than free space losses. 2) DO-235C propagation losses are higher than free space and P-528 losses as long as 𝑟 > 𝑅 1 . As a reminder, P-528 model is the recommended ITU propagation channel model for aeronautical for air-to-air, air-to-ground or ground-to-air UHF communications. 3) When 0 ≤ 𝑟 < 𝑅 1 , DO-235C model presents some oscillations. Indeed, when 𝑟 ≤ 𝑅 1 , DO-235C propagation uses a tworays model, and the multipath is successively constructive and destructive depending on 𝑟. When 𝑟 tends toward 0, the free space and P-528 propagation losses are higher than DO-235C propagation losses by 3dB.

Eventually, Figure 12 shows P-528 and free space losses for intermediate and high-altitude operations.

Figure 12

Propagation losses for P-528 and free space models for intermediate and high-altitude operations

Figure 12 shows a similar relative behavior between P-528 and free space loss models at intermediate and high altitudes compared to Figure 11. Free space and P-528 losses are very similar up to a breakpoint. As P-528 losses are then higher than free space losses, 𝐼 0,𝑡𝑒𝑟𝑟 is expected to be lower using P-528 model than free space losses.

4.3

Results: 𝑰 𝟎,𝒕𝒆𝒓𝒓 calculation This section provides the values of 𝐼 0,𝑡𝑒𝑟𝑟 for the different cases identified in TABLE 4. For each operation (taxiing, Cat I and Cat II precision approach, intermediate and high altitude operations), 𝐼 0,𝑡𝑒𝑟𝑟 is computed. Moreover, this section discusses and highlights the conclusions of the results on the impact of terrestrial emitters on GNSS receivers.

Taxiing TABLE 6 provides the value of 𝐼 0,𝑡𝑒𝑟𝑟 for taxiing operations. The results are given for both the L1/E1 band and L5/E5a band, in Frankfurt EDDF and Los Angeles KLAX. In both cases, the DFMC Cat II antenna model is used, and despite not being given here, results with DMFC Cat I antenna model are very similar. The aircraft is supposed to be at an airport gate. -156.5 -156.8 -148.5 -154.4 -154.8 -150.8 -158.6 -161.1 -148.3 -156.9 -158 𝐼 0,𝑡𝑒𝑟𝑟 estimated with free space losses is higher by at least 5.5 dB than 𝐼 0,𝑡𝑒𝑟𝑟 values obtained with P-528 and DO-235C propagation models. This observation can be explained by the fact that free space losses are higher than propagations losses computed from P-528 and DO-235B models, in particular when the horizontal distance between the terrestrial RFI source and the aircraft antenna is larger than about 600m-700 (see Figure 11).

Moreover, 𝐼 0,𝑡𝑒𝑟𝑟 in EDDF has the same magnitude than 𝐼 0,𝑡𝑒𝑟𝑟 in KLAX. This observation can be justified by the fact that the distribution of terrestrial emitters around the analyzed point is similar in EDDF and KLAX, as illustrated in Figure 13. In Figure 13, the red points represent the aircraft location at a boarding gate. The green and red points are respectively circles with radius 100 m and 1000 m centered on the aircraft position. Figure 13 shows that the distribution of terrestrial emitters is similar in EDDF and KLAX. In the immediate vicinity of the aircraft, the aircraft is mostly surrounded by a no-emission zone. Then, terrestrial emitters from the airport terminals impact the GNSS antenna. Terrestrial emitters outside the airport area are more distant from the aircraft.

Figure 13

Terrestrial emitters distribution around aircraft location (taxiing operations)

Cat II DH TABLE 7 provides the value of 𝐼 0,𝑡𝑒𝑟𝑟 for Cat II operations. The results are given for both the L1/E1 band and L5/E5a band, in Frankfurt EDDF and Los Angeles KLAX. In both cases, the DFMC Cat II antenna model is used. The aircraft location is determined so that the aircraft, following a 3° approach slope, land at the 07L and 25L runway touch point in EDDF and KLAX respectively. 6 -157.7 -181.1 -149.1 -157 -178.6 -148.8 -151.2 -156.5 -146.3 -149.1 -153.9 𝐼 0,𝑡𝑒𝑟𝑟 values are strongly dependent on the propagation channel model. DO-235C propagation channel model provides the lowest 𝐼 0,𝑡𝑒𝑟𝑟 value. That can be explained by the higher propagation losses in DO-235C as long as the horizontal distance exceeds 100 m, which is the case for the majority of terrestrial emitters. Similarly, free space model provides a higher 𝐼 0,𝑡𝑒𝑟𝑟 value than P-528, because long range propagation losses, which concerns most of terrestrial sources as can be seen in Figure 14 for the evaluated aircraft position, are lower with the free space losses propagation channel model. Indeed, it can be seen from Figure 14 that the aircraft is place in the zero-emission zone and thus the main terrestrial emitters contributors to 𝐼 0,𝑡𝑒𝑟𝑟 are situated at a distance equal to or higher than 800m-1km.

Finally, conversely to the taxiing analysis, 𝐼 0,𝑡𝑒𝑟𝑟 is here strongly dependent on the choice of the airport. Indeed, 𝐼 0,𝑡𝑒𝑟𝑟 is higher in KLAX than in EDDF. That can be explained by the distribution of terrestrial emitters and the configuration of the Cat II approach illustrated in Figure 14. Cat I DH TABLE 8 provides the value of 𝐼 0,𝑡𝑒𝑟𝑟 for Cat I operations. The results are given for both the L1/E1 band and L5/E5a band, in Frankfurt EDDF and Los Angeles KLAX. In both cases, the DFMC Cat I antenna model is used. The aircraft location is determined so that the aircraft, following a 3° approach slope, land at the 07L and 25L runway touch point in EDDF and KLAX respectively. -154.7 -174.9 -148.6 -153.1 -172.4 -148.3 -150 -154.3 -145.8 -147.7 -151.6 Comments similar to Cat II analysis can be done: -𝐼 0,𝑡𝑒𝑟𝑟 depends on the choice of the propagation model and free space losses propagation channel model provides the highest 𝐼 0,𝑡𝑒𝑟𝑟 , DO-235C model is the less conservative model. -𝐼 0,𝑡𝑒𝑟𝑟 is location dependent and depends on the distribution of terrestrial emitters which is illustrated in Figure 15. In particular, there are more terrestrial emitters in the immediate vicinity of the aircraft in KLAX runway 25L than in EDDF runway 07L, leading to a higher 𝐼 0,𝑡𝑒𝑟𝑟 value. As P-528 attenuation is stronger than free space losses, 𝐼 0,𝑡𝑒𝑟𝑟 obtained with P528 propagation channel model is higher than the value obtained considering free space loss propagation channel model. The difference of results between the L1/E1 and L5/E5a frequency bands are around 2.5 dB. This difference is due to better propagation characteristics of the L5/E5a band and is equal to the square ratio of carrier frequencies.

High altitude TABLE 10 provides the value of 𝐼 0,𝑡𝑒𝑟𝑟 for high altitude operations. The results are given for both the L1/E1 band and L5/E5a band, at the EDDF TIXAK FAF. In both cases, the DFMC Cat I antenna model is used. Similar to the intermediate altitude analysis, P-528 propagation channel model provides a lower 𝐼 0,𝑡𝑒𝑟𝑟 value as it is less conservative than free-space loss model. In addition, 𝐼 0,𝑡𝑒𝑟𝑟 is higher in the L5/E5a band than in the L1/E1 frequency band.

Discussion and conclusions

From the results presented above, main conclusions on the impact analysis of the terrestrial emitters on GNSS receivers can be extracted:

-The choice of the propagation channel model is significant. Indeed, free space loss model, which is used in DO-235B to estimate the impact of terrestrial emitters on GNSS receivers is pessimistic at low altitude in comparison to P-528 and DO-235C model, since it provides a higher 𝐼 0,𝑡𝑒𝑟𝑟 value. Conversely, DO-235C model, which has been elaborated from empirical data, is optimistic in comparison to the two other models as it considers strong propagation attenuation when the horizontal distance between the aircraft and the terrestrial emitter source exceeds several hundreds of meters. This overestimation of the propagation channel attenuation results in a maximum difference of almost 30 dB of 𝐼 0,𝑡𝑒𝑟𝑟 (case of Cat II approach in EDDF in the L1/E1 band). At intermediate and high altitude, DO-235C model is no longer applicable. In this case, the difference between P-528 and free space models is reduced to less than 1 dB.

-Density of terrestrial emitters function is also important to be precisely modeled in order to have a precise estimation of 𝐼 0,𝑡𝑒𝑟𝑟 . In particular at low altitude, 𝐼 0,𝑡𝑒𝑟𝑟 is higher if the aircraft is surrounded by terrestrial sources in the immediate vicinity of the aircraft (instead of being in a zero-emission zone inside the airport environment). Moreover, at intermediate and high altitude, the environment below the aircraft has an influence on 𝐼 0,𝑡𝑒𝑟𝑟 , irrespective of the considered propagation model (but note that the estimated 𝐼 0,𝑡𝑒𝑟𝑟 still depends on the selected propagation channel model). 𝐼 0,𝑡𝑒𝑟𝑟 is indeed more significant when the aircraft flies above urban areas.

-Propagation losses increase with the aircraft altitude. However, comparison of TABLE 6, TABLE 7 and TABLE 8 shows that 𝐼 0,𝑡𝑒𝑟𝑟 increases with altitude up to Cat I DH. Indeed, 𝐼 0,𝑡𝑒𝑟𝑟 does not only depend on the propagation losses which decrease with altitude (distance between the RFI source and the aircraft), but also on the environment below the aircraft. In particular, the number of terrestrial sources visible by the aircraft is indeed higher at Cat I DH than Cat II DH or taxiing, which is the dominant factor rather than the propagation losses. For intermediate and high-altitude cases, 𝐼 0,𝑡𝑒𝑟𝑟 decreases with the aircraft altitude meaning that in these situations the dominant factor has changed and it is now the propagation losses. This observation also highlights the necessity to precisely describe the terrestrial emitters environment when estimating the impact of terrestrial sources on GNSS receiver.

-Eventually, the level of terrestrial emitters equivalent noise is higher in the L5/E5a frequency band than in L1/E1 band. Propagation losses are indeed lower in the L5/E5a band due to the lower central frequency (𝑓 𝐿5 < 𝑓 𝐿1 ), leading to a higher received power from terrestrial emitters at the aircraft GNSS antenna in the L5/E5a band.

CONCLUSION

This article analyzes the impact on civil aviation GNSS receiver of various electronic terrestrial emitters which act as RFI sources transmitting spurious emission. The impact of terrestrial emitters RFI on GNSS receiver capability to meet minimum signal processing requirements is equivalently modeled as an increase of the GNSS RF module noise floor. The equivalent noise 𝐼 0,𝑡𝑒𝑟𝑟 induced by terrestrial emitters characterizes this increase of the noise floor. This article proposes a 𝐼 0,𝑡𝑒𝑟𝑟 calculation from a statistical approach, considering that:

-The number of emitters in view from a RF point of view is random and follows a Poisson distribution.

-The location of each terrestrial emitter is random, and its distribution is derived from density of population.

This article also describes and discusses main 𝐼 0,𝑡𝑒𝑟𝑟 calculation inputs. The propagation channel model is included among these inputs. Three relevant propagation channel models are discussed and compared: free space losses (FSL) model, the model used in DO-235C, issued from empirical propagation attenuation analysis, and ITU P-528 propagation model which is applicable for signal transmission in the VHF and UHF frequency bands. Among the three propagation channel models, FSL presents the lowest attenuation (propagation losses) and therefore leads to the highest 𝐼 0,𝑡𝑒𝑟𝑟 . DO-235C model presents stronger attenuation values when the horizontal distance between the aircraft and the terrestrial emitter exceeds some hundreds of meters. Therefore, DO-235C provides lower 𝐼 0,𝑡𝑒𝑟𝑟 values than FSL values. One drawback of DO-235C is its limited applicability as it can only be applied to low altitude operations. ITU P-528 propagation model differs from FSL model on the long-range segment.

𝐼 0,𝑡𝑒𝑟𝑟 is also linked to the density of terrestrial emitters. Two distributions of terrestrial emitters are proposed in this article. For low altitude operations around airports, as the RFI environment has a key influence on 𝐼 0,𝑡𝑒𝑟𝑟 , airport zone is supposed to have no emitter outside airport buildings. For higher altitude operations, terrestrial emitters distribution is derived from population density. From the analysis performed in this article, it appears that 𝐼 0,𝑡𝑒𝑟𝑟 increases with the altitude up to Cat I DH as the number of terrestrial emitters visible from the aircraft also grows. Above Cat I DH, 𝐼 0,𝑡𝑒𝑟𝑟 decreases along the altitude as propagation losses increase (due to the higher distance between terrestrial emitters and the victim aircraft) although the number of terrestrial emitters in view from a RF perspective also increases.

For standardization consideration, it is recommended to use the P-528 model when computing 𝐼 0,𝑡𝑒𝑟𝑟 because this model is validated at an ITU level and represents a good compromise between the model of DO-235C and free space loss. The terrestrial emitters density derived from population density (for high and intermediate altitude analysis), or from the modeling of the airport environment (at low altitude), appear to be relevant for standardization purpose. Concerning the L5/E5a frequency band, 𝐼 0,𝑡𝑒𝑟𝑟 reaches a maximum of -147.7 dBW/MHz with P-528 (Cat I DH, LAX). 𝐼 0,𝑡𝑒𝑟𝑟 is then reduced at intermediate and high altitude, and is equal to -161.6 dBW/MHz at high altitude with P-528 model. It is recommended to consider these values as part of standardization of L5/E5a receivers (DO-292).

Some uncertainties remain in the calculation of 𝐼 0,𝑡𝑒𝑟𝑟 . Even though these uncertainties are covered in the GNSS performance assessment by the consideration of a 6 dB safety margin, it might be interesting to improve the precision of 𝐼 0,𝑡𝑒𝑟𝑟 estimation. In particular, additional test on spurious emission of recent electronic devices would be appreciated to validate the 𝑃 0 value which is used in this article and that is issued from testing performed in 2003 by NASA.

  𝐼 0,𝑡𝑒𝑟𝑟 = 𝔼(𝑃 𝑡𝑒𝑟𝑟 ) = 𝔼 (𝑃 0 ∑ 𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 )𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 ) 𝑁 𝑘=1 ) = 𝑃 0 𝔼 𝑁 (𝔼(∑ 𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 )𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 ) 𝑁𝑃 0 𝔼 (𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 )𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 ))) = 𝑃 0 𝔼(𝑁)𝔼 (𝐺 𝑅𝑋 (𝑟 𝑘 , 𝜑 𝑘 )𝐿 𝑝𝑟𝑜𝑝 (𝑟 𝑘 , 𝜑 𝑘 )) = 𝑃 0 𝑁 ∫ ∫ 𝐺 𝑅𝑋 (𝑟, 𝜑)𝐿 𝑝𝑟𝑜𝑝 (𝑟, 𝜑)𝑓(𝑟,

Figure 1

 1 Figure 1 Maximum Lower Hemisphere DO-301 and DO-373 Antennas Gain

Figure 2

 2 Figure 2 represents the position of the ground RFI source with respect to the aircraft. The position of the emitter with respect to the victim receiver is characterized by the horizontal distance 𝑟 (slant range projection on the ground) and the elevation angle 𝜑. Let us define 𝐻 𝑎 as the height of the aircraft GNSS antenna above ground level (AGL), and ℎ 𝐵 the AGL height of the terrestrial source.

Figure 2

 2 Figure 2 Position of the terrestrial emitter with respect to the aircraft GNSS antenna Under flat Earth assumption, the elevation angle of the terrestrial emitter signal inputting the GNSS antenna, which is needed to compute de GNSS receiver antenna gain from Figure 1, can be expressed as in Equation (7).

  Fading:DO-235C model also considers fast and slow fading. Fast and slow fading are modeled by the introduction of a multiplication term. DO-235C propagation model including fading is expressed as in Equation (20) (see(RTCA, 2011)).𝐿 𝑝𝑟𝑜𝑝𝑎 = 𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 𝛼 exp ( (𝑙𝑛(𝜎)) 2 2 )(20) 𝛼 is the mean value of a non-central Chi-square fast fading distribution and 𝜎 belongs to the slow fading distribution. Consistently with the values used in DO-235C (RTCA, 2022), 𝛼 and 𝜎 are numerically represented on Figure 4.

Figure 4

 4 Figure 4 Fast and slow fading coefficients Note that the parameters 𝛼 and 𝜎 used here are consistent with (FAA, 2012). 𝑅 1 and 𝑅 2 determination: Breakpoints 𝑅 1 and 𝑅 2 between the different zones must be determined. They depend on the domain of validity of each propagation model. Moreover, 𝑅 1 and 𝑅 2 values depend on the medium range path loss model: -If the Erceg Greenstein model if chosen (𝐻 𝑎 ≤ 80 𝑚), then 𝑅 1 and 𝑅 2 are determined to ensure the continuity of the complete propagation losses, 𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 , between zones 1 and 2. In other words, 𝑅 1 and 𝑅 2 are respectively the intersection points between two-rays propagation loss model and Erceg Greenstein model, and between Erceg Greenstein model and Hata Okumura model. -If 𝐻 𝑎 > 80 𝑚, Erceg Greenstein model cannot be applied anymore, and a linear logarithm slope is used between 𝑅 1 and 𝑅 2 .In this case, 𝑅 1 and 𝑅 2 are computed as follows. 𝑅 1 is determined to target a grazing angle verifying tan(𝜃(𝑅 1 )) ≈ 0.5. Some adjustments can be done to avoid strong propagation loss slope breakout. 𝑅 2 is determined such that the elevation angle toward the antenna is around 4°: 𝐸(𝑅 2 ) ≈ -4°, where 𝐸(𝑅 2 ) is computed using Equation (7).

Figure 5

 5 Figure 5 Zones with Different Density of Emitters in DO-235B En-Route OperationsZone A represents the urban areas above which the aircraft if flying.(RTCA, 2008) and[START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] consider that a reasonable estimate for the density of emitters in a urban area is 10 -4 𝑚 -2 . Indeed,[START_REF] Peterson | Analytic Statistical Model for Aggregate Radio Frequency Interference to Airborne GPS Receivers from Ground-Based Emitters: Analytic Statistical Model for Aggregate Radio Frequency Interference[END_REF] mentions that cellphones activity in Los Angeles does not exceed 0.78 active emitter per hectare. Therefore, 10 -4 active terrestrial emitter per square meter seems to be a reasonable upper bound. Zone A radius is set to 30 km in DO-235B. Zones B and C respectively refers to suburban and rural areas. Density of emitters in suburban areas is set to

Figure 7

 7 Figure 7 Density of population across Western Europe

  , 𝜑) = 𝛼 𝑡𝑒𝑟𝑟 𝑑 𝑝𝑜𝑝 (𝑟, 𝜑) (21) 𝛼 𝑡𝑒𝑟𝑟 = 5 10 -2 is the linear regression slope issued from Figure 8 data. Unit of 𝑑(𝑟, 𝜑) is emitter per square meter if 𝑑 𝑝𝑜𝑝 (𝑟, 𝜑) is in inhabitant per square meter. Eventually, Figure 9 represents the density of terrestrial emitter in Western Europe, computed following Equation (21).

Figure 9

 9 Figure 9 Density of terrestrial emitters in western Europe

Figure 10

 10 Figure 10 Distribution of terrestrial emitters in an airport environment, at EDDF (left hand) and KLAX (right hand)

Figure 14

 14 Figure 14 Terrestrial emitter distribution around aircraft location (Cat II operation) According to Figure 14, the distribution of terrestrial emitters is different between EDDF and KLAX. The immediate vicinity of the aircraft is a zero-emission zone in EDDF runway 07L, whereas terrestrial emitters outside the airport area are relatively close to the aircraft position in KLAX runway 25L. Terrestrial emitters close to the aircraft has a stronger impact because of lower propagation losses. As a consequence, 𝐼 0,𝑡𝑒𝑟𝑟 is higher in KLAX runway 25L than in EDDF runway 07L.

Figure 15

 15 Figure 15 Terrestrial emitter distribution around aircraft location (Cat I operation) Intermediate altitude TABLE 9 provides the value of 𝐼 0,𝑡𝑒𝑟𝑟 for intermediate altitude operations. The results are given for both the L1/E1 band and L5/E5a band, at the EDDF TIXAK FAF. In both cases, the DFMC Cat I antenna model is used.

TABLE 1Coefficients

 1Coefficients

	of Erceg-Greenstein propagation model in airport area
	𝑎	𝑏	𝑐
	3.6	0.005	20
	Zone 3:

  The propagation losses for the proposed propagation model in a suburban area such as an airport environment can be expressed as in Equation (14) (for the originally model where 𝑟 < 20 km).𝐴 and 𝐵 are empirically determined coefficients. 𝐾 is a correcting factor to translate the model from urban to suburban areas. Coefficients 𝐴, 𝐵 and 𝐾 are given in Equation (15).

	𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 = 10 -	𝐴+𝐵 log( 10 1000 𝑟	)-𝐾	(14)
	𝐴 = 69.55 + 26.16 log(𝑓 𝑐 ) -𝑎(ℎ 𝐵 ) -13.82 log(𝐻 𝑎 )
	𝐵 = 44.9 -6.55 log(𝐻 𝑎 ) 𝐾 = 5.4 + 2 (log ( 28 )) 𝑓 𝑐 2

is applied. The Hata-Okumura propagation model aims to model urban environments and was established through empirical measurements. It can be applied to GNSS receiver antennas whose altitude is up to 200 m and for horizontal distance 𝑟 below 20 km; for 𝑟 > 20 km a modification recommended by ITU is additionally applied.

  When 𝑟 ≥ 20 𝑘𝑚, the Hata Okumura model is no longer applicable. In this case, propagation losses expressed in Equation (17) are slightly modified according to ITU P-529 recommendations(ITU, 1999) and are given by Equation (18).

	𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 = 10 -69.12+26.16 log(𝑓 𝑐 )-13.82 log(ℎ 𝐵 )-3.2(log(11.75ℎ 𝐵 )) 2 -2(log( 10	𝑓 𝑐 28 ))	2	+(44.9-6.55 log(ℎ 𝐵 )) log(	𝑟 1000	)	(17)
	𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 = 10 -69.12+26.16 log(𝑓 𝑐 )-13.82 log(ℎ 𝐵 )-3.2(log(11.75ℎ 𝐵 )) 2 -2(log( 10	𝑓 𝑐 28 ))	2	+(44.9-6.55 log(ℎ 𝐵 ))(log(	𝑟 1000	))	𝐹(𝐻 𝑎 )	(18)
	𝐹(𝐻 𝑎 ) = 1 + (0.014 + 1.87 10 -4 𝑓 𝑐 + 1.87 10 -3	𝐻 𝑎 1 + 7 10 -6 𝐻 𝑎	) (log (	𝑟 20 10 3 ))	0.8	(19)

TABLE 3

 3

	Density of emitters function characteristics of zones B, C and D in DO-235B low altitude scenarios
		Zone A	Zone B	Zone C	Zone D
	Environment	Rural and suburban	Urban	Hot spot	No-emission zone
	Density of terrestrial emitters (𝑚 -2 )	1 3	10 -4	10 -4	10 -3	0

TABLE 4

 4 summarizes the assumptions of 𝐼 0,𝑡𝑒𝑟𝑟 calculations for the different flight phases. This section compares the propagation models at the different altitudes identified in TABLE4. DO-235C propagation model breakpoints (zone borders) for the taxiing, Cat II DH and Cat I DH altitudes are summarized in TABLE5. ℎ 𝐵 = 1.8 𝑚 is a representative height for terrestrial emitter.

	TABLE 4
	Assumptions on 𝐼 0,𝑡𝑒𝑟𝑟 calculation for the different flight phases
			Taxiing		Cat II decision	Cat I DH	Intermediate	En-route
					height (DH)		altitude
		𝐻 𝑎	5.64 m		25.94 m		53.34 m	520 m	FL400
	Location	Frankfurt EDDF	Frankfurt EDDF	Frankfurt EDDF	Frankfurt EDDF	50.5°N, 5.5°E
			Los	Angeles	Los	Angeles	Los	Angeles	FAF (TIXAK)
			KLAX		KLAX		KLAX
	Propagation model	Free space		Free space	Free space	Free space	Free space
			P528		P528		P528	P528	P-528
			DO-235C		DO-235C	DO-235C
	Density of Terrestrial	Airport		Airport		Airport	Population	Population
	emitters function	environment	environment	environment	density	density
	Antenna model	Cat II/III		Cat II/III		Cat I	Cat I	Cat I
	4.2	Comparison of propagation models

TABLE 5

 5 Figure 11 represents propagation losses for DO-235C, P-528 and free space propagation models 𝐿 𝑝𝑟𝑜𝑝𝑎 𝑤𝑖𝑡ℎ𝑜𝑢𝑡 𝑓𝑎𝑑𝑖𝑛𝑔 for taxiing, Cat II and Cat I operations.

	DO-235C propagation model characteristics for Taxiing, Cat I and Cat II operations
		Taxiing		Cat II DH		Cat I DH
	𝐻 𝑎 (𝑚)	5.64		25.94		53.34
	𝑅 𝐿𝑂𝑆 (𝑘𝑚)	9.78		20.98		30.09
	Frequency band	L1/E1	L5/E5a	L1/E1	L5/E5a	L1/E1	L5/E5a
	𝑅 1 (𝑚)	118.51	98.8	114.7	110.9	101.21	95
	𝑅 2 (𝑚)	286.59	288.06	2480	2543.57	11263	11782.1

TABLE 6

 6

	𝐼 0,𝑡𝑒𝑟𝑟 values for taxiing operations
				Frankfurt EDDF					Los Angeles KLAX
	Aircraft location		50.047120°N 8.572768°E			33.939747°N 118.406580°W
			L1/E1			L5/E5a			L1/E1			L5/E5a
		Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C
	𝐼 0,𝑡𝑒𝑟𝑟 (dBW/MHz)	-151

TABLE 7

 7

	𝐼 0,𝑡𝑒𝑟𝑟 values for Cat II operations
				Frankfurt EDDF					Los Angeles KLAX
	Aircraft location		50.032031°N 8.531156°E			33.937649°N 118.379546°W
			L1/E1			L5/E5a			L1/E1			L5/E5a
		Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C
	𝐼 0,𝑡𝑒𝑟𝑟 (dBW/MHz)	-151.

TABLE 8

 8

	𝐼 0,𝑡𝑒𝑟𝑟 values for Cat I operations
				Frankfurt EDDF					Los Angeles KLAX
	Aircraft location		50.031272°N 8.527074°E			33.938272°N 118.372765°W
			L1/E1			L5/E5a			L1/E1			L5/E5a
		Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C	Free space	P-528	DO-235C
	𝐼 0,𝑡𝑒𝑟𝑟 (dBW/MHz)	-151.1

TABLE 9

 9

	𝐼 0,𝑡𝑒𝑟𝑟 values for intermediate altitude operations
			TIXAK FAF
		L1/E1		L5/E5a
		Free space	P-528	Free space	P-528
	𝐼 0,𝑡𝑒𝑟𝑟 (dBW/MHz)	-162.3	-163.1	-159.8	-160.8

TABLE 10

 10

	𝐼 0,𝑡𝑒𝑟𝑟 values for high altitude operations
		High altitude European L5/E5a hot spot (50.5°N, 5.5°E)
		L1/E1		L5/E5a
		Free space	P-528	Free space	P-528
	𝐼 0,𝑡𝑒𝑟𝑟 (dBW/MHz)	-163.3	-164.1	-160.8	-161.6

Philippe ESTIVAL is a GNSS expert at the French civil aviation service provider (DGAC/DSNA/DTI) currently involved in EGNOS mission requirements and standardization activities on future multi-constellation GNSS receivers Eurocae WG62 and within EU/US cooperation agreement Working Group C. He has been in charge of CCF (EGNOS Central Control Facility) operational support and evolutions for PACF department until 2009 and then for ESSP SAS System Operation Unit up to 2015. then was involved at DSNA/DTI, the French Air Navigation


Simon Le Berre
email: simon.leberre@cea.fr

Isabelle Ramière
email: isabelle.ramiere@cea.fr

David Ryckelynck
email: david.ryckelynck@minesparis.psl.eu

Error analysis of the Hybrid Hyper-Reduction method for frictionless contact problems

Keywords: Hyper-Reduction, reduced order models, frictionless contact problems, Lagrange multipliers method, Error estimators

The solution of contact problems is often time consuming because of the non-linearity and nonregularity of the contact constraints. For contact problems treated by Lagrange multipliers (ensuring the Signorini contact conditions to be respected), applying reduced-order modeling remains challenging due in particular to the non-negativity constraint on the Lagrange multipliers [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Kollepara | On the limitations of low-rank approximations in contact mechanics problems[END_REF].

We focus here on the Hybrid Hyper-Reduction (HHR) approach for frictionless elastostatic contact problems proposed in [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Le Berre | Condition Number and Clustering-Based Efficiency Improvement of Reduced-Order Solvers for Contact Problems Using Lagrange Multipliers[END_REF]. This approach is based on a reduced integration domain (RID) built through a discrete empirical interpolation method as described in [START_REF] Ryckelynck | Hyper reduction of mechanical models involving internal variables[END_REF]. A reliable reduced dual basis is obtained by restricting the full order dual basis to the RID. The HHR model preserves the saddle-point formulation of the Lagrange multipliers method, and allows the verification of the Signorini conditions inside the RID (without projection). Numerical experiments proved the capacity of the HHR method to provide accurate results both for the displacement and the Lagrange multipliers solutions [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Le Berre | Condition Number and Clustering-Based Efficiency Improvement of Reduced-Order Solvers for Contact Problems Using Lagrange Multipliers[END_REF], even for large parametric spaces [START_REF] Le Berre | Condition Number and Clustering-Based Efficiency Improvement of Reduced-Order Solvers for Contact Problems Using Lagrange Multipliers[END_REF]. In this work, we introduce a priori and a posteriori estimators of the HHR method, providing error bounds to quantify the discrepancy between the HHR and FE (Finite Element) solutions. Thanks to the hybrid formulation of the HHR model, a Cea-like primal bound is obtained on the RID (see Figures 1a and2a). Indeed, the FE-like dual reduced basis hybrid enables the corresponding bounding term to vanish. Classical Brezzi-like analysis leads to bound the dual approximation error by the primal one (see Figure 1b). Hence, the final a priori error estimators for both primal and dual errors, depend on the distance between the primal FE solution and the HHR primal approximation space only. Furthemore, an a priori primal error bound can be derived on the whole computational domain, following Cea's lemma for Petrov-Galerkin formulation (see Figure 2b). Primal and dual a posteriori error estimators (depending only on HHR solutions) are also introduced to propose certified HHR error bounds, useful in particular for a efficient greedy snapshot selection. These estimators are based on the HHR equilibrium and constraints equations residuals and will be introduced in the light of those already present in the literature for ROM methods on mixed problems [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF][START_REF] Haasdonk | A reduced basis method for parametrized variational inequalities[END_REF][START_REF] Gerner | Certified Reduced Basis Methods for Parametrized Saddle Point Problems[END_REF].

  Dual error bounded by the primal error

Figure 1 :

 1 Figure 1: Error bounds on the RID for the Hertz's half disks test case

  (a) On the RID. (b) On the whole domain.

Figure 2 :

 2 Figure 2: A priori primal error bounds for an industrial test case with 2 parameters.