Non-Stationary 3D M2M Channel Modeling and Verification With Aircraft-to-Aircraft, Drone-to-Drone, Vehicle-to-Vehicle, and Ship-to-Ship Measurements

Mobile-to-mobile (M2M) propagation channels have gained significant attention over the last years with the development of advanced communication systems for all kind of mobile stations such as aircraft, drones, cars, and ships. However, most available channel models do not account for the environment where the stations are located, but are defined for either average or worst-case conditions, not being able to predict the channel behaviour in specific scenarios. This is especially true for the scattering components of the channel, which are generally either ignored or defined as a rough extrapolation of the scattering components observed in other scenarios. In this work, we propose a geometry-based channel modeling technique that can be applied to any M2M scenario and that can calculate the channel accurately based on the environment around the stations. We first use finite and infinite planes to model the environment. Then, we use the proposed channel modeling technique to obtain analytically the contributions of each plane to the delay-dependent and joint delay Doppler probability density functions of the channel, as well as its squared delay/Doppler-spread function. Our technique focuses mainly on the scattering components but it also addresses the line-of-sight and specular reflection components. We apply the proposed channel modeling technique to different aircraft-to-aircraft, drone-to-drone, car-to-car, and ship-to-ship scenarios where channel measurements are available. In all scenarios, the channel estimated using the proposed channel modeling technique matches the channel measurements very accurately. Specifically, we observe that the scattering components are recreated very faithfully, and that we can even estimate how the channel evolves over time as the stations move and are affected differently by the environment.


Non-Stationary 3D M2M Channel Modeling and
Verification With Aircraft-to-Aircraft, Drone-to-Drone, Vehicle-to-Vehicle, and Ship-to-Ship Measurements Miguel A. Bellido-Manganell and Michael Walter , Senior Member, IEEE Abstract-Mobile-to-mobile (M2M) propagation channels have gained significant attention over the last years with the development of advanced communication systems for all kind of mobile stations such as aircraft, drones, cars, and ships.However, most available channel models do not account for the environment where the stations are located, but are defined for either average or worst-case conditions, not being able to predict the channel behaviour in specific scenarios.This is especially true for the scattering components of the channel, which are generally either ignored or defined as a rough extrapolation of the scattering components observed in other scenarios.In this work, we propose a geometry-based channel modeling technique that can be applied to any M2M scenario and that can calculate the channel accurately based on the environment around the stations.We first use finite and infinite planes to model the environment.Then, we use the proposed channel modeling technique to obtain analytically the contributions of each plane to the delay-dependent and joint delay Doppler probability density functions of the channel, as well as its squared delay/Dopplerspread function.Our technique focuses mainly on the scattering components but it also addresses the line-of-sight and specular reflection components.We apply the proposed channel modeling technique to different aircraft-to-aircraft, drone-to-drone, car-tocar, and ship-to-ship scenarios where channel measurements are available.In all scenarios, the channel estimated using the proposed channel modeling technique matches the channel measurements very accurately.Specifically, we observe that the scattering components are recreated very faithfully, and that we can even estimate how the channel evolves over time as the stations move and are affected differently by the environment.

I. INTRODUCTION
M OBILE-TO-MOBILE propagation channels have gained increasing attention over the last years with the de-velopment of mobile-to-mobile (M2M) data links in 5G and 6G systems.The peculiarity of such channels is that both, transmitter and receiver, are mobile.Compared to the fixed-to-mobile channels of traditional cellular systems, where only one station is mobile, M2M channels present rapidlychanging, non-stationary channel characteristics because of the movement of both stations.A general M2M channel model shall be able to adequately describe all known M2M channels, including the aircraft-to-aircraft (A2A), drone-todrone (D2D), vehicle-to-vehicle (V2V), 1 and ship-to-ship (S2S) channels.
Since the introduction of the wide-sense stationary uncorrelated scattering (WSSUS) assumption by Bello in [1], traditional narrowband fixed-to-mobile channels were generally assumed to behave like a WSSUS system and were often characterized by purely stochastic models like in [2].For uniformly distributed scatterers around the transmitter, and considering a stationary receiver, Clarke derived in [2] the wellknown Jakes power spectral density [3].For M2M channels, the Doppler spectral density was provided by [4], where they assumed that the scattering around the transmitter and receiver is uncorrelated, leading to a convolution of two Jakes spectra.More M2M channel models were derived in [5], [6], [7], [8] for various 2D and 3D scenarios, as well as for multiple input multiple output (MIMO) scenarios in [9], [10], [11], [12], [13].
The problem for most M2M channels, however, is that the WSSUS assumption is violated due to the movement of the transmitter and receiver.This has been observed in different M2M channels, including V2V [14], [15] and A2A [16] channels.The non-stationarity in channel modeling has been addressed in multiple ways.In [17], Matz generalized Bello's model to account for a non-WSSUS description of the channel consisting of fourdimensional channel correlation functions.The authors in [18] address the non-stationarity of the channel by accounting, among others, for the probabilistic presence and movement of clusters of scatterers.Based on their work in [19] and [20] for V2V and A2A channels, respectively, Walter et al. proposed in [21] a general way of analytically deriving the delay-dependent and joint delay Doppler probability density functions (pdfs) of M2M channels.The approach is based on defining an infinite arbitrarily-oriented scattering plane where the effective scatterers are located.This, together with the demonstration that the time-variant squared delay/Doppler-spread function remains proportional to the timevariant joint delay Doppler pdf [22], allows for an analytical description of the time-variant channel in any 3D M2M scenario where the scatterers can be considered to be uniformly distributed over a plane.The potential of such an approach has already been shown with measurement data in V2V [23], [24], [25] and A2A [26], [27] scenarios.However, the technique proposed in [21] has two key constraints: the scatterers can only be assumed to be distributed over a single plane, which must also be infinite.This prevents the technique to be applicable to more complex M2M scenarios where the scatterers or reflectors are present in multiple surfaces, which might also be limited in space.
The main new contributions of this work are the following.We adapt and extend the analytical technique proposed in [21] to 1) enable the use of any number of arbitrarily-oriented planes containing the scatterers, and to 2) allow the planes to be either infinite or finite, i.e., limited in space.This enables us to obtain the channel between two mobile stations by first recreating the environment around them using infinite and finite planes, and then by using the proposed analytical model to calculate the time-variant channel.This leads us to the next main contributions of this work.We apply our channel modeling technique to multiple M2M scenarios of interest where channel measurements are available.These include A2A, D2D, V2V, and S2S scenarios.We then compare the available channel measurements with the channel obtained using our channel modeling technique.Although the main focus of this work are the scattering components of the channel, we also account for the line-of-sight (LOS) and specular reflection components, which are studied more frequently in the literature compared to the scattering components.An excellent match between our results and the channel measurements can be observed, which validates our approach and motivates future work on this topic.We additionally show examples of how the computed pdfs can be used to derive some interesting characteristics of the channel, such as the channel auto-correlation and the time-variant delay and Doppler means and spreads.Among others, this work can additionally contribute to a better understanding of the M2M channels, which is key for realistic performance evaluations, e.g., [28], and for the development of 6G systems, e.g., [29], [30].It can also contribute to the research on intelligent reflecting surfaces, e.g., [31], as they could be recreated using finite planes like in our approach.
The remainder of this article is structured as follows.In Section II, we describe the proposed channel modeling technique and indicate how it can be implemented.The M2M scenarios of interest and the available channel measurements are discussed in Section III.We validate our approach by comparing its results with the channel measurements in Section IV.Finally, our conclusions are drawn in Section V.

A. Components of the Channel
We consider the channel to be composed of three type of components as proposed by Bello in [32]: the LOS component, the specular reflection components, and the scattering components.
First, the LOS component represents the signal propagated directly from the transmitter to the receiver.The LOS component, when present, will always arrive first, i.e., with the shortest delay, and with a Doppler frequency shift.In addition to the LOS component, the channel response is composed of multiple multi-path components (MPCs) arriving after the LOS component and caused by the reflection of the transmitted wave off different objects.We can distinguish between MPCs caused by specular reflections and by scattering.When a wave hits an object, a part of the wave is reflected in multiple directions.The wave component reflected in the specular direction, where the angle of reflection coincides with the angle of incidence, is denoted as the specular reflection component.In the other directions, the wave is scattered.Depending on the wavelength and on the roughness and electric characteristics of the reflector, the wave may be more focused in the specular direction, leading to a stronger specular reflection component, or scattered in all directions, leading to stronger scattering components.For example, a strong specular reflection component can be expected for a 1 GHz signal reflecting off calm water, e.g., a lake, whereas the same wave would be mostly scattered off vegetation, e.g., grass.
From the receiver point of view, the difference between the specular reflection components and the scattering components is clear.On the one hand, only few specular reflection components reach the receiver.They are, however, generally stronger than the scattering components and have discrete delays and Doppler frequency shifts.On the other hand, most objects around the receiver contribute to the scattering components received by it.Given the different geometries between the receiver and the surrounding objects, the scattering components are received with a wide range of delays and Doppler frequency shifts.

B. Using Planes to Recreate the Environment
In order to calculate the channel, we propose using planes to recreate the environment around the transmitter and the receiver, and to estimate analytically the contributions of the planes to the channel.Compared to previous publications where only one infinite plane could be used to obtain the channel, we can use any number of finite and infinite planes.This allows us to recreate any 3D M2M environment much more accurately and, thus, to obtain a better estimation of the channel.Depending on the M2M scenario to be considered, a different number of planes are required to recreate the environment accurately.For example, we could use a single plane to calculate the channel between two aircraft flying at a high altitude above a relatively flat field.However, more planes are needed if the aircraft fly at a low altitude with mountains nearby, as we show later in Section III.
We define the planes using a local Cartesian coordinate system (CCS) relative to the position of both stations, such that the n-th plane is defined as [21] A where l is half of the distance between the stations.Note that the planes can be defined initially in any arbitrary coordinate system and then be converted to our local CCS to calculate the channel.For example, the planes used to recreate the A2A scenario (Section III-A) are based on topographical data.The positions and velocity vectors of both aircraft are derived from the global navigation satellite system (GNSS).Thus, the data in the A2A scenario are initially defined in geographical coordinates, which we then convert to our local CCS in order to run our channel modeling technique and to calculate the channel.

C. Prolate Spheroidal Coordinate System
The description in Cartesian coordinates is useful to define the planes, but it turns out to be non-tractable to calculate the channel, since in this coordinate system the symmetries of the geometric setup are not exploited.Thus, we use a more advantageous coordinate system, which is called prolate spheroidal coordinate system (PSCS).The orthogonal surfaces of this particular coordinate system are ellipsoids, hyperboloids, and a half-plane.With the help of the ellipsoids centered around the transmitter (TX) and the receiver (RX), we can now exploit the geometry of the delay of single-bounce scattering.The coordinate transformation is given by [33] as We can now obtain the Doppler frequency of a component scattered off any arbitrary point.The Doppler frequency is interpreted as the spatial derivative of the total distance d sc (t, ξ) times the velocity vectors v t (t) = [v tx , v ty , v tz ] T and v r (t) = [v rx , v ry , v rz ] T of TX and RX, respectively.Note that the components of v t (t) and v r (t) are also time variant.However, we drop the time parameter to simplify the notation.The derivation of the Doppler frequency in PSCs can be found in [34].It is given by Thus, we obtained a full 3D Doppler frequency description in PSCs, where the normalized delay ξ is included in the description.By fixing ξ = ξ * a delay-dependent Doppler frequency can be easily realized.

D. Scattering Components
Given that the estimation of the LOS and specular reflection components is a well-studied field, we first focus on calculating the scattering components caused by the planes recreating the environment.
Our way of computing the scattering components is based on the methodology published in [21], where the authors propose to use one infinite plane to obtain the delay-dependent Doppler pdf and the joint delay Doppler pdf.For each delay, the authors compute the intersection ellipse q(η, ϑ) = 0 with the scattering plane, and obtain the density of scatterers either along the intersecting ellipse (for the delay-dependent Doppler pdf), or circumscribed by it (for the joint delay Doppler pdf).Assuming then that all scatterers are identical and uniformly distributed, their spatial distribution is transformed to obtain their Doppler frequency, leading to the computation of the pdfs.The transformation from the spatial domain to the frequency domain leads to ambiguities, which are solved by using the algebraic curve theory.
The approach proposed in [21] has two clear limitations: only one plane can be considered and it has to be infinite.This restricts the number of scenarios that can be recreated realistically.We generalize the methodology proposed in [21] by modifying and extending it to account for any arbitrary number of planes, which can now be either finite or infinite.Using more than one plane leads to the problem that the different planes might block each other entirely or partially.In other words, not all points of a plane might contribute to the scattering components received by RX, as either the TX-to-scatterer path or the scatterer-to-RX path might be blocked by another plane.A similar problem is faced when considering finite planes.The analytical expressions proposed in [21] assume an infinite plane where all points of the plane contribute to the scattering components.This allows for the analytical computation of the intersection ellipses, which are then used to normalize the delay-dependent and joint delay Doppler pdfs.However, in the case of a finite plane, the intersection ellipses (obtained analytically assuming an infinite plane) might not contribute fully to the scattering components, as some parts might actually be outside of the limits of the actual finite plane.
We adapt and extend the expressions given in [21] in the following way.In order not to make this manuscript excessively long, we do not provide all the unmodified expressions and explicitly refer the reader to [21] for some of them.The density of scatterers along the intersection ellipses is considered to be uniform and given by where is the accumulated effective length of all intersection ellipses, L n is the theoretical maximum length of the ellipse with the n-th plane, obtainable as in [21], and r l,n ∈ [0, 1] models the ratio of the intersection ellipse with the n-th plane that actually contributes to the scattering components, such that L n = r l,n L n is the effective length of the intersection ellipse with the n-th plane.
The delay-dependent Doppler pdf can then be obtained as where m 1,n represents the contribution of each plane as where and models if a specific point on a plane, defined by its (η n,k=j , ξ)or (ϑ n,k=i , ξ)-coordinates, contributes to the scattering components or not.First, δ in n,k ∈ {0, 1} indicates if a point is actually located within the n-th plane.This allows us to consider finite planes while maintaining the analytical expressions valid for infinite planes.Second, δ TX n,k ∈ {0, 1} and δ RX n,k ∈ {0, 1} model if a possible scatterer from the n-th plane is blocked by other planes either in the TX-to-scatterer path or in the scatterer-to-RX path, respectively.The factor | dη n,j df d |, the offset frequency f o , and the limiting frequency f lim can be obtained for each plane separately as in [21].Note that |H| ≤ 4 represents the up to four branches of the η (f d ) function.Likewise, when A n = 0 and B n = 0, the intersection ellipse becomes a circle and i = {1, 2} accounts for its two symmetric branches.
The joint delay Doppler pdf is computed as where m 2,n is weighted by the path loss and is obtained as and This is to be consistent with the notation used in [21], where the different notation accounts formally for the fact that, for the joint delay Doppler pdf, ξ is not fixed.Despite the different notation, the resulting factors used for the computation are the same.
The accumulated effective weighted area within all intersection ellipses is computed as where Y 1,n is the theoretical total area of the n-th intersection ellipse obtained as in [21], and r Y,n ∈ [0, 1] is the ratio of the area of the n-th intersection ellipse that contributes to the scattering components, such that Y n = r Y,n Y 1,n is the effective weighted area within the intersection ellipse.
Computing L T and Y T analytically is actually very complex, given that one has to obtain, for each plane, the ratio of each intersection ellipse that contributes to the scattering components, i.e., that is neither shadowed by other plane nor outside of the actual bounds of the plane (for finite planes).However, we can use the fact that L T and Y T are actually the normalization factors that guarantee that the delay-dependent and joint delay Doppler pdfs, respectively, integrate to 1.Then, one can simply estimate L T and Y T as the factors that normalize (7) and (12), respectively, taking into account the sampling interval used for ξ and f d .
In Algorithm 1, we show how the delay-dependent and joint delay Doppler pdfs can be computed efficiently.Note that depending on the parametrization that can be used for each plane, the computation of the contributions of each plane to the scattering components is based either on η (if A n = 0 or B n = 0) or on ϑ (otherwise).For simplicity, Algorithm 1 only shows the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.for f d do 6: for n-th plane do 7: Obtain end for 15: end for 16: Compute the delay-dependent and joint delay Doppler pdfs by using ( 7) and ( 12), respectively computation for the η-based parametrization, as it is the most common case.The computation for the ϑ-based parametrization follows the same steps but using the expressions for ϑ instead of those for η.Note that a more efficient implementation can be achieved if f d is defined as an array and array operations are performed.
In order to fully describe the channel in terms of its squared delay/Doppler-spread function, one must also account for the received power of the scattering components.Computing the power of the scattering components is, however, not trivial, and will be investigated for our channel modeling technique in future work.Nevertheless, according to [22], the joint delay Doppler pdf is proportional to the squared delay/Doppler-spread function.Thus, we can estimate the latter by scaling the calculated joint delay Doppler pdf.In this work, we use an arbitrary scaling factor Γ to compute the squared delay/Doppler-spread function from the joint delay Doppler pdf, in order to be able to compare the results of our channel modeling technique with the channel measurements.
It is important to see that the delay-dependent and joint delay Doppler pdfs are obtained for specific positions and velocity vectors of TX and RX.Thus, we calculate the channel for a specific instant, i.e., a snapshot of the channel.We can then change the positions of the stations according to their velocity vectors and calculate the channel again.This allows us to see how the channel changes as TX and RX move and the geometry between both stations and their environment varies, as we show later in Section IV.

E. LOS and Specular Reflection Components
The characteristics of the LOS component can be directly derived from the geometry between TX and RX, which are separated by a distance d los = 2 l.Thus, if no plane blocks the direct path between TX and RX, the LOS component is received after a delay τ los = d los /c with a Doppler frequency shift f d los and suffering an attenuation α los caused by the free-space path loss.
Each plane may cause a specular reflection of the transmitted wave.The relative delay of the specular reflection component off the n-th plane can be obtained as leading to an absolute delay τ n = d n /c = 2lξ n /c.However, RX receives a specular reflection component off the n-th plane only if certain conditions are met: ξ sr n ≥ 1, the point of reflection is within the actual plane boundaries (always the case for infinite planes), and neither the TX-to-reflector path nor the reflectorto-RX path is blocked by another plane.The specular reflection component off the n-th plane would be received with a Doppler frequency shift f d n and suffering an attenuation α n caused by the free-space path loss α p n over the travelled distance d n and by the reflection.The reflection coefficient can be obtained as where X = X h = e g n − cos 2 (θ n ) for horizontal polarization and X = X v = X h e gn for vertical polarization.The angle θ n between the n-th plane and the reflected signal can be obtained from the geometry.The complex relative permittivity e g n of the n-th plane is given by the ITU Recommendation P.527-4 for different surfaces and frequencies.
Each LOS or specular reflection component has a deterministic delay, Doppler frequency shift, and attenuation coefficient that only depend on the geometry between TX, RX, and the planes.Thus, each LOS or specular reflection component should be represented in the channel impulse response as a Dirac delta centered at its delay.Equivalently, they should be modeled as discrete points in the delay/Doppler-spread function.In the measurements shown in Section IV, we see that these components do not appear as discrete points in the measured squared delay/Doppler-spread function of the channel, but as sinc functions centered at the expected delay and Doppler frequency shift and stretching in the delay and Doppler directions.As we show in Appendix A, this representation is caused by the time-and bandwidth-limited sampling of the channel.For an easier comparison between the results of our channel modeling technique and the channel measurements, we also account for this effect by modeling the LOS and specular reflection components using (18).

F. Deriving Additional Properties of the Channel
The computed joint delay Doppler pdf provides a very complete description of the channel, as it shows how the channel components are distributed in the delay and Doppler domains and how such distribution evolves over time.Thus, the joint delay Doppler pdf can be used to directly obtain many other parameters of the channel, such as the mean delay and mean Doppler, as well as their spreads.Thanks to the time-variant description of the joint delay Doppler pdf, these parameters can also be computed for specific time instants in order to assess how they evolve over time as the transmitter and receiver move.They can also be obtained for specific channel delays, i.e., delay-dependent, for specific Doppler frequencies, i.e., Doppler-dependent, or for the entire channel.This enables a deep understanding of the channel that can be used to optimally design wireless links.For example, the delay-dependent Doppler mean and spread can be used for the design of an equalizer where the different delay taps are updated at different speeds.Given that the computation of the mean and standard deviation, i.e., spread, of a variable from its pdf is straightforward, we do not discuss the relations here but provide several examples in Section IV-E.
However, there are some properties of the channel that cannot be derived directly from the obtained joint delay Doppler pdf.For example, in order to estimate the fading properties of the channel, one would need both amplitude and phase information of each delay/Doppler component.This way, the constructive and destructive interference between the different MPCs can be recreated.Although we can estimate the amplitude information from the joint delay Doppler pdf because of its proportionality with the squared delay/Doppler-spread function, our method does not directly provide the phase information.Estimating it is not trivial, as one has to account for the phase shift caused by the pure wave propagation, as well as the one caused by the scattering for each component separately.In addition, the delay and Doppler resolution must be taken into account, and a stochastic distribution might be needed to compensate for it.Thus, extending our model to also obtain the phase information would require multiple assumptions and might still be cumbersome.In addition, its verification with measurement data would not be trivial because of the complexity of obtaining the fading properties from channel measurements, e.g., see [35].In order not to extend this manuscript excessively, we decided to leave such extension of our method as future work.

III. M2M SCENARIOS AND MEASUREMENTS
In order to verify the proposed M2M channel modeling technique, we consider different A2A, D2D, V2V, and S2S scenarios where channel sounding measurements are available.This way, we can recreate the different scenarios and use the proposed channel modeling technique to calculate the expected channel in each scenario.Each expected channel can then be compared with the measured channel obtained using the data collected in the different measurement campaigns.
The main parameters of the channel measurements can be found in Table I.In all cases, a Medav RUSK channel sounder [36] was employed to measure the time-variant channel every 1.024 ms (D2D and V2V measurements) or 2.048 ms (A2A and S2S measurements).Time and frequency synchronization was guaranteed through pre-and post-calibration, as well as by using GNSS receivers and atomic clocks when necessary.The accurate timestamps of the channel sounding data allow us to map the time-variant TX and RX positions and velocity vectors to the channel measurement data.
We initially use an arbitrary three-dimensional CCS to define each scenario.The CCS used to define the A2A and S2S scenarios corresponds to the Earth-centered, Earth-fixed (ECEF) coordinate system.This is very useful given that the GNSS-derived positions and velocity vectors of the stations can be easily converted to this coordinate system, as well as the topographical data used to define the scenarios.For the D2D and V2V scenarios, we use a local CCS for simplicity.Once the scenario is defined in the initial coordinate system, it is transformed to the local coordinate system defined in Section II in order to apply our channel modeling technique.Note that this enables adding new M2M scenarios without having to adjust the channel modeling technique.The only overhead of adding a new M2M scenario is the initial definition of the planes recreating the specific environment of interest in an arbitrary coordinate system, and the transformation from that coordinate system to the CCS defined in Section II.From that point on, the channel modeling technique is the same for any scenario and does not need to be adjusted.Thus, one can easily model other interesting scenarios such as aircraft-to-drone and drone-to-vehicle scenarios and obtain their channel using our technique.We do not include these scenarios in this work as no measurements are available to verify the results.

A. Aircraft-to-Aircraft (A2A) Scenarios
We consider three A2A scenarios of interest: flying over a field, flying through a valley, and flying over a lake.These scenarios were covered in the 2009 A2A channel measurement campaign [37] 1, where the aircraft positions in the three scenarios are highlighted.In all cases, both aircraft flew at an altitude of approximately 600 m above the ground and tried to keep a constant speed in order to maintain a distance between them of approximately 1.5 km.The Cessna aircraft transmitting the channel sounding signal followed the Dornier aircraft receiving the signal.The transmitting and receiving antennas were located below the aircraft fuselage.
1) Field: First, we focus on the flight segment where the aircraft flew over a field with a small hill nearby (see Fig. 1).As explained in Section II, we can recreate any scenario by employing finite and infinite planes to model it.In this case, we use one infinite plane to model the ground, i.e., the field, and 14 finite planes to recreate the hill located nearby.In order to accurately define the planes modeling this scenario, we use the topographical data available for the region, e.g., from Google Earth.The topographical data are sampled such that each scattering surface, e.g. the ground or a hillside, is well represented by a set of data triplets (latitude, longitude, and elevation).If the plane is finite, which is the case for all planes but the ground plane in this scenario, the set of data triplets must also enclose the surface.Then, the plane for each scattering surface is computed by finding the best-fitting plane for the set of data triplets representing the scattering surface.The resulting model of the field scenario can be seen in Fig. 2(a).Note that the process of modeling the landscape in the surroundings of the aircraft can be fully automatized as only the positions of the aircraft and the topographical data are required.However, as it is not the focus of this work, we have not fully automatized this process and we use the tools provided by the Google Earth application to identify which mountainsides might cause scattering, taking into account the expected area of operation of the aircraft, and to manually sample the landscape around the aircraft positions.
2) Valley: We consider the valley located in the Alps mountain range system and marked in Fig. 1.We use the topographical data to define a set of 38 finite planes modeling the mountains surrounding the valley and an infinite plane modeling the ground below the aircraft.The resulting model can be seen in Fig. 2(b).
3) Lake: We additionally study the scenario where two aircraft fly over a lake.As highlighted in Fig. 1, we consider a lake located within the Alps mountains.Again, we use topographical data to define 71 finite planes modeling the mountains and the ground around the lake.The resulting model can be seen in Fig. 2(c).Note that the reflection off a completely calm lake should only be specular.However, there were waves on the lake surface caused by the wind during the measurements, which we learnt to be more realistically modeled using a narrow plane below the aircraft and stretching in the direction of the water waves, i.e., in the wind direction.

B. Drone-to-Drone (D2D) Scenarios
We consider two D2D scenarios based on the channel measurement campaign conducted in 2019 between two drones in the premises of the German Aerospace Center (DLR) near Munich [38].One can see in Fig. 3 the model used to recreate the buildings (70 planes shown in gray) and grass fields (7 planes shown in green) around the location where the drones were flown.One can see that this environment is mostly urban, given that the drones flew near the surrounding buildings and below their rooftops.Compared to the A2A, V2V, and S2S model, the D2D model required significantly more planes to recreate the environment.This is understandable given that the stations are much closer to the reflecting objects in this case.Thus, the model of the environment must be able to recreate it more faithfully and to account for the irregularities of the reflectors, i.e., of the buildings, which could be neglected in the other scenarios because of the higher distance.
We consider the two following urban D2D scenarios.

1) Urban LOS:
The drones fly in the urban scenario (see Fig. 3) and maintain direct LOS between them.They fly at an altitude of 18 m (transmitter) and 15 m (receiver) above the ground, remaining below the altitude of the rooftops of the surrounding buildings.The receiver flies at 0.5 m/s towards the transmitter, which tries to maintain a constant position and altitude throughout the experiment.
2) Urban Without LOS: In this case, the LOS is blocked by a building as both drones fly away from the road intersection (see Fig. 3) at 0.5 m/s with an altitude of 10 m (transmitter) and 15 m (receiver) above the ground.

C. Vehicle-to-Vehicle (V2V) Scenarios
In the V2V scenario, the two vehicles drive on a road with forest on both sides.Three scenarios are considered: driving on the same lane in the same direction, driving on opposite lanes approaching each other in opposite directions, and driving on opposite lanes receding from each other in opposite directions.The three scenarios are recreated using the model shown in Fig. 4, which consists of two vertical finite planes modeling the forest lines at each side of the road.Note that we do not model the ground beneath the vehicles given that the antennas are located on their roofs and, thus, the ground reflections are blocked by the body of the vehicles.
The channel measurements were conducted on a road near Munich using a transmitting SUV G-Class Mercedes-Benz and a receiving van [39].Given that the vehicles drove in a straight line and maintained an approximately constant speed, we define the V2V scenarios using the coordinate system shown in Fig. 4, where the vehicles move only along the y-axis, i.e., v t (t) = [0, v ty , 0] T and v r (t) = [0, v ry , 0] T .This way, the definition of the V2V scenarios is easier compared to the A2A, D2D, and S2S scenarios.

D. Ship-to-Ship (S2S) Scenario
We also consider three S2S scenarios where two ships communicate with each other while sailing in the Heligoland archipelago in the North Sea at a speed of approximately 5 m/s.More information on the S2S channel measurement campaign can be found in [40].The model used to recreate the S2S scenarios is shown in Fig. 5.The mountain of the main island Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. of Heligoland is modeled using 13 finite planes (light green in Fig. 5), whereas the elevated buildings and constructions of both islands are modeled using 3 finite planes (light gray in Fig. 5).Note that neither the sea surface nor the flat fields are considered in the model, given that their contributions are expected to be negligible in the considered S2S scenarios compared to the reflections from the protuberant objects on land.Three scenarios are analysed: when the ships are sailing toward each other, when they are about to pass each other, and when they are receding from each other.The positions and velocity vectors of the ships for these cases are shown in Fig. 5.

IV. RESULTS
We now compare the results of our channel modeling technique with the measured channel in each scenario.We apply our channel modeling technique to each one of the scenarios defined in Section III and compute the expected time-variant squared delay/Doppler-spread function of the channel as described in Section II.Separately, we use the data of the channel measurement campaigns to obtain the measured time-variant squared delay/Doppler-spread function of the channel.In the channel measurement campaigns, the channel impulse response was measured every 1.024 ms (D2D and V2V measurements) or 2.048 ms (A2A and S2S measurements).In order to compute the delay/Doppler-spread function from the measurements, we must use multiple consecutive channel impulse responses (or snapshots of the channel).Logically, the resulting delay/Dopplerspread function is an average of the channel over the measuring window, which is not desirable in scenarios where the channel changes rapidly.However, it is a compromise between the desired Doppler frequency resolution and the averaging of the channel.We consider a measuring window of 2.097 s for the A2A, D2D, and S2S measurements and of 1.049 s for the V2V measurements.This entails using either 1024 or 2048 consecutive snapshots depending on the measurement time grid (see Table I).We use a shorter window for the V2V measurements given that the geometry in the approaching and receding V2V scenarios changes very quickly, and so does the channel, compared to the other scenarios.The D2D scenarios are also highly dynamic, but we consider a longer measuring window of 2.097 s to enhance the Doppler resolution.This is needed because of the low speed of the drones.
For the A2A and S2S scenarios, the averaging is unnoticeable, as the geometry practically does not change in 2.097 s.By contrast, the highly dynamic V2V and D2D scenarios lead to a quickly changing channel, and thus to a noticeable averaging effect.In order to be able to compare the results of our technique with the channel measurements, we recreated the averaging of the channel for the V2V and D2D scenarios.Given that our channel modeling technique calculates the channel at a specific time instant, i.e., snapshot, we apply it at multiple equally-spaced time instants within each measuring window and average the results.In other words, we apply our channel modeling technique to the different positions of the stations within each 2.097 s measuring window, and average the resulting snapshots of the channel.

A. Aircraft-to-Aircraft (A2A)
Let us first present the results for the A2A scenarios under test.Fig. 6 shows a comparison between the channel obtained using the proposed channel modeling technique in the A2A field scenario (upper figure) and the channel measurements obtained during the flight campaign (lower figure).In both cases, we see that the channel is composed mainly of two distinguishable components: the LOS component centered at 0 Hz and 3 μs and forming a cross in Doppler and delay directions, and the scattering components arriving after 5 μs.
As discussed in Section II, the LOS component should theoretically appear as a discrete point in the squared delay/Dopplerspread function.However, as shown in Appendix A, the timeand bandwidth-limited channel sampling entails a spreading both in delay and Doppler directions around the expected Doppler and delay shifts.This effect is mainly noticeable for the strong components, such as the LOS component.As we The scattering components arrive after a delay of approximately 5 μs and widen in the Doppler direction as the delay increases, forming the parabolic shape commonly seen in the literature.We can see that our channel modeling technique is also capable of accurately recreating the scattering components, including its starting delay, its shape, and its approximate distribution within its limiting frequencies.It is also interesting to see how our model accurately predicts that the power of the scattering components concentrates on the extreme Doppler frequencies as the delay increases, evolving toward the well-known Jakes spectrum.It is to be noted that the scattering components are mainly caused, in this scenario, by the field, i.e., the ground plane, given that the hill located nearby is not elevated enough compared to the aircraft altitude.
Note that some additional components can be seen in the channel measurements (lower figure of Fig. 6) before the LOS component, i.e., before 3 μs, especially at the limiting Doppler frequencies (±100 Hz).We neglect these components as they are artifacts caused by the periodic correlation of the channel sounding signal during the channel measurements.Thus, our channel modeling technique does not recreate such artifacts, as it can be seen in the upper figure of Fig. 6.
We have simulated the A2A valley scenario for different positions of the aircraft within the valley.Fig. 7 shows our results, compared with the measurements obtained in the flight campaign, for three different points within the valley.While Fig. 7(a) is obtained shortly after both aircraft enter the valley, Fig. 7(b) and (c) are obtained approximately 30 s later, with both aircraft flying in the middle of the valley.In the three figures we can recognise the strong LOS component centered at 5 μs and with a Doppler shift slightly above 0 Hz.After the LOS component, we can see the strong specular reflection component from the ground.This component, which was not observed in the results of the A2A field scenario because of the irregularity of the terrain, is centered at approximately 6.5 μs and presents a Doppler shift slightly below the LOS Doppler shift.In addition, it also spreads in delay and Doppler directions because of the time-and bandwidth-limited effect.We see that both, the LOS and the specular reflection components, are well recreated by the proposed analytical model.
After the LOS and specular reflection components, we see the scattering components coming from the ground and from the mountains surrounding the aircraft.Overall, the channel calculated by our channel modeling technique matches the measurements very accurately.Specifically, one can identify many interesting effects that are well recreated by our channel modeling technique.First, the external shape of the scattering components is practically identical to the one observed in the A2A field scenario and can be perfectly recreated by our technique.Second, the power density of the scattering components is particularly high shortly after the specular reflection, i.e., between 6.5 μs and 8 μs.This is actually caused by the geometry of the valley.The lower parts of the mountains and the area around the specular reflection point lead to many scattering components with very similar Doppler shifts and delays, which increases the power density in a narrow Doppler shift range immediately after the specular reflection component.Third, one can see that the proposed channel modeling technique is capable of accurately predicting the Doppler and delay regions where no scattering components are expected, as well as the isolated clusters of scattering components appearing at higher delays.In Fig. 7(a), the scattering components span most Doppler shifts (within the limiting Doppler frequencies) until a delay of approximately 15 μs.Afterwards, the scattering components are mainly focused on the limiting Doppler frequencies, i.e., around ±100 Hz, and only some isolated clusters can be seen, like the one observed in Fig. 7(a) between 23 μs and 25.6 μs with a Doppler shift between 20 Hz and 60 Hz.These isolated scattering components are also accurately recreated by our channel modeling technique, as well as the wider distribution of the scattering components at −100 Hz, compared to its narrower distribution at 100 Hz.
One can also see in Fig. 7(b) and (c) that our channel modeling technique recreates the channel very accurately.Not only the LOS and specular reflection components are well recreated, but also the shape of the scattering components and their distribution Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. in the delay/Doppler domain.For example, it could predict that the scattering components below 12 μs span most Doppler frequencies, whereas they are mainly concentrated around the limiting frequencies (±100 Hz) at higher delays.We can also highlight that the model is capable of recreating most of the isolated clusters of scattering components.For instance, the clusters of scattering components spanning from −40 to −80 Hz at 15-17 μs, and from 20 to 90 Hz at 17-22 μs, as well as the two small clusters spanning from 20 Hz down to −60 Hz at 21-25.6 μs.It is important to understand that the presence of the scattering components at certain Doppler frequencies and delays, or the absence thereof, is a deterministic process mainly driven by the geometry between the aircraft and their surroundings.Thus, our channel modeling technique is not only capable of predicting the channel accurately as shown here, but it can also identify which parts of the channel response correspond to specific reflectors, e.g., a specific mountainside, which can be very helpful for many applications.
Let us now briefly discuss the channel behaviour at different time instants.We can see some similarities and differences in the channel experienced by the aircraft when entering the valley (Fig. 7(a)) and in the middle of the valley approximately 30 seconds afterwards (Fig. 7(b) and (c)).First, the LOS and specular reflection components did not change radically, as the distance between both aircraft and to the ground was mostly maintained.In addition, the overall shape of the scattering and the limiting frequencies did not change significantly between both positions.By contrast, the distribution of the scattering components within its limiting frequencies changed noticeably from the beginning (Fig. 7(a)) to the middle of the valley (Fig. 7(b) and (c)).Comparing now Fig. 7(b) and (c), one can see that the channel did not change significantly, which is explainable by the fact that the geometry did not vary drastically between both points.Nevertheless, we can also notice that the isolated clusters of scattering components moved slightly in Doppler and delay directions.For example, one can see that the cluster at around 0 Hz and 22 μs moved toward lower delays and frequencies from Fig. 7(b) and (c).This movement of scattering components in the delay/Doppler domain, which can be seen for the other clusters as well, is to be expected because of the changing geometry between the aircraft and the surrounding reflectors as the aircraft fly.In fact, the slow movement of the scattering components led to the noticeable change in the channel between Fig. 7(a) and (b), given that the geometry changed slowly but consistently between both time instants.Importantly, we verified that the proposed channel modeling technique is capable of recreating the channel at any time instant and, thus, of tracking the channel changes as the transmitter and receiver move.
The channel obtained using the proposed channel modeling technique and through measurements in the A2A lake scenario is shown in Fig. 8 for three different positions of the aircraft above the lake: shortly after arriving to the lake (Fig. 8(a)), at the middle of it (Fig. 8(b)), and just before leaving it (Fig. 8(c)).In the three positions, the channel presents strong LOS and specular reflection components, which we could recreate very accurately.The presence of a strong specular reflection component was actually to be expected, given that the flat water reflects the signals mainly in the specular direction.In fact, the presence of the water explains most of the channel response seen in Fig. 8. First, we see that the power of the scattering components is high immediately after the specular reflection component, but decreases very rapidly.Afterwards, the scattering components off the lake are mainly concentrated around the limiting frequencies and are only sparsely distributed at intermediate Doppler frequencies.The fact that the scattering components off the lake are concentrated around the specular reflection component and on the limiting frequencies can be explained by the direction of the aircraft and of the lake waves, and is well represented in the channel generated by our channel modeling technique.After a certain delay, we see again some strong and uniform scattering components at intermediate frequencies, which are caused by the lake shore and the mountains surrounding the lake.This effect is especially noticeable in the upper figures of Fig. 8.In fact, it is reasonable to see that the scattering components appear first at negative Doppler frequencies in Fig. 8(a), given that the aircraft just entered the lake and left the lake shore immediately behind them, leading to negative Doppler frequencies.The scattering components with positive Doppler frequencies arrive from the opposite lake shore that the aircraft are flying toward.The opposite case can be seen in Fig. 8(c), where the aircraft are about to leave the lake.Likewise, one can see that the aircraft are halfway through the lake in Fig. 8(b) because the scattering components off the terrain appear at negative and positive Doppler frequencies from the same delay approximately.Overall, the proposed channel modeling technique is also capable of recreating quite accurately the scattering components, including their shape and distribution, off the terrain and mountains around the lake.
We now briefly discuss the way to model the lake surface.From the measurements observed in Fig. 8, we can see that the lake surface contributes to the channel response not only with the specular reflection component, but also with additional scattering components near the specular reflection component as well as near the limiting frequencies.In addition, we see that there is a higher concentration of scattering components at negative intermediate frequencies than at positive ones.We believe this to be caused by localised water turbulence flow in the southern part of the lake.Thus, given that the lake surface does not behave neither as a perfect mirror nor as a perfect scatterer, one might have to use a more specific model to recreate its waves.In this work, we opt for a more general solution, where we model the lake surface as a narrow plane in order to recreate its main contributions, i.e., the specular reflection component and the scattering components near the specular reflection component and at the limiting frequencies.We can observe in Fig. 8 that this model of the lake recreates the main lake contributions very accurately without being too scenario-specific.We leave a more complex modeling of the water surface for future work.

B. Drone-to-Drone (D2D)
Fig. 9 shows the squared delay/Doppler-spread function obtained using our channel modeling technique (upper sub-figures) and through channel measurements (lower sub-figures) for the two D2D urban scenarios of interest, i.e., with LOS (Fig. 9(a)) and without LOS (Fig. 9(b)).
It is important to first highlight some limitations of the channel measurements conducted in the D2D scenarios.First, the drones had to fly for safety reasons at a very low speed, i.e., below 1 m/s, due to the proximity to the buildings.Thus, the channel response is restricted to a low Doppler frequency range, as it can be seen in the measurements shown in Fig. 9. Given the low Doppler frequency range, we had to use a 2.097 s averaging window to enhance the Doppler resolution in order to distinguish the components of the channel and compare them with the ones calculated by our channel modeling technique.Second, the drones did not fly at a constant speed, but they adjusted their speed by accelerating and decelerating sporadically.This leads to spurious scattering components that are only present at certain delays and Doppler frequencies within very short periods of time.Not all of these components are recreated by our channel modeling technique, given that we use the instantaneous positions and velocity vectors logged by the sensors and GNSS receivers located on the drones, which are not capable of perfectly tracking the movement of the drones and suffer from the poor GNSS accuracy between the buildings.As a result, the measured D2D channel is generally blurrier than the one calculated using our channel modeling technique.Let us now compare the channel measurements with the channel obtained using our channel modeling technique in the D2D scenarios.In the channel measurements shown in Fig. 9 The LOS component is not present in the non-LOS scenario, as it is blocked by a building.We can see that our channel modeling technique was capable of predicting and recreating these phenomena accurately.In the LOS scenario (Fig. 9(a)), the scattering components arrive immediately after the LOS component, and spread in the negative Doppler direction as the delay increases, reaching a minimum Doppler frequency of −10 Hz.The scattering components decrease in power as the delay increases, becoming almost negligible after 400-500 ns for most Doppler frequencies.However, in the non-LOS scenario (Fig. 9(b)), some scattering components arrive before the specular reflection component and drift toward higher Doppler frequencies while losing power as the delay increases.One can observe that our channel modeling technique recreates the scattering components in both scenarios faithfully, as it is capable of reproducing their starting delay, their shape, their limiting frequencies, and how the power is distributed in the delay/Doppler domain.As discussed before, a blurrier channel response from the channel measurements was expected given the limitations of the D2D measurements.In addition, we have not modeled the smaller objects that were present during the D2D measurements, such as measuring equipment, bushes, or benches.We opted not to model them in order to show that our channel modeling technique can already recreate most of the channel effects faithfully, even if a simple model of the environment is used.

C. Vehicle-to-Vehicle (V2V)
Fig. 10 shows the channel calculated using the proposed channel modeling technique (upper row) and through channel measurements (lower row) in the three V2V scenarios.
In the first V2V scenario (Fig. 10(a)), where the vehicles follow each other on the same lane at a constant speed, the channel measurements show a strong LOS component at 155 ns.The scattering components arriving after the LOS follow the parabolic shape that we already saw in the A2A measurements, i.e., they widen symmetrically in Doppler domain as the delay increases.In this case, however, the scattering components are mainly present at the limiting frequencies, and are mostly negligible for intermediate frequencies after a short delay.The scattering components are caused by the forest lines on both sides of the road.Given that the vehicles are closer to one of the forest lines, one can actually see that the scattering components are distributed in two overlapping parabolic forms arriving at slightly different delays, each one caused by one forest line.We can see that the channel modeling technique manages to recreate the channel very accurately, including the LOS and scattering components.The characteristic distribution of the scattering components is also well recreated, including the contributions of both forest lines starting at different delays and with slightly different shapes.However, we can also see that the scattering components are also present in intermediate frequencies for delays of 200-400 ns, which is not reproduced by our model.This is actually to be expected, given that we only modeled the two forest lines and did not consider other potential scatterers such as the terrain between the forest lines and the road.In addition, we consider a very simple model for the forest line and neglect its many irregularities, such as the different trees and their branches.Modeling the environment in such depth is not worth the additional effort, as the simple model considered here already allows us to recreate the main channel effects accurately.
In the second V2V scenario (Fig. 10(b)), both vehicles are driving toward each other on opposite lanes.Thus, the geometry between both vehicles changes within the measuring window as they get closer, which leads to the averaging over time that we see in the results shown in Fig. 10(b).One can see the very good match between the expected and the measured channels.In both cases, we observe that the LOS component spans in delay from 300 ns up to 350 ns because of the movement of the vehicles.The scattering components decrease in frequency as the delay increases, converging into two separate tails close to 0 Hz, each one caused by a forest line.One can see that this is faithfully recreated by our channel modeling technique.
Finally, we see in Fig. 10(c) the channel in the third V2V scenario, where both cars are receding from each other on opposite lanes shortly after passing by.We can actually see that the channel is practically a mirrored version for negative frequencies of the channel observed in the second V2V scenario.Again, an excellent match between the measurements and the channel calculated using the proposed channel modeling technique is achieved.

D. Ship-to-Ship (S2S)
We see in Fig. 11 the channel obtained using our channel modeling technique and through channel measurements in the three S2S scenarios shown in Fig. 5.
In all three cases we observe a very strong LOS component at 4.2 μs, 0.6 μs, and 2.4 μs respectively in each scenario.As expected, the Doppler frequency shift of the LOS component decreases from positive frequencies in the first S2S scenario (Fig. 11(a)), where the ships sail toward each other, down to negative frequencies in the third S2S scenario (Fig. 11(c)), where the ships sail away from each other.In all three cases we see that the LOS component, as well as its spreading in delay and Doppler directions, is perfectly recreated by our channel modeling technique.
One can also notice the scattering components arriving 1-2 μs after the LOS component.In all three S2S scenarios, there is a very good match between the scattering components estimated by our channel modeling technique (upper row) and the measured scattering components (lower row).As expected, the scattering components change as the ships sail in the Heligoland archipelago, given that the geometry between the ships and the islands also changes.
It is to be noted that we did not model the Heligoland islands entirely (see Fig. 5), but only the mountains (only present in one of the islands) and the regions where objects stand out, such as elevated constructions.Thus, we neglected both the sea around the ships and the non-prominent terrain, given that their reflections are expected to be mainly specular and not reach the receiver.This decision has proven correct in the considered S2S scenarios given the excellent match between the results obtained using our model and the channel measurements.Therefore, we can conclude that most scattering components, if not all, came from the protruding objects of the islands, and were recreated very accurately by our channel modeling technique using the model shown in Fig. 5.

E. Exemplary Channel Characteristics Computed From the pdf
As discussed in Section II-F, many characteristics of the channel can be derived from its joint delay Doppler pdf obtained using our channel modeling technique.It is to be noted that computing all possible parameters and comparing them with measurements in each scenario would make this work too lengthy, and shall be done in separate publications.For example, we analyze in [35] the fading characteristics of the A2A channel based on measurements, which can be used as basis for a later comparison.However, we provide here several examples of the characteristics of the channel that can be computed using the obtained joint delay Doppler pdf.
First, we consider the A2A valley scenario and use one of the obtained joint delay Doppler pdf of the scattering components to compute the delay-dependent Doppler mean and spread, shown in Fig. 12 together with the pdf of the scattering components.The mean delay and mean Doppler for all delays and Doppler frequency shifts are also shown in Fig. 12 as a red asterisk, with the bounds around it depicting the computed delay and Doppler spreads.Note that we do not consider the LOS and the specular reflection components here for simplicity, but can likewise be considered as shown in Fig. 7(a).One can clearly see in Fig. 12 that the main scattering components are present at low delays with a moderate Doppler spread, which increases with the delay as the scattering components concentrate around the limiting frequencies.
We can also obtain the normalized channel auto-correlation, related to the characteristic function of the channel, by doing an inverse Fourier transform in Doppler direction of the joint delay Doppler pdf shown in Fig. 12.The resulting auto-correlation function, shown in Fig. 13, gives information about the channel stability.For example, Fig. 14 depicts different cuts of the computed auto-correlation function at different delays, which show how the channel becomes less stable as the delay increases.At high delays, the spectrum is practically identical to a theoretical Jakes Doppler spectrum, which leads to an auto-correlation function matching the Bessel function.This well-known relation [2],    15.Evolution of the mean delay, delay spread, mean Doppler, and Doppler spread computed using the joint delay Doppler pdf of the scattering components in the V2V crossing scenario (Fig. 4) as both vehicles approach each other (positive distances) and recede from each other (negative distances).[41] was actually demonstrated analytically in [16] for a single infinite plane building the environment.
We now consider the V2V crossing scenario and obtain the joint delay Doppler pdf, again only for the scattering components for simplicity, at different time instants as the vehicles move toward each other, cross each other, and finally recede from each other.Using the computed pdfs, we derive the time-variant mean delay and time-variant mean Doppler shift of the scattering components, as well as their time-variant spread in delay and Doppler directions.The obtained results can be seen in Fig. 15.As expected, the mean channel delay decreases as the vehicles approach each other and increases as they recede from each other, reaching a minimum at the crossing point.Likewise, the mean Doppler shift drifts from positive values toward negative values, as expected already from the channel shown in Fig. 10(b) and (c).The Doppler spread reaches a minimum when the vehicles cross each other because of the shift of the spectrum from positive to negative frequencies.

V. CONCLUSION
We propose a channel modeling technique for M2M propagation that employs a model of the environment to obtain the time-variant channel between two mobile stations.We apply the proposed channel modeling technique to multiple A2A, D2D, V2V, and S2S scenarios, and compare the obtained results with the channel measurements collected in different channel measurement campaigns.In all scenarios, we observe that the proposed channel modeling technique is capable of recreating the time-variant channel very accurately, including the LOS, specular reflection, and scattering components of the channel.As expected, the best match is observed in the scenarios where the scatterers are well-described by planes, i.e., when their irregularities are negligible compared to the overall environment.In addition, this technique allows us to see how the channel evolves as the stations move.Moreover, we can use this technique to identify which reflectors or scatterers cause specific responses in the delay/Doppler domain.

APPENDIX A EQUIVALENT LOW-PASS CHANNEL MODEL
The time-variant equivalent low-pass channel transfer function for L MPCs is given by α r l (t)α p l (t)e −j2πf c τ l (t) e −j2πf τ l (t) .(17) In the channel measurements, the transfer function is sampled in time domain at t = mΔt with m = 0, 1, . .., M − 1 for an observation window T o = M Δt, and in frequency domain at N points f = nΔf within a bandwidth B such that Δf = B/N .If we use a Taylor series expansion of the delay and approximate it until the linear term in the first exponential and until the constant term in the second exponential in (17), we can separate t and τ .Then, by considering that ∂τ l (t)/∂t| t=t 0 f c = −f d,l and using t 0 = 0 without loss of generality, one can demonstrate that the delay/Doppler-spread function, simplifying for constant αl and τ l , is given by with n = − N −1 2 , . .., 0, . .., N −1 2 and αl (t) = α r l (t)α p l (t)e −j2πf c τ l (0) .One can notice in (18) the sinc functions affecting the delay/Doppler-spread function in both k (delay) and q (Doppler) directions, which can be seen for the strong MPCs in the measurements shown in Section IV.

Algorithm 1 :
Compute the Delay-Dependent and Joint Delay Doppler Pdfs.1: Compute positions and velocity vectors of TX and RX 2: Define planes, i.e., (A n , B n , C n , D n ) for n = 0, 1, . .., N − 1, including the boundaries of the finite planes 3: Define range of ξ and of f d to be investigated 4: for ξ do 5:

Fig. 1 .
Fig. 1.Flight route of the A2A channel measurement campaign in southern Germany.The three A2A scenarios of interest are highlighted.Copyright of background map: Map data ©2021 GeoBasis-DE/BKG (©2009), Google.
using a Cessna C-208B (D-FDLR) and a Dornier 228-101 (D-CODE) aircraft to measure the A2A channel at 250 MHz.Part of the flight route can be seen in Fig.

Fig. 2 .
Fig. 2. Models of the A2A scenarios.The topographical data are sampled and then used to calculate the best-fitting infinite and finite planes modeling the topography.We show the planes here, together with the instantaneous aircraft positions and velocity vectors in the middle of each scenario.One can see the multiple finite planes (dark green) modeling the hills or mountains.An infinite plane (light green) models the ground in (a) and (b).The lake is modeled in (c) using a narrow finite plane (blue).(a) A2A field scenario.(b) A2A valley scenario.(c) A2A lake scenario.

Fig. 3 .
Fig. 3. Model of the urban D2D scenarios.The instantaneous drone positions and velocity vectors in the middle of the scenarios are shown.WGS84 coordinates of the site: 48.084559 • latitude and 11.278297 • longitude.

Fig. 4 .
Fig. 4. Model of the V2V scenarios.The instantaneous vehicle positions and velocity vectors in the middle of the second scenario, i.e., opposite lane approaching, are shown as an example.

1 )
Same Lane Following: In this case, the transmitting vehicle follows the receiving vehicle maintaining a distance of d los = 46 m.Both vehicles drive at an approximately constant speed of v t = v r = [0, 32.4,0] T km/h.2) Opposite Lanes Approaching: The vehicles drive now on opposite lanes and approach each other, as shown in Fig. 4. The distance between the vehicles decreases in this scenario from d los = 105 m down to d los = 90 m.In this case, v t = [0, 32.1, 0] T km/h and v r = [0, −29.3, 0] T km/h.3) Opposite Lanes Receding: Finally, the vehicles are receding from each other on opposite lanes with a speed of v t = [0, 32.1, 0] T km/h and v r = [0, −29, 0] T km/h and separated by a distance increasing from d los = 85 m up to d los = 100 m.

Fig. 5 .
Fig. 5. Model of the S2S scenarios.One can see the instantaneous ship positions and velocity vectors for the three cases of interest, i.e., approaching, shortly before passing by, and receding.

Fig. 6 .
Fig. 6.Channel obtained using the proposed M2M channel modeling technique (upper figure) and through measurements (lower figure) in the A2A field scenario shown in Fig. 2(a).The LOS and scattering components are recreated very accurately by our channel modeling technique.

Fig. 7 .
Fig. 7. Channel obtained using the proposed M2M channel modeling technique (upper row) and through measurements (lower row) at three positions of the aircraft in the A2A valley scenario shown in Fig. 2(b).The channel changes as the aircraft move throughout the valley ((a) vs. (b)) but it does not change very quickly ((b) vs. (c)) because of the slowly changing geometry.Our M2M channel modeling technique can track the channel changes.(a) 14:19:00.00to 14:19:02.10UTC.(b) 14:19:29.36to 14:19:31.46UTC.(c) 14:19:31.46to 14:19:33.55UTC. can see in Fig. 6, our channel modeling technique is capable of recreating the LOS component accurately, as well as the timeand bandwidth-limited effect, which is clearly visible both in the expected (upper figure) and in the measured (lower figure) squared delay/Doppler-spread functions.The scattering components arrive after a delay of approximately 5 μs and widen in the Doppler direction as the delay increases, forming the parabolic shape commonly seen in the literature.We can see that our channel modeling technique is also capable of accurately recreating the scattering components, including its starting delay, its shape, and its approximate distribution within its limiting frequencies.It is also interesting to see how our model accurately predicts that the power of the scattering components concentrates on the extreme Doppler frequencies as the delay increases, evolving toward the well-known Jakes spectrum.It is to be noted that the scattering components are mainly caused, in this scenario, by the field, i.e., the ground plane, given that the hill located nearby is not elevated enough compared to the aircraft altitude.Note that some additional components can be seen in the channel measurements (lower figure of Fig.6) before the LOS component, i.e., before 3 μs, especially at the limiting Doppler frequencies (±100 Hz).We neglect these components as they are artifacts caused by the periodic correlation of the channel sounding signal during the channel measurements.Thus, our channel modeling technique does not recreate such artifacts, as it can be seen in the upper figure of Fig.6.We have simulated the A2A valley scenario for different positions of the aircraft within the valley.Fig.7shows our results, compared with the measurements obtained in the flight campaign, for three different points within the valley.While Fig.7(a) is obtained shortly after both aircraft enter the valley, Fig. 7(b) and (c) are obtained approximately 30 s later, with both aircraft flying in the middle of the valley.In the three figures we can recognise the strong LOS component centered at 5 μs and with a Doppler shift slightly above 0 Hz.After the LOS component, we can see the strong specular reflection component from the ground.This component, which was not observed in the results of the A2A field scenario because of the irregularity of the terrain, is centered at approximately 6.5 μs and presents a

Fig. 8 .
Fig. 8. Channel obtained using the proposed M2M channel modeling technique (upper row) and through measurements (lower row) for three positions of the aircraft in the A2A lake scenario shown in Fig. 2(c): Shortly after arriving to the lake (a), at the middle of it (b), and just before leaving it (c).The channel changes as the aircraft fly over the lake and our M2M channel modeling technique can recreate the channel accurately at all positions.(a) 14:24:26.79to 14:24:28.89UTC.(b) 14:24:39.37 to 14:24:41.47UTC.(c) 14:24:51.95to 14:24:54.05UTC.
(a) and (b) (lower sub-figures), one can observe a strong component standing out above the scattering components.In the LOS scenario (Fig. 9(a)), the dominant component at 120 ns corresponds to the LOS component, whereas in the non-LOS scenario (Fig. 9(b)), the dominant component at 190 ns is a strong specular reflection component off a building facade.

Fig. 12 .
Fig. 12. Computed joint delay Doppler pdf of the scattering components in the A2A valley scenario (Fig. 7(a)).The joint delay Doppler pdf is used to obtain the delay-dependent Doppler mean μ Doppler (τ ) and spread σ Doppler (τ ) shown here.The mean delay and mean Doppler for all delays and Doppler frequency shifts are shown as a red asterisk ( * ) with the bounds depicting the spread in the delay and Doppler directions.

Fig.
Fig.15.Evolution of the mean delay, delay spread, mean Doppler, and Doppler spread computed using the joint delay Doppler pdf of the scattering components in the V2V crossing scenario (Fig.4) as both vehicles approach each other (positive distances) and recede from each other (negative distances).

TABLE I PARAMETERS
OF THE CHANNEL MEASUREMENT CAMPAIGNS