Application of Path-Based Multipath Component Tracking on Air-Ground Channel Measurement Data

The operation of unmanned aircraft is unthinkable without reliable wireless communication links: Despite a comparatively high expected level of autonomy of unmanned aircraft, monitoring and remote controlling are required to some degree. As the design of every communication link requires good knowledge on the characteristics of the communication channel, we have performed a measurement campaign to collect channel sounding data for the wireless air-ground channel at C-band. While we have focused on the campaign description and the analysis of the dominant signal component in previous publications, we now concentrate on the detection and tracking of multipath components. In this paper, we present our data processing chain that allows a fast parallel processing of the measurement data, as data dependencies are reduced as much as possible. We furthermore introduce a path-based multipath component tracking approach and apply it to our measurement data. This tracking allows us to estimate the location of reflectors causing multipath component signals. We apply our processing chain to data recorded during take-off at a small airport and compare the results of the reflector localization to a satellite image of the airport to successfully verify our approach.


Application of Path-Based Multipath Component
Tracking on Air-Ground Channel Measurement Data potential communication partners. The type and amount of data that needs to be transmitted has been widely discussed in [3], where both satellite-based and terrestrial (i.e. air-ground) data links are investigated.
For the development of reliable, high-throughput wireless data links, a good understanding of the physical characteristics of the wireless channel between the communication partners, in our case of a terrestrial system the Ground Station (GS) and the UA, is necessary. Wireless channels are characterized by reflection, diffraction, and scattering; these effects are summarized as multi-path propagation [4]. Channel measurements are a common procedure to gain knowledge on wireless channels. This knowledge can later be used for the development and validation of a channel model. The resulting channel model may cover several scenarios -in case of the air-ground channel, those scenarios can be the different phases of a flight like take-off, en-route, and landing. A further differentiation, e.g. based on the ground surface, is also considered in literature, e.g. [5], [6]. The findings from the channel measurements and the resulting channel model can be used to design and evaluate waveforms for a new wireless communication system.
Inherently, channel measurements for aeronautical channels are quite complex and costly as they require at least one air vehicle. Both costs and effort may rise significantly when jet aircraft operating at high speeds at several kilometers of altitude are involved. Additionally, there was no high demand for sophisticated models for aeronautical channels in the past since modern digital communication links did not play a role in civil aviation for decades. These two aspects might be the reasons, why aeronautical channels in general and the aeronautical air-ground channel in particular have not received as much attention as other channels like those used in cellular networks, e.g. LTE and 5 G. However, we want to provide a brief overview of the available literature in the field of aeronautical (air-ground) channel models. A comprehensive review on air-ground channels, with a special focus on UAs can be found in [7].
In [8], [9], the application of simple tap-delay line models on aeronautical channel modeling has been discussed. The authors have already stated, that the different phases of a flight require different models or at least an individual parametrization.
The development of digital data links for civil aviation motivated the development of more sophisticated channel models. surface data link called Aeronautical Mobile Airport Communication System (AeroMACS) [10]; however, according to the application of Aeronautical Mobile Airport Communication System (AeroMACS), the model focuses on the airport surface only. An alternative approach for modeling the airport surface channel, also in C-band, has been introduced in [11].
In contrast to AeroMACS, the L-Band Digital Aeronautical Communication System (LDACS) is a communication link designed for all flight phases of manned aircraft, predominantly developed by the German Aerospace Center (DLR). As part of its development, a channel measurement campaign in L-band was performed and a new approach of channel modeling was proposed [12], [13]. The channel model is based on the evaluation of the tracked Multipath Components (MPCs) and does not only incorporate statistical element but also integrates geometric aspects.
In [14], air-ground channel measurements with a bandwidth of 2 MHz in the Ultra High Frequency (UHF) band have been performed. The authors also provide a model on the path loss and discuss the observed effects of multipath propagation by a statistical evaluation of sets of Channel Impulse Responses (CIRs).
Comprehensive channel measurements in both C-and L-band have been performed by a team of the University of South Carolina and NASA. The evaluation has been published in a series of articles: While [5] described the campaign setup and proposed models for the air-ground channel in over-water scenarios, the team focuses on the modeling of the air-ground channel in hilly and mountainous terrain in [6]. In [15], the focus lies on the modeling of the air-ground channel in suburban and near-urban environments. The described models concentrate on the evaluation of sets of CIRs and not on the tracking of individual MPCs over time.
The German Aerospace Center (DLR) also performed a C-band channel measurement campaign with a jet aircraft that covered multiple flight scenarios. This measurement campaign was described in detail in [16], where we have also provided a close analysis of the dominant signal component of the received channel sounding signal. In our present paper we focus on the signal components besides the dominant component, the MPCs, since they, as stated above, strongly define the physical properties of a wireless channel. We address the task of identifying and tracking the MPCs in our measurement data by the application of a multi-stage processing chain. In contrast to [17] we do not apply a filter-based algorithm but a concept, where the MPCs are first identified and are later tracked by a novel path-based algorithm in a second step. This two-staged approach allows a parallelization of the computational expensive processing of the raw measurement data.
The paper is structured as follows: In Section II we briefly describe the setup of the measurement campaign the processed data has been collected. The pre-processing of this data is explained in Section III. Sections IV and V provide information on how the individual MPCs are extracted from the measurement data. In Section VI, we introduce our approach used for the tracking of the detected MPCs over time before we apply it to our measurement data in Section VII. We conclude our paper in Section VIII. Fig. 1. Block diagram of the processing described in this paper. The section where the corresponding processing step is discussed in detail is given next to the box. Fig. 1 also provides a graphical overview on how the part of the paper explaining the processing is structured.

II. MEASUREMENT CAMPAIGN
The DLR performed a measurement campaign in 2018 where the air-ground/ground-air channel in C-band was measured. The applied channel sounding signal had a bandwidth of about 50 MHz and used a carrier frequency of f c = 5.2 GHz. A detailed description of the campaign can be found in [16]; this section just intends to give a brief overview on the campaign setup.

A. Hardware Setup
The campaign setup consisted of a Ground Station (GS) and an Airborne Station (AS) aboard DLR's Falcon 20E aircraft. Through all experiments, the GS was the transmitter of the channel sounding signal and the Airborne Station (AS) was the receiver of the signal.
The most relevant devices of the GS are the Arbitrary Waveform Generator (AWG) and the High Power Amplifier (HPA) that amplifies the output signal of the AWG. The GS also contains a Global Navigation Satellite System (GNSS) receiver and an atomic clock. The latter two devices act as a GNSSdisciplined oscillator that is used as time base for the AWG.
The most important part of the AS setup is the IQ-recorder that samples the received signal, and counts and stores these IQ-samples. Additionally, the setup contains an Inertial Measurement Unit (IMU) and also a GNSS-disciplined oscillator that is used as a time base for the IQ-recorder. As the applied IQrecorder is not capable of handling the applied carrier frequency, a frequency mixer ("downconverter") is used to shift the signal to an intermediate frequency. This mixer also provides an Adaptive Gain Control (AGC), that allows an adaptive amplification or attenuation, respectively, of the received signal. It is therefore possible to adjust the receiver's operating point during recording e.g. depending on the Line of Sight (LOS) distance to the transmitter.

B. Measurement Signal and Procedure
The channel sounding signal consists of a so called channel sounding sequence of length 40.96 μs that is gaplessly repeated in an infinite loop on the GS's AWG. The channel sounding sequence is an iteratively filtered multitone signal whose initial phases are distributed as Newman Phases. A more detailed description of the signal generation can be found in Section II-B of [16].
As the receiver samples the incoming signal at f sr = 50 MHz, it is assured that N = 2048 consecutive samples of the received signal consist exactly one complete instance of the channel sounding sequence.
Before and after all measurement flights, the output of the GS's HPA is directly connected to the receiving hardware aboard the AS using attenuators and a cable to perform a reference measurement of the channel sounding sequence. The transmission and the recording of the channel sounding signal is then started for a few seconds. During post-processing, N consecutive samples of this recorded signal are cut out and called reference signal x ref ∈ C N . We also store the IQ-recorder's sample counter value ρ ref for the first sample of x ref . Additionally, the GNSS-disciplined oscillators of the GS and the AS are synchronized during this process.
During the measurement flight, the GS is in its transmission location on the rooftop of DLR's IKN building and the output of the HPA is connected to the transmitting antenna. The AS's receiving hardware is connected to the receiving antenna mounted at the bottom fuselage of the Falcon aircraft. During the flight, the AS's IQ-recorder stores the received samples in s Rx .
It is ensured, that the GS's AWG and the AS's IQ-recorder are continuously running during and between the pre-flight reference measurement, the actual measurement flight, and the post-flight reference measurement. Thus, we can simply use the sample counter ρ of the IQ-recorder as a reliable time base during offline processing, when the measurement data recorded during flight are evaluated using the reference signal x ref , see Sections III to V.

C. Flight Tracks and Maneuvers
All flights started and ended at the EDMO airport close to Munich, Germany. The campaign involved four flights where different flight scenarios were covered. The flight scenarios not only included typical en-route flight patterns at cruising speed at multiple cruising altitudes, but also take-offs, landings, go-arounds, and flights with extreme banking. The data collected during the different flight scenarios allows us to develop individual channel models for each scenario.

A. Ground Truth
Based on the data recorded by both GNSS-receivers and the IMU, a ground truth is computed. The ground truth provides information on the aircraft's position, orientation, and heading at a certain time instant. This allows us to estimate the LOS component's Free Space Path Loss (FSPL), the LOS component's delay, denoted by τ LOS , and the LOS component's Doppler shift, denoted by ν LOS .

B. Block Processing
The processing of the measurement data is performed block wise. A block b of BN consecutive samples is cut out of the received signal: where ρ (b) denotes the IQ-recorder's sample counter corresponding to block b. We call B ∈ N the block size, as it defines how many consecutive channel sounding sequence instances of length N are processed together in one block. A larger B increases the Doppler resolution during processing, however, this increase comes at the cost of time resolution: If B is chosen large, it is possible that relevant fast fading channel effects are vanished out and thus are not modeled appropriately later. This effect becomes even more likely considering the comparatively high velocity of an aircraft. Based on the ground truth data, x (b) is shifted such that the current LOS delay and LOS Doppler shift are compensated: where F ν τ {·} denotes a function shifting a signal by delay τ and frequency ν.
In a last step, we reshape the vector x (b) ∈ C BN to a matrix X (b) ∈ C N ×B , where each column corresponds to one instance of the channel sounding sequence, thus the items of X (b) are set according to: In the following, we assume that the power levels of x ref and x (b) are matched, such that the attenuation used during the recording of the reference signal and the FSPL of the LOS component during the recording of x (b) are compensated as suggested in [16, (11)].

IV. COARSE MULTIPATH COMPONENT DETECTION
In this section, we describe our approach for detecting Multipath Components (MPCs) in our measurement data. As shown in [17], an MPC l can be described by a quadruple where τ l denotes the delay shift, ν l the Doppler shift, and α l the complex weight of the MPC. 1 The weight α l is split into its absolute value and its phase to achieve ξ l ∈ R 4 , which not only allows a direct access to the MPCs' amplitudes, but also a more memory efficient processing. If not denoted otherwise, we assume that both τ l and ν l are given with respect to the LOS component.
The set of MPCs that have been detected in a processing block b is denoted by M (b) having cardinality L (b) . To avoid ambiguities, the block index b is added to the MPC notation according to ξ We now describe the processing steps for the coarse detection of MPCs or clusters of MPCs, respectively. The detected set of MPCs is therefore denoted by M (b,crs) . Please note that the described approaches have no data dependencies on previously processed blocks, which allows a parallel processing of an arbitrary amount of blocks at the same time.

A. Based on Impulse Response/Power Delay Profile
The basic idea of this approach is to detect peaks in the absolute value of the Impulse Response (IR) or the Power Delay Profile (PDP), respectively, of a data block. Each of these peaks represents either a single MPC or a cluster of MPCs; the position of the peaks can be used to estimate the underlying MPCs' delay τ l . However, no Doppler information ν l can be extracted (directly) following this approach.
1) Impulse Response/Power Delay Profile: The first step is to compute either the Impulse Responses (IRs) or the complex Power Delay Profile (PDP) of the current data block X (b) by applying the processing described in Appendix A or Appendix B, respectively. While the (coherent) PDP has the advantage of a lower noise floor, short term MPCs may vanish out if the block size B is chosen too large. On the other hand, if B is chosen too small, signals of some reflectors may vanish in the noise floor and will remain undetected.
As the further processing is the same for both cases, in the following, we refer to both the individual IRs of the current block and the PDP, respectively, by y (b) ∈ C Nf up . f up denotes the upsampling factor.
We also define the corresponding logarithmic vector 2) Noise Floor: The power of the noise floor for the current block σ (b) 2 is estimated based on areas of y (b) log where the appearance of measurable strong MPCs or measurement artifacts is very unlikely. Thus, the average power for the part of y (b) log where 1 Note that MPCs that cannot be resolved due to limited resolution are handled as a single MPC, although the more correct term would be MPC cluster.
{τ σ |τ max < τ σ < N/f sr } is computed, where τ max denotes the threshold delay up to which significant MPCs are expected.
3) Peak Detection: The next task is to detect local maxima ("peaks") in y (b) log . A peak must fulfill the following two conditions: r The peak's value must exceed a value of σ (b) 2 | dB + P thresh . r The peak must have a certain prominence defined by a minimum distance to other detected peaks. All detected L (b) peaks' positions n l inside of y (b) log are stored in P (b) , where n l denotes the index of the l-th peak.
4) Index Translation: Based on the resolution of the IRs given by (22), the indices found in the previous step can be translated to a delay according to The complex weight is given by In order to complete the four elements of the MPC as defined in (4), the Doppler shift of MPC l is set to ν l = 0. Then, ξ l is added to M (b,crs) .

B. Based on Doppler-Delay Spreading Function
The basic idea of this approach is very similar to the previous one, however, instead of using the IR/PDP, the Doppler-Delay Spreading Function (DDSF) is used. The Doppler-Delay Spreading Function (DDSF) provides not only information on the delay, but also on the Doppler shift of an MPC. Again, peaks in the DDSF represent either a single MPC or a cluster of MPCs.

1) Doppler-Delay Spreading Function:
The detection process starts with the computation of the DDSF according to Appendix C of the current data block b as defined above and take its logarithmic absolute value: where f up denotes the upsampling factor along the delay axis i.e. the columns of A (b) ; the rows correspond to the Doppler axis.
2) Noise Floor: The power of the noise floor for the current block σ (b) 2 is estimated based on areas of A (b) where the appearance of measurable strong MPCs or measurement artifacts is very unlikely. Thus, the average power for the areas 3) Peak Detection: The next task is to detect local maxima ("peaks") in A (b) . A peak must fulfill the following two conditions: r The peak's value must exceed a value of σ (b) 2 | dB + P thresh . r The peak must have a certain prominence defined by an elliptical footprint. This allows a more flexible search than a simple distance measure. All detected L (b) peaks' coordinates inside of A (b) are stored pairwise (n l , m l ) in P (b) , where n l denotes the row index and m l denotes the column index of the l-th peak.

4) Index Translation:
Based on the resolution of the DDSF given by (22), and (23), the index pairs found in the previous step can be translated to a delay and Doppler value. For the delay, this mapping is done according to (6), for the Doppler, the mapping is done according to

5) Least Squares Optimization:
To find the complex weights for the MPCs, we first create a matrix describing a synthesized version of the signal based on the delays and Doppler shifts of the MPCs detected above using the shift function F ν τ {.} introduced in Section III: ref denotes a vector with B concatenated instances of the reference signal. Thus, the l-th row of Y (b) contains a modified version of the reference signal, shifted in time and frequency according to τ l and ν l , respectively. Now, an approximation for the vector α (b) ∈ C L (b) needs to be found for the equation We apply the Least Squares (LS) algorithm for this task. Besides the desired approximation for α (b) , the Least Squares (LS) algorithm also returns a residual (b) that we understand as a measure of precision of the estimated MPCs' parameters. The residual is computed according to which corresponds to the sum over all elements of the vector that results from the squared Euclidean distance between the measured signal block and the weighted synthesized signal. By combining τ l , ν l , and the absolute value and phase of the l-th entry of α (b) , all elements of MPC ξ l , as defined in (4), are given. ξ l is now added to M (b,crs) .

A. Fine Detection
The potential problem of the approach described in Section IV-B is, that its precision is limited by the resolution of the discrete DDSF. Therefore, the exact position of an MPC inside of a tile of the DDSF matrix -and consequently its exact delay and Doppler shift -cannot be determined any further.
We address the this issue by applying an optimization algorithm that tries to adjust the delay and Doppler shift of every detected MPC within the bounds defined by the size of a tile of the DDSF. 1) Optimization Algorithm: As suggested in [17], we apply the Bound Optimization by Quadratic Approximation (BOBYQA) algorithm [18] to the given problem. The BOBYQA algorithm can be applied to multi-dimensional problems, does not require a derivative of the optimization objective, and allows the usage of optimization bounds.
2) Initialization: The optimization is initialized by the delay and the Doppler shifts of the MPCs given in M (b,crs) .
3) Dimensions: The optimization algorithm tries to achieve an optimal result by adjusting the delay and Doppler shift of every MPC. As the amount of assumed MPCs does not change during the processing of block b, the optimization problem is solved along 2|M (b,crs) | = 2L (b,crs) dimensions.

4) Optimization Bounds:
The optimization bounds for each MPC are given by the DDSF's tiles' dimensions. Thus, the lower bound (lb) and the upper bound (ub) for the delay optimization of MPC l are given by: and, correspondingly, for the Doppler optimization of MPC l, the bounds are given by

5) Optimization Process:
The algorithm basically executes the LS optimization described in Section IV-B5 in an iterative loop indexed by i ∈ N: i is generated using the delay and Doppler shifts of the MPCs given in the current (hypothetical) set of MPCsM i similar to (9). Then, an equation similar to (10) is solved forα i -again using the LS approach. Thus, the overall objective of the optimization problem is to minimize the residual (b) i returned by the LS algorithm applied to (14). Using the definition given in (11), this leads to

6) Termination:
The execution of the optimization is terminated when the result does not improve more than a certain threshold from one iteration to the next:

B. Cluster Resolution
Although the approach presented in Section V-A has the ability to determine the parameters of the MPCs even below the resolution of the DDSF, it cannot resolve multiple MPCs whose mutual distances in delay and/or Doppler are lower than the corresponding DDSF resolution. Thus, they can only be described as clusters and not get detected individually. In case this is not acceptable, we suggest the following approach to resolve the individual MPCs.
It is assumed that the coarse detection process as defined in Section IV-B has been completed; thus there is a set M (b,crs) . For each item in M (b,crs) , Ψ ∈ N copies are created. A random delay and Doppler shift is added to each of these copies according to where τ l and ν l correspond to the delay and Doppler shift of . So far, this approach alleges, that every item in M (b,crs) is in fact a MPC cluster. However, the idea is that due to the further processing, i.e. by applying the optimization algorithm from Section V-A, the weights of the "unnecessary" MPCs copies will tend to 0: The algorithm given in Section V-A is initialized using M (b,crs) Ψ . For each ξ l,ψ , the optimization bounds are the same as for the corresponding "parent" MPC ξ l .
The amount of optimization dimensions is 2ΨL (b,crs) . Once the termination criterion is fulfilled (see above), the resulting MPCs are given in M

VI. TRACKING
Tracking of MPCs describes the behavior of the components over time. This is often realized by filter-based approaches, e.g. in [17]. However, path-based approaches have also been discussed, e.g. in [19]. We now describe our path based approach to track MPCs over time, i.e. multiple consecutive blocks b.

A. Graph Representation
where D : R 4 × R 4 → R is a function expressing the distance between two MPCs by a real value. The entry of the l s -th row and the l e -th column of the distance matrix We furthermore assume, that all MPCs of block b are sorted in descending order based on their absolute weight: |α

B. Path Detection
We use a path, represented by a sequence γ = ((b, l) i ) i=0,1,...,β−1 , to describe the behavior of an MPC over time. Each element of a path is a tuple (b, l) that sufficiently describes the position of each node inside of G. To avoid ambiguities, we add a superscript to our notation (e.g. γ (a) ) when necessary. The set of all paths is denoted by P.
For the sake of simplification, we have split the description of our path detection algorithm into three parts: the initialization, an inner part, and an outer part.
1) Initialization: For each block b, a vector Γ (b) ∈ N L b is defined and all of its items are initially set to zero. The purpose of this vector is to indicate whether an MPC has been assigned to a path: If MPC l of block b gets assigned to a path, Γ (b) [l] = 1 is set.
A threshold d thresh ∈ R is defined that provides a maximum distance two MPCs are allowed to have within a path. If this threshold is exceeded, the respective MPC cannot be part of the given path.
For the sake of completeness, an ordered list (corresponding to P) is created in Listings 1, line 2 where all detected paths are stored.
The purpose of the outer part (see Listing 1) of the algorithm is to detect the root note of a new path. Once a root element has been found, the inner part, represented by the subroutine detect_path, is called to find the remaining elements of the path. In the following, we explain the outer part of the algorithm in detail: The loop defined in line 5 iterates over the blocks starting at the first block. Considering the graph in Fig. 2, this loop iterates over the blocks from left to right.
The loop in line 8 ensures that the algorithm keeps processing on the current block b, until all of its MPCs have been assigned to a path.
The loop defined in line 11, finally, iterates over the MPCs of the current block b. Unless the current MPC l has been assigned to a path yet (check in line 14), the actual path detection subroutine detect_path, corresponding to the inner part of our approach, is called in line 20. The arguments passed to this call are the current block b and MPC l as they are used as the root of the new path.
Once detect_path returns a new path γ, it is appended to the global path pool in line 26.
Please note that the update of Γ (b) is performed inside of the called subroutine. Due to this side-effect we do not call detect_path a function to highlight the lack of idempotence.
3) Inner Part: The inner part of the algorithm is given in the subroutine detect_path in Listing 2. As the arguments passed to this subroutine describe the root of a new path inside of the graph, this tuple is added as a first element (line 10) to the list representing the path γ.
During this subroutine, the MPC index of the last element of the path is represented by l; in the beginning, this corresponds to the MPC index of the root node (line 15).
Similar to the block loop in Listing 1, the block loop of the subroutine iterates from "left" ro "right" over the blocks, starting with the root block (line 18).
Within the loop, the distances from the current MPC l to all MPCs of the next block that are not yet assigned to a path are computed and the minimum is determined (line 30). Please note, that the actual call of D can be substituted by a lookup in the corresponding distance matrix D (b) .
In case either no unassigned MPCs can be found (check in line 24) or the determined minimum distance exceeds the threshold d thresh (check in line 33), the loop stops -corresponding to a termination of the current path.
If none of these conditions is fulfilled, the index of the MPC with the minimum distance is assigned to l (line 39) and the new node is appended to the path γ (line 42). Finally, the detected MPC is marked as assigned (line 47).
Once the loop has terminated, no matter what criterion caused the loop's termination, the subroutine returns the detected path γ.

C. Path Improvements
The algorithm presented in Section VI-B has properties that may result in a degraded MPC tracking under certain conditions. We suggest the following approaches to address these issues.

1) Path Merging:
In case an MPC remains undetected for the duration of a single (or more) block(s), e.g. because it does not exceed the detection threshold, an actual path is detected as multiple individual paths. To detect and connect these path segments, we evaluate the distances between the last element of all paths and the root elements of all other paths: We first set a threshold β max ∈ N that defines the maximum length of interruptions given in blocks that we consider as acceptable. We then iterate over all permutations of all paths γ (m) and γ (n) to detect all path combinations, where holds. The paths γ (m) and γ (n) are merged if (20) is fulfilled. In case this condition holds for multiple path pairs, the one with the lowest distance is chosen. So far, the stability of a path γ (m) was equal to its length β (m) . However, the process of path merging motivates the definition of a new measure of the path stability, as the path now may skip some blocks. We therefore define the path stability η (m) as the delta of the block index of the last element of the path end the block index of the root element of a path: 2) Short Path Elimination: MPCs that appear for the duration of just one block are likely to be the results of a mis-detection. In this case, the algorithm in Section VI-B assigns these single MPCs to individual paths of length one. To only keep track of MPCs that have been detected for a longer duration, a threshold β min ∈ N can be set, that defines a minimum length of paths. Consequently, all paths that do not fulfill are then eliminated. If the proposed MPC tracking approach is used for channel modeling, those MPCs that have been assigned to distinct paths which have not been eliminated can be modeled by a process that generates paths. However, the MPCs assigned to one of the eliminated paths should not be ignored: We interpret these MPCs as a result of diffuse scattering and suggest to model them by a separate process. Thus, β min is understood as a parameter that has an influence on the amount of MPCs that are modeled by the process that generates paths. The higher this amount, the smaller the amount of MPCs that are modeled by the separate process.

VII. RESULTS AND DISCUSSION
In this section we present some results from the different processing steps described above. Due to the power normalization assumed in Section III, all amplitudes are given in dB with respect to the transmitted power and compensated FSPL. Delays and Doppler shifts are given with respect to LOS if not denoted otherwise.
A. Detection 1) Take-Off: Fig. 3 shows a color-coded representation of the DDSFs of three consecutive blocks b ∈ {0, 1, 2} computed according to the steps described in Section IV-B. The underlying data were recorded during take-off and processed with a block size of B = 480, corresponding to a block duration of BN/f sr = 19.6 ms. The detected MPCs are highlighted by red markers.
All three figures show two main reflector areas: One that corresponds to the LOS path and its close environment and one that consists of multiple MPCs with a Doppler shift of around −1.1 kHz and a delay of 2.1 μs to 3.5 μs. A third, less intense reflector area can be seen for a Doppler shift of around 0.8 kHz and a delay of 4.0 μs to 4.3 μs. Only very few MPCs with a positive Doppler shift can be observed.
It can be observed, that most of the MPCs detected in block b = 0 are also detected for b = 1 and b = 2 and vice versa. However, this does not apply for some MPC: Especially those having a delay of more than 5.8 μs tend to be missed in some blocks during detection. At the same time, it can be observed that these MPCs not only have a comparatively huge delay, but are also comparatively weak in power. This behavior motivates the concept of Path Merging (see Section VI-C) during the MPC tracking process.
2) Flyover: Fig. 4 also shows a color-coded representation of the DDSFs of three consecutive blocks b ∈ {0, 1, 2}. The underlying data were recorded just seconds before a direct flyover of the AS over the GS in an altitude of around 3.2 km. All processing parameters were the same as for the take-off scenario described above.
First, it can be observed that the noise floor is about 7.3 dB stronger compared to the take-off scenario. We explain this by the applied AGC setting of the mixer, see Section II-A, during recording: As the LOS distance -and consequently the FSPLwas higher than during the recording of the data shown in Fig. 3, a higher pre-amplification of the received signal was required, including an increase of the noise floor.
It can also be observed, that the LOS component is blurred along the Doppler axis. We explain this by the fact, that the aircraft was moving, while the received data are processed in blocks. Theoretically, block processing as applied here requires a perfect snapshot of the received data where all parameters that affect the measurement remain exactly the same for the length of one block. For obvious reasons this is not possible in a real world scenario, especially when an aircraft is moving at cruising speed. However, to extract any kind of Doppler information from measurement data, some sort of block processing is required as the Doppler can be described as the time derivative of the delay.  Thus, a trade-off between available Doppler resolution and time precision has to be resolved.
The envelope of the detected MPCs shows a shape similar to a parabola. This is not a coincidence, but can be explained by the fact, that in the air-ground channel, most of the reflectors causing the MPCs are located on the ground comparatively close to the GS. At the same time, each reflector is located on the surface of a prolate spheroid having the transmitter and receiver in its focal points [20]. The size of the spheroid a given reflector is located on is determined by the corresponding MPC's delay. As most of the reflectors that cause the MPCs are distributed on the earth's surface, which can be modeled by a plane, the intersection of the prolate spheroid and the surface results in a parabola-like shape. Fig. 5 shows a 5 s cutout of the evolution of the MPC components represented as paths. The MPCs were detected using the approach shown in Section IV-B; the paths were found using the algorithm presented in Section VI using the delay-Doppler-only metric as defined in Appendix D. The underlying data were recorded while the aircraft was on the runway and accelerating for take-off. For the sake of clearness, short paths having a length less than 50 blocks have been excluded from the plots.

B. Path-Based MPC Tracking
1) Interpolated Data Points: Fig. 5(a) shows the paths as they were detected by the algorithm including the application of the path improvements presented in Section VI-C. Due to the path merging, paths with missing blocks are part of the detected path set. Here, a cubic spline interpolation was applied to estimate the values of these missing blocks.
Most of the Doppler shifts increase in value over time by becoming more negative.. This is caused by the increasing speed (and therefore higher absolute Doppler shift) of the aircraft during take-off and by the relative movement of the aircraft with respect to the reflectors like airport buildings, trees and fences.
The figure shows no paths with a positive Doppler shift. This matches our observations from Fig. 3, where no strong MPCs with a positive Doppler shift were detected. In both cases, the underlying data was recorded during the same take-off; however, the data shown in Fig. 3 was recorded shortly after the data presented in Fig. 5.
The Doppler shifts, especially the shifts of the paths with the lowest Doppler shifts, show a step-wise unsteady evolution. From the plot this step size can be determined as 25.4 Hz which corresponds to the Doppler resolution of the DDSFs the underlying MPCs got extracted from. We therefore consider these steps not as a natural effect, but a processing artifact that will be addressed below.
The amplitudes of the detected components, including the LOS component given in orange, show a typical small scale fading behavior. We assume this fading of the individual components is caused by non-resolvable MPCs that add up either constructively or destructively.
2) Filtered Data Points: Fig. 5(b) shows low-pass filtered versions of the paths presented in Fig. 5(a). The purpose of the filtering is to smooth the curves to compensate the effects introduced by the processing, e.g. mis-detection of MPCs in the DDSF or, most prominent, the unsteady evolution of the Doppler shifts due to the low resolution. Another cause could be erroneously assigned MPCs: If d thresh is chosen very large, some MPCs may be assigned to the same path although they are of different origin and therefore should be assigned to separate paths.
As neither the aircraft nor the reflectors are expected to do sudden movements that would explain these effects, it is a reasonable assumption that these effects are just processing artifacts in most of the cases.
The impact of the filtering can be clearly seen in the plot of the Doppler shifts: The step-wise leaps that dominated the plotted paths in the unfiltered case are now gone and the curves look much smoother.
The effect on the delay and the amplitudes, respectively, is not as strong as for the Doppler shifts; however, the paths now look smoother than before.

C. Reflector Localization
To validate the processing described above, the results of the MPC processing and tracking have been used to perform a reflector localization according to Appendix E. The estimated reflectors' locations are then compared to the actual location of potential reflectors given in a map. Fig. 6 shows a map of the EDMO airport. The position of the transmitter of the channel sounding signal is marked by a red cross. The track the aircraft has traveled while the data that got evaluated for the reflector localization has been recorded is denoted by a blue line; the aircraft started heading south-west.
The path based reflector localization works by superimposing the results of the reflector localization for each individual MPC of an MPC path. The longer a path is, i.e. the more persistent an MPC is, the better will be the result of the estimation, as outliers and even ambiguities are averaged out due to the superimposition. In the figure, the certainty about the position of a reflector is represented by the data point's opacity. A color code is used to let the reader distinguish, what data points belong to a common path. However, as the amount of detected paths exceeds the amount of available colors by far, the colors are not exclusively mapped to a specific path.
While the comparatively high delay resolution allows a pretty precise estimation of the ellipse the reflector is located on, the estimation of the angle suffers from the comparatively low Doppler resolution. Although the superimposition helps to improve the precision of the angle estimation, the shape of the data point clusters shows the lack of precision of the angle estimation. However, in most of the cases the algorithm was able to resolve the ambiguities caused by the Doppler-to-angle conversion: The actual reflector's position has been identified, while the mirrored reflector's position got averaged out.
In the following, we want to discuss a few of the most prominent reflectors highlighted in the map.
The reflector at A in the south has been identified as a shelter with a metal roof. The reflectors around B have been identified as a fence that is part of the airfield barrier.
Many reflectors have also been found close to the aircraft along the runway: The reflectors both in C and D are very persistent and have been identified as airport service buildings.
The reflectors in E, F, and G are most likely related to the railway lines along which they have been detected. Their origin might be the power lines along the railway, the rails or even a train that passed by during recording.
A few reflectors have been detected in the mining area north of the highway at H; the reflector in I is located between a parking lot and an intersection with traffic lights.
The reflectors in J and K are located in the suburbs north of the highway and are most likely some higher buildings.
East of the runway we have highlighted the building in L, as it is also the origin of a quite persistent reflector.

VIII. CONCLUSION AND OUTLOOK
In this paper we have described a processing chain for the evaluation of channel sounding data recorded during flight. We have explained our pre-processing and the actual MPC extraction process before introducing a path-based MPC tracking approach. In a last step, we have applied our processing chain to our measurement data and performed a reflector localization based on the detected and tracked MPCs. We have shown that our approach can be used to locate reflectors that are causing MPCs that have a significant impact on the wireless channel.
In a next step we want to evaluate more of our measurement data recorded during other flight scenarios and use the results for a reflector localization -depending on the scenario also in three rather than in just two dimensions -and finally for channel modeling.
The authors are also planning to apply the presented MPC tracking approach to measurement data collected during other campaigns, not necessarily limited to the aeronautical air-ground channel.

A. Calculation of the Impulse Response
The Impulse Response (IR) is computed by correlating the transmitted signal, here denoted by the reference signal x ref ∈ C N , and one instance of the channel sounding sequence in the measurement data, either denoted by y ∈ C BN or Y ∈ C N ×B . We implement the correlation in frequency domain, and assume a block processing of B consecutive impulse responses. 2 The computation of the IR either maps C BN → C Nf up ×B or C N ×B → C Nf up ×B , depending on the input data. In case of a vector input, y ∈ C BN is converted to a matrix Y ∈ C N ×B where the item in the n-th row and the m-th column y n,m is set according to y n,m = y[mB + n].
After this optional reshaping, the processing continues for both cases as follows. Assuming that x ref and y or Y , respectively, are sampled at a rate of f sr , the delay resolution of the IRs is given by

B Calculation of the Power Delay Profile
We denote the function computing the Power Delay Profile (PDP) by PDP The computation of the PDP is based on B consecutive Impulse Responses (IRs), thus we assume the IR matrix K has been computed according to Appendix A.
The complex-valued PDP is then computed by taking the mean of matrix K along its rows: a C :=K → .
The real-valued PDP is given by the complex-valued PDP's absolute value: a R := |K → |.
Please note that the definitions above correspond to the coherent PDP. Its delay resolution is given by (22).

C Calculation of the Doppler-Delay Spreading Function
We denote the discrete Doppler-Delay Spreading Function (DDSF), sometimes also called Cross Ambiguity Function (CAF), using reference signal x ref  First, the B Impulse Responses (IRs) are computed according to Appendix A, resulting in the IR matrix K.
Then, an FFT along the rows of K is performed: A := FFT B,→ {K}. We assume that this FFT includes an FFT-shift, such that the 0-th frequency is shifted to column B/2.
The columns of the resulting matrix A correspond to the delay axes, the rows correspond to the Doppler axes. While the resolution along the delay axis is the same as for the IRs as given in (22), the resolution along the Doppler axis is given by (23)

D. Distance Metrics
The function D : R 4 × R 4 → R returns the distance metric between two Multipath Components (MPCs) ξ n and ξ m .
The first step is to perform a feature scaling, as the typical values and units of the items of an MPC vary. If not denoted otherwise, we apply the popular min-max scaling to all elements x of an MPC: where x {min,max} denote a minimum/maximum value, respectively, that are used for the scaling of the specific item. They are chosen such that all appearing values of the specific item are mapped to a value within [0, 1).

1) Delay-Only Metric:
In this case, only the MPCs' delay is considered. This metric is useful in case no Doppler information is available and huge variations in the MPCs' amplitude are expected. We define it as where τ {n,m} are feature scaled (24) versions of τ {n,m} ∈ ξ {n,m} .
2) Delay-Doppler-Only Metric: In this case, the MPCs' complex weight is ignored. This metric is useful in case huge variations in the MPCs' amplitude are expected. We define it as D {ξ n , ξ m } = (τ n − τ m ) 2 + (ν n − ν m ) 2 ,

E. Two-Dimensional Reflector Localization
The delay and Doppler information of an MPC can be used to estimate the corresponding reflector's location. Here, we briefly describe the implementation of the reflector localization for the two-dimensional case (i.e. no height information is provided).
1) Exploiting Delay Information: All (potential) reflectors with the same delay are located on an ellipse whose focal points are given by the position of the transmitter (Tx) and the receiver (Rx). It is straight forward to compute this set of points based on a given delay τ : E τ = {P | TxP + P Rx = c air τ }, where c air denotes the speed of light in air.
If we assume an MPC-detection with a delay resolution of Δ τ , we can expand this to an inner and an outer ellipse. The reflector's position is then lying in the set: 2) Exploiting Doppler Information: The angle of arrival α with respect to the receiver's heading and speed v Rx can be estimated based on the Doppler shift ν according to: The reader may note the non-resolvable ambiguity caused by the nature of the cosine function. In case of a Doppler resolution of Δ ν , the set of points that arrive at the receiver within the given angle range is given by 3) Localization: The final position of the reflector is lying in the intersection of the sets E = E τ ∩ E ν . By combining the information gained through the processing of multiple blocks recorded over time, the uncertainty over the reflector's position can be reduced. 4) Exploiting Evolution: As it is assumed that an MPC's signal is caused by a reflector, its localization can be improved by observing the evolution of the corresponding MPC represented by a path γ. This can be done by performing the steps described above for each node of the path and finally joining these results: E γ = ∀b∈γ E (b) . 5) Implementation: In practice, the sets that are used to locate the reflector are represented by a two-dimensional matrix, where the rows correspond to the east-west axis and the columns correspond to the north-south axis of an East-North-Up (ENU) coordinate system with the transmitter in its origin. For each node of an MPC's path, the matrix is determined as given above.
Instead of computing the intersection of the sets as given in the last step, the exploitation of the evolution to reduce ambiguities is realized by computing the mean of the matrices of each path node. This approach is more robust against outliers as erroneous detections are averaged out.