Distributed Communication and Sensing System Co-Design for Improved UAV Network Resilience

The progress in wireless technology over the past decade led to the rapid adoption of unmanned aerial vehicles (UAVs) for various applications. As the interest in UAVs is accelerating, increased attention is paid to the reliability and resilience of UAV-based systems with respect to the collision avoidance. One of the ways to improve this aspect is to utilize the RADAR functionality. In this work, we consider a cellular network employed for communication jointly with RADAR operation. The critical parameter that affects the RADAR algorithm is the radar cross-section (RCS). Since the task of obtaining the RCS of a complex-shaped object is extremely challenging, we first propose a novel, accurate, and fast method of scattered field assessment. We further perform radio network planning for the cellular deployment, as well as link budget estimations for the RADAR system that co-exists with it. Under this system model, we carry out detailed Monte-Carlo simulations of the RADAR detection process to obtain reliable statistical results and answer the question of how the actual bistatic RCS model affects the detection algorithm. We then apply mathematical modeling based on stochastic geometry to estimate the collision probability without the need to simulate an extensive number of flight-hours. Our numerical results confirm the robustness to RCS pattern nulls, which is crucial for safety-centric applications such as collision avoidance.

However, prior to the mass adoption of UAVs within densely populated urban areas, their safety must be ensured. To this end, regulators attempt to assure safety of UAV usage in terms of individual reliability as well as their ability to avoid collisions with each other.
By design, the payload of UAVs is limited. Hence, any sensing equipment such as LIDARs or cameras for collision avoidance cuts into payload significantly. Therefore, the steps toward assurance of UAV safety in an urban environment should, if possible, be based on the integration of sensing and communication capabilities [2]. Employing communication signals from terrestrial cellular network jointly with sensing allows to save on the mass, cost, control circuitry, and power consumption that would be required to build a separate sensing system.
There are various proposals reported in recent literature to codesign and analyze the converged sensing and communication functionalities within city infrastructure that target application of RADAR for UAV collision avoidance. As a feasible approach to this challenge, studies address the deployment of multiple monostatic RADARs on the facades of buildings and/or street lamps to ensure the coverage of the area between buildings [3]. As another more complex example, there are proposals to exploit a dynamic RADAR network (DRN) composed of UAVs able to intelligently adapt and share the sensed information with their neighbors via multi-hop networks [4], [5]. Apparently, both of these options require dedicated spectrum resources, as well as massive investments in either infrastructure or communications protocol design.
The coexistence of communication and RADAR systems has been extensively studied over the past decade, with a focus on developing efficient interference management techniques, such that the two individually deployed systems can operate concurrently [2], [6]. However, by utilizing a converged sensingand-communication approach, one should be able to arrive at a more balanced solution. We expect that in a city with dense cellular infrastructure, there are plentiful signals transmitted by the base stations (BSs) that could, in principle, be utilized for sensing. It has been demonstrated that it is indeed possible to use conventional 5G waveforms for localization [7]. In our study, we consider a system where BSs on the ground form a distributed RADAR transmitter (TX), while UAVs are receivers (RXs) only. Such a system has multiple advantages, as it requires no extra infrastructure other than more advanced RX on the UAVs. By combining signals from several TX sites over time, one can perform multiple perspective observations to improve This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ spatial resolution [8]. Despite these benefits, as a standalone system, distributed RADAR may be prohibitively expensive. In our envisaged solution, however, it is seamlessly integrated with the cellular infrastructure at minimal additional expenses.
While simple in principle, numerical evaluation of such system's performance requires the detailed knowledge of radar cross-section (RCS) [9], [10] for accurate results. Hence, in addition to the system design, we propose a new method for accurate estimation of the RCS for curved objects with complex shapes, such as drones. This further allows us to answer the question of how the accuracy of the RCS model affects the converged sensing and communication system's performance.

B. Delivered Contributions
First, in Section II, we consider system design and model co-design for an integrated cellular/RADAR scenario. We perform radio network planning for the cellular deployment, as well as aim to conduct link budget estimations for the RADAR system that co-exists with it.
Further, to perform qualitative analysis, we need an accurate and robust method for high-precision calculation of the reflected electric field from objects with a complex shape (i.e., UAVs). This is necessary to compute the bistatic RCS, which affects how well does a UAV reflect a signal toward a certain RX point when illuminated from a certain TX point. It is important to note that most works on UAV detection tend to avoid the need for accurate RCS by focusing on Doppler patterns of the UAVs [11], [12].
As the calculation of a bistatic RCS pattern for a complex object is extremely challenging [13], many researchers resort to using a constant RCS value for their simulations [5], [14], or abandon the bistatic RADAR setup entirely. In Section III, we demonstrate a new method that achieves highly accurate bistatic RCS computation with adequate computational complexity, and thus becomes well-suited for large-scale simulation studies. We further illustrate how it can be used to compute 3D channel impulse response (CIR), or other desired parameters of interest.
Further, in Section V, we perform a system-level mathematical modeling based on stochastic geometry to obtain the estimates of collision probability without simulating the extensive volumes of flight-hours. This helps understand which parameters of our model are the most impactful w.r.t. the overall performance in terms of the final collision probability. Further, in Section VI, we embed the proposed scattering model into the system model and conduct detailed Monte-Carlo simulations of the RADAR detection process to obtain reliable statistical results. Finally, we assess the results and discuss possible extensions of the proposed methodology. We present the essential acronyms in Table I. The main parameters and the description of employed notation are collected in Table II.

II. SYSTEM DESIGN CONSIDERATIONS
In this section, we introduce our key system design considerations. Overall, we seek to capture the case where multiple UAVs are aiming to leverage the existing cellular network deployment  I  SUMMARY OF ACRONYMS   TABLE II  NOTATION OF THIS WORK for sensing purposes. To this end, we introduce the basic assumptions of the considered system. First, we assume that a 4G/5G network is deployed in the area where UAVs are operating. For this first-order evaluation, we require the spectrum to be non-overlapping because we use the sounding reference signal (SRS), which features Zadoff-Chu sequences [15] as a reference signal for RADAR processing: each of the BSs has its own orthogonal sequence. We further assume that the interference between the adjacent BSs is on that level that there is a possibility to distinguish different BSs using code division. Because of the high robustness of these sequences to radio interference, one does not need to take into account the interference in the simulation of sensing operations. Our proposed system may operate in either 4G or 5G environment, as long as BS antenna has acceptable amount of radiation toward the UAVs and does not employ beacon beamforming. For the sake of exposition, however, below we only consider 5G systems.
Importantly, the flight of UAVs is restricted in altitude to a narrow range. This is justified by the existing regulations on the maximum UAV altitude, which are present e.g., across the EU [16], and a reasonable expectation that drones would need to reside sufficiently above the buildings to be able to maneuver in emergencies. Hence, we do not consider buildings as possible collision threats. The UAVs in question are expected to be large, bulky devices incapable of abrupt turns; for analysis purposes, we assume that the UAVs fly along straight lines, unless performing collision avoidance, and that their flight paths are therefore random.

A. Detailed Network Model
We now discuss the 5G network model in detail, and focus on specific parameters that impact the numerical results. For reference, an overview of the considered deployment is illustrated in Fig. 1.
The terrestrial BSs are assumed to follow cellular structure in 2 with inter-site distance of L, thus, resulting in the deployment density of λ B units/km 2 . Regarding the height of the BSs, we argue that there are plenty of buildings whose height is around 30 m for BS deployment even in the cities seeing an average height of 10 m. Moreover, the LTE BS heights range from approximately 23 m to 32 m, with the average height being 27.81 m. Therefore, the height of all the BSs h A is required to be similar and, on average, equal to 30 m. The corresponding "BS plane" in Fig. 1 is denoted by A 0 . There are also uniformly distributed square buildings on the ground with the density of λ b units/km 2 . The height and width of the buildings are denoted as h b = 12 m and w b = 30 m, correspondingly. We utilize the same object model for all buildings in the scene, however, they have random orientation.
The altitude of all the UAVs is constrained within the allowed limits (see above), thus, placing them between the planes A 1 and A 2 . The height under A 0 of the lower UAV flight level is set to h L = 10 m. The height of the permitted UAV flight region is set to h U = 25 m. The approximate speed of all the UAVs (V UAV ) is 20 mps. However, this value may vary because of the changing direction of the UAVs. Therefore, the approximate maximum Doppler shift for our system is equal to ±100 Hz. We assume that the speed of the UAVs is much lower compared to the measurement window of any RADAR. This allows to consider a fixed time instant t in our model for analysis purposes. The trajectories of the UAVs are assumed to form a Poisson line process in 2 . Hence, at any given time instant t, the UAVs constitute a homogeneous Poisson point process (PPP) in 2 with the density of λ U units/km 2 .
To assess the system performance, we randomly choose one UAV ("UAV A" in Fig. 1), at which the RX is located. All other UAVs are thus designated as targets, and in what follows, we are interested in characterizing the detection performance of UAV A w.r.t. all possible targets around it.

B. Radio Part Specifics
In a RADAR-centric system, to improve the detection probability, one may prefer to use the lowest frequency band possible. Hence, we consider the lowest available cellular frequency of F c = 800 MHz, which would normally be used in the networks to ensure that there are no gaps in coverage and/or MTC.
1) Antenna Model: Where possible, we follow the 3GPP recommendations from TR 37.840 [17]. In particular, a crucial component of our system is the BS antenna model, G A ( − → d ), which specifies the gain as a function of the signal propagation direction − → d . In our scenario, the direction of arrival − → d OA and the direction of departure − → d OD are illustrated in Fig. 1. The antenna patterns are modeled as where w m,n , v m,n are the weighting factors, while N H and N V are the numbers of antenna elements in horizontal and vertical dimensions. The first term in (1), , is a single element pattern given by where G E is the maximum directional gain, which is set to 8 dBi, while A E H and A E V are the values of attenuation in horizontal and vertical planes, and A m = 30 dB is the front/back ratio. The reception pattern of the UAV RX antenna is modeled similarly.
To capture the radio propagation for radio network planning purposes, we utilize the 3GPP urban micro (UMi) model [18]. With the above considerations in mind, and for the purposes of a numerical evaluation, we select a generic target inter-site distance of L = 350 m, TX power P A of 40 dBm, antenna downtilt of 10 degrees, and frequency reuse scheme 1/3/3. The values for w m,n , v m,n , N H , and N V are then adjusted through the parameter search to match the desired average SINR of 20 dB within each cell subject to reasonable fluctuations. The resulting distribution of the RX power is demonstrated in Fig. 2. However, our key findings are expected to hold across a wide range of system parameters.
The UAV antennas affect the performance of the target system, and an antenna with high directivity is preferable to reduce the parasitic reflections from the ground clutter, as well as from objects that are behind the UAV. Therefore, for the considered in this section antenna model, we use such parameters that make it highly directional. However, an improved UAV antenna may also provide better performance in the proposed system.
2) Propagation Model: For UAV RADAR path loss calculation, one can utilize the free-space model, since the events of interest occur above the rooftop level, and thus multipath is insignificant. In our analysis, we employ the free-space model in its linear form ( 4π·F c ·x c ) −γ = Ax −γ , where A ≈ 11229416.5 and γ = 2.0. For collision avoidance, the UAV RX antenna needs to provide the horizontal field of view (FoV) of 90 degrees, with the vertical FOV of 10 degrees (to minimize the impact of buildings). Maximum isolation is required otherwise, thus offering the gain of G RX ( − → d ) toward the UAVs in front of the RX, while attenuating all other possible targets by at least 20 dB.
However, in practice, for the carrier frequency of interest such antenna may not realistically be expected to yield a significant gain (or high front/back ratio), and thus may be approximated with an isotropic antenna for the first-order analysis purposes. By contrast, for the numerical assessment, a realistic antenna designed under the above constraints is employed. Another concern is the dependence of the TX antenna gain on the UAV altitude (due to cellular antenna downtilt), but that turns out not to be an issue for low-flying UAVs, as illustrated in Fig. 3. Specifically, at larger distances R, the BS antenna does not have sufficient angular "resolution" to significantly attenuate the signals toward UAVs flying near the ground. The illuminating signal is reflected by all the UAVs according to the 3D reflection function follows for a given direction of incident and reflected rays.
The overall power available at the RX for a link is determined according to the radar equation where r 1 and r 2 are the distances from the target to the RX and from the target to the TX, respectively, and A ef f = λ 2 4π is the effective antenna area. Its size depends on the relative BS locations, as well as the form and structure of the objects themselves as discussed in Section III.

C. Bistatic RADAR Operation
The most conventional RADAR system is monostatic, wherein TX and RX are colocated in one device, thus making them conceptually simpler. In hardware terms, however, most monostatic RADARs require a specialized transciever and RF components, and therefore are not suitable candidates for integration into concentional communications RX. In contrast, our distributed RADAR system is a superposition of multiple bistatic RADAR configurations, which require at least two nodes to operate, one transmitter and one receiver, and is thus more complex. On the other hand, bistatic configurations can be implemented with current communication hardware, as no individual node needs to act as TX and RX at the same time.
Due to bistatic operation of our system, the RX receives an original pulse from the BS and sees reflections from all the surrounding objects. Hence, if UAV A aims to specifically detect other UAVs, it needs to contend with: (i) direct leakage signal from the BS to UAV, which can in principle be removed by signal processing, (ii) desired UAV reflections, which in our case are the objects we actually need to avoid collisions with, (iii) reflections from buildings, which are treated as clutter.
To perform separation between buildings and UAVs, there are two components of information available for each potential target: (i) Doppler shift and (ii) power level. While the power level is straightforwardly computed using (3), the Doppler shift for bistatic RADAR is obtained as where λ is the wavelength, R T X is the distance between the TX and the target, and R RX is the distance between the RX and the target.

D. Metrics of Interest
Our parameters of interest follow from the safety concerns, in particular, UAV collision probability. Specifically, we seek the probability P d that the target is detectable by the RX, given that the distance between the RX and the target is less than a certain value R min = 20 m. This probability is affected by two key factors: (i) the probability that the received signal power is sufficient P RX > S, where S = −120 dBm is the RADAR sensitivity level, and (ii) the probability that the RX can separate the target from the background, i.e., that it has a sufficiently different Doppler shift. These two can then be combined to determine the successful detection probability P d .
Based on the successful detection probability, one can then infer the likelihood of the situation that two UAVs remain unaware of each other while approaching dangerously close, thus creating a collision risk. From the system analysis point of view, we are further interested in the following considerations: 1) Understand the impact of the quality of a UAV RCS model on the predicted system reliability. 2) Quantify the probability of a false detection, i.e., the chances that a certain UAV is reported as a potential collision threat, whereas in reality it is not (e.g., it is farther away than estimated, or it is not a UAV at all).

III. CALCULATION OF BISTATIC RCS FOR UAV
In this section, we present a novel method for calculating the scattered electromagnetic field under the assumptions of physical optics (PO). Then, we demonstrate that our approach converges to selected closed-form solutions, while offering superior speed w.r.t. GO. Our method is thus shown to be most helpful when the existing solutions cannot be applied due to their slow speed or poor accuracy.

A. Related Studies
There are three main approaches to obtaining the RCS of an object. First, it can be produced using direct measurements. This method requires an anechoic chamber and fairly expensive equipment [19], [20], but remains in use due to its accuracy.
The second approach relates to "full-wave" simulation tools that include Method of Moments (MoM), Finite Element Method (FEM), and Finite-Difference Time-Domain (FDTD). For MoM and FEM, the idea is to apply numerical analysis to Maxwell's equations in the form of ordinary differential equations (ODEs). This, in turn, requires a subdivision of the environment surfaces into smaller elements (cylindrical, rectangular, or triangular) [21], [22]. While being reasonably accurate, these methods face restrictions on scene complexity. The number of the corresponding ODEs needed to represent complex or non-homogeneous scenes increases significantly, thus causing growth in computational and memory burden [23]. This makes MoM and FEM solvers challenging to use where the scene is large relative to the wavelength.
For the FDTD method, the algorithms are based on a discretetime solution to Maxwell's equations, and thus call for a voxel (aka "Yee lattice") scene representation [24]. Similarly to MoM and FEM, a typical performance bottleneck of FDTD is in the voxel size, which must be made much smaller than the wavelength used in the modeling [23]. Therefore, similar to MoM/FEM, FDTD is best suited for smaller scenes (or longer wavelengths). A key feature of all full-wave methods is that once the solution is found for a given driving wave, the entire reflected field is known. Hence, if the TX location is fixed, their computational complexity may be acceptable.
The third alternative approach is analytical. Closed-form expressions for the RCS of simple targets, such as a sphere or a corner reflector, have long been available [25]. However, when the shape of an object is more complex, it becomes nearly impossible to obtain a closed-form solution. Because of this, much effort is devoted to studying a set of approximations named physical optics (PO). PO takes into account wave effects, such as interference, diffraction, and polarization, but does not consider the impact of induced fields. The rationale behind the PO approach is that the modeling outcome is close to that of the full-wave solutions [26], but is much less computationally expensive.
The PO approximation allows calculating the transmitted or scattered field [27] by solving the PO integral over the surfaces illuminated by the incident wave. The latter is the most computationally intensive part of the PO procedure. There are different techniques for the calculation of the said integral. Some authors propose an analytical solution for the PO integral, such as transformation of the surface integral to line integrals under certain constraints [28]. In other works, authors consider analytical expressions to calculate the PO scattered field with plane triangular patches [29]. Others put forward discretization of a surface as quadratic patches in such a way that the amplitude and phase terms of the PO integral change in quadratic form [30]. In [31], the authors offer a way for converting these quadratic patches to square patches.

B. Description of Proposed Method
The key idea behind the PO is that one only needs to model possible transformations of a wavefront at the interfaces between different media, which are represented by integral over all surfaces. We now proceed by describing our numerical method for its evaluation.
According to the proposed method, we replace the integration over a surface with the sum of the fields of tiny segments of the surface. For this, we associate each segment with its own RCS. For simplicity, it can be represented by a square shape. The RCS of a square plate is expressed as [26] where a is the side of the square segment that the point represents and λ is the wavelength in the current environment.
To accurately capture the geometry of a region on an object's surface, each segment carries a normal vector of that fragment of the surface. To take into account the normal direction, we introduce a directivity coefficient where θ ref is the angle between the point's normal vector and the direction to the RX. In practice, any macroscopic object can be readily represented as a uniformly distributed point cloud, wherein every point is the center of one segment as exemplified in Fig. 4. The power carried by a reflected spherical wave at the RX location using the radar equation for each segment [32] can be estimated as follows where P T X denotes the power of radiation at the TX, Dt i is the distance between the TX and i-th point, Dr i is the distance between the RX and i-th point, and A e = λ 2 4π is the effective antenna area at the RX (assuming ideal isotropic RX).
Finally, the overall electric field at the RX point can be derived as a coherent sum of the fields from each point as where φ i = 2π d λ is the phase of i-th component at the RX. Hence, surface subdivision followed by summing over multiple elements allows us to approximate the final electric field at the RX location, thus transforming the evaluation of the PO surface integral into a sum over multiple points, which is well-suited for vector processing, such as AVX and GPU. Our following objective is to verify that the results obtained with this approximate method remain sufficiently accurate for RCS calculations.

IV. VALIDATION OF PROPOSED METHOD
In this section, we report the validation scenarios, the selected modeling conditions, and the comparison of our method with the theoretical expressions. In most examples, we are interested in validating the results for the RCS, as those are key for RADAR processing.

A. Validation Scenarios and Conditions
We compare the performance of our proposed method with the analytical expressions for which closed-form RCS solutions are known [26]. There are certain constraints imposed by the closed-form solutions: (i) the incident wave has to be a plane wave, (ii) the calculations are performed in the far-field and for this we ensure the separation of 200λ, (iii) the sizes of objects should be more than ≈ 2λ, (iv) all objects are perfect electric conductors.
We consider the following validation scenarios: flat plane of different size, sphere, cylinder, and concave object. One can estimate the RCS plane for the bistatic RCS of a flat plane [26] via the following where c and d are the sides of the plane, θ inc is the angle of incidence (the angle between the surface normal and the vector to the TX point), and β = 2·π λ is the wavenumber. The monostatic RCS of a metal sphere and a cylinder, respectively, can be assessed as [26] RCS sphere = πr 2 s , where r s is the radius of the sphere, and where r c is the radius of the cylinder and h c is its height.

B. Validation Results
First, we assess the reflected power as a function of the angle of reflection for a square flat plate, since the theoretical expression for the latter is known (9) [26]. We fix the grid step at 4 mm, and the carrier frequency is 300 MHz. As can be observed in Fig. 5, for the plate size of 5λ, our method demonstrates adequate agreement with the exact solution. When the plate size is set to 20λ (see Fig. 6), the obtained result for the proposed approach remains similar to the exact solution. Therefore, for a range of sizes, we demonstrate that the proposed method captures both the main reflection lobe as well as all the side lobes correctly.
Second, we compare the results for the spheres of different radii [26]. We set r s = 3, 4, 5 m, and consider different distances between the sphere and the TX/RX point. For a sphere, only monostatic results are known in closed-form; hence, TX and RX have to be collocated. Here, our approach continues to demonstrate excellent agreement with theory, as seen in Fig. 7. It is important to note that it is extremely difficult to obtain similarly accurate results with a geometric optics (GO) solver   (such as a ray tracer) due to smooth convex shape of the object in question.
Further, we consider a cylinder. Its height is fixed at 8 m, while the radius is being varied (r c = 3, 4, 5 m). The dependence of the reflected power on the distance from the TX/RX to the object surface is displayed in Fig. 8. The curves demonstrate a perfect match. The final validation is done with a concave object,   specifically, a parabolic reflector with the radius of 5 m. As we can see in Fig. 9, the maximum of power is at the focus point for various shapes of the reflector, which matches the initial expectations. In Table III, we present the theoretical values of the gain for a parabolic reflector at the focus point together with our numerical results. These match the analytical expression. In summary, for a wide variety of shapes, our method produces results that are very close to those obtained via exact solutions. It has to be noted, however, that while our approach is far superior to GO (i.e., ray tracing and its variants), it does have its own limitations [33]. Specifically, for those cases where the objects are smaller than the wavelength, one can expect induced currents to affect the applicability of PO in general, thus making the proposed approach unreliable.
Similarly, for extremely high frequencies, the density of the point cloud needed for accurate calculations becomes rather high, thus straining the memory budget. Despite these shortcomings, the proposed method is observed to be superior to the MoM in terms of scalability, by approaching that of typical GO solutions, while providing accuracy far superior to what is possible with the GO, especially for smooth, organic shapes common for UAVs.

C. Computational Complexity
Here, we consider the computational complexity of the proposed method. First, one needs to assess how many points per object are required. As an example, the well-known NEC++ [34] engine may represent a sphere of radius 5λ with about N = 1000 elements using the MoM, where each element is internally an ODE. However, to solve the system of N equations, one needs at least O(N 2 ) operations, and a substantial amount of memory to store the intermediate results.
In our case, a sphere of the same size has N = 20, 000 points, but the solution speed is O(N ), with no extra memory requirements. In addition, with the proposed approach, one does not need to represent free space in the scene, unlike in FDTD; therefore, this reduces the computational cost as well. Further, in our method, while more points are required to represent an object, the volume of computations per point is reduced, and the code has no branching. Overall, there is a difference in the runtime speeds on the order of 1000 for the RCS of a sphere as compared to e.g., NEC++.

V. STOCHASTIC GEOMETRY ANALYSIS
In this section, we continue by outlining a stochastic geometry framework tailored to the performance assessment of our UAVbased communications-and-RADAR system.

A. Analysis Principles
The proposed analysis is based on stochastic geometry considerations. The overall framework comprises of thee stages. In the first one, we determine the probability that the target is actually distinguishable from the clutter. To this aim, we assess whether the difference between the Doppler shift from the target at the RX is sufficiently distinguishable as compared to the Doppler shift from all the scatterers located in the environment.
Distinguishing the target from the scatterers alone is not sufficient for a successful target detection. Hence, at the second stage, we determine the probability that the received signal is actually higher than the detection threshold. To this end, we account for the reflection profile of the target to characterize the probability density function (PDF) of the received signal strength and then obtain the sought probability of signal detection by numerical integration.
Finally, at the last stage, we deliver the overall target detection probability. Accordingly, we exploit the independence property for the events that (i) the received signal is greater than the threshold and (ii) the target is distinguishable from the clutter. In the end, the target probability of signal detection is obtained as a product of the probabilities of these events.

B. Target Visibility Probability
With the metrics of interest specified in Section II, we consider the top view of the scenario geometry illustrated in Fig. 10. Here, we assume that the RX positions are uniformly distributed across the coverage area of a cell with radius R. The targets of interest are uniformly distributed within a circle around the RX having radius r < R. To determine the target detection probability, one needs to ensure that: (i) the target is distinguishable from the clutter and (ii) the received power is sufficient to perform a detection. For the former to hold, the Doppler shift from the target at the RX, D RX , has to be sufficiently dissimilar from all the Doppler shifts D Sc from the scatterers located in the environment, i.e., buildings.

1) Doppler Spread From Stationary Scatterer:
The Doppler spread from a scatterer can be written as where f c is the carrier frequency, v T X,Sc rel is the relative speed of the TX and the scatterer, v Sc,RX rel is the relative speed of the scatterer and the RX. Since both TX and scatterers are assumed to be stationary, we may simplify (12) as Observing the structure of (13), one may notice that D Sc is the function of a single variable v Sc,RX rel . Due to our assumption on random UAV trajectories, the latter is a random variable (RV), even for a constant speed of the UAVs, v. By utilizing Fig. 11(a), we determine this component as where α is uniformly distributed in (0, 2π) and (13) reads as Note that (15) implies that D Sc is the function of a single RV. Recall that the PDF of RV Y , w(y), expressed as a function y = φ(x) of another RV X with the PDF f (x) is given by [35] is i-th branch of the inverse. Observe that applying the transformation directly to (15) is difficult. However, one can apply it initially to the component U = v sin α and then to the overall function. Under this transformation, each value of U located in the interval (−v, v) corresponds to multiple values of α, i.e., By employing the transformation in question, we arrive at As α is distributed uniformly in (0, 2π), one may notice that there are only two non-zero components in (18): (i) k = 0, x < v and (ii) k = 1, 0 < x < v. Hence, the final PDF of U = v sin α takes the following form With the distribution of U = v sin α at hand, we are now able to determine the PDF of D Sc . Since the transformation at hand is a bijection, and observing that the derivative of the inverse function is given by 2 , we arrive at the following PDF 2) Doppler Spread From Actual Targets: We now determine the distribution of the Doppler spread from any of the actual targets. Accordingly, we have Contrarily to the Doppler spread from a scatterer, here both components v T X,T g rel and v T g,RX rel are RVs. The former one can be determined similarly to (14), as v T X,T g rel = v sin α, where α is uniformly distributed in (0, 2π).
To obtain the component v T g,RX rel , we need to account for the relative speed between the target and the RX, which are both moving. By considering the projection of the RX trajectory onto the line connecting the RX and the target, see Fig. 11(b), we obtain the speed of the RX with respect to the target as where β is uniformly distributed in (0, 2π). The resultant relative speed is thus the vector sum of v and u T g,RX rel . To operate with scalar variables, one needs to supplement u T g,RX rel with an appropriate sign. Furthermore, since the trajectories of the target and the RX are assumed to be independent, α and β are mutually independent. With these observations in mind, the Doppler spread from the target can be written as thus implying that it is a function of two independent RVs. By following the same technique as the one we utilized to establish the Doppler spread from a scatterer, we first determine the PDFs of two independent RVs, U = v sin α and Z = v sin β, where both α and β are uniformly distributed in (0, 2π). We thus have Note that by applying the RV transformation technique, no closed-form expression for w D T g (x) is feasible. However, one can readily determine the cumulative distribution function (CDF) of D T g by numerically estimating the following integral (23), while w U (x) and w Z (y) are available in (24).
As an alternative solution, one may first perform linear transformations of the numerator and the denominator as U = 1 + U and 1 − v). Then, to estimate the ratio U /V , one needs to numerically solve the following integral: By differentiating (26), one can obtain which is easier to estimate than (25). To complete the derivation, one needs to perform a linear transformation of U /V according

3) Probability of Target Visibility:
After determining the PDFs of D Sc and D T g , we aim to establish D T g ∈ {D Sc − Δ D , D Sc + Δ D }, where Δ D is the Doppler threshold, which can be straightforwardly derived from the RX UAV groundspeed. Mathematically, we seek the probability By utilizing (28), the task of quantifying p V thus reduces to determining the CDF of the modulo of the difference between D T g and D Sc . The PDF of the difference D − is readily available by performing a convolution of the RVs D T g and −D Sc , i.e., The PDF of the absolute value |D dif | is made available by utilizing a transformation of the RV D dif with respect to the modulo function. In this case, the inverse function has two branches φ 1 (y) = −y and φ 2 (y) = y, thus leading to the sought PDF in the following form To obtain the CDF of |D dif f |, one needs to integrate (30), and then the sought probability in (28) readily follows. Recall our assumption that the scatterers follow a homogeneous PPP in 2 . Hence, by employing the properties of the PPP, the target visibility probability in a field of scatterers with density λ C is given by where the first component provides the probability that there are no scatterers in the region covered by the BS antenna, A, (i.e., the void probability of the PPP), while the second component specifies the probabilities that there are k scatterers in the region A but none of those occlude the visibility of the target according to (28). The only unknown in (31) is the mean number of scatterers in the region covered by the BS antenna, Λ. This parameter is related to the density of scatterers as S A λ C , where S A is the area of the region covered by the BS antenna. In general, S A depends on the type of the BS antenna, the operating frequency, and the frequency reuse plan. For example, for a tri-sector antenna, the area S A can be approximated by a triangle with 120 • angle. If the carriers utilized in all the sectors are the same, then the region of interest can be well approximated by a circle.

C. Power Detection Probability
Observe that the visibility of the target does not guarantee that the latter is detected successfully. Even when the target is visible in a field of scatterers, the received power may be insufficient to successfully detect it. Hence, here we derive the probability that the received power at the RX is higher than a certain threshold Δ P , i.e., q P = P r{P RX ≥ Δ P }.
The received power at the RX can be written as where P T X is the BS transmit power, G T X (ζ, θ) is the BS transmit gain toward the current target location, G cs (ζ, θ, ζ , θ ) is the RCS, G RX (ζ , θ ) is the receive gain at the RX, γ is the propagation exponent, A is the frequency-dependent attenuation factor, L 1 and L 2 are the distances shown in Fig. 10, ζ, ζ and θ, θ are the elevation and azimuth angles, respectively, in the appropriate direction. These parameters are illustrated in Fig. 12.
By analyzing the structure of (32), we observe the following. The RCS is derived in Section III as a function of the azimuth and elevation angles. The parameters G T X (ζ, θ) and G RX (ζ , θ ) are available as part of the problem formulation. It is important to note that by assuming fixed BS and UAV heights, the elevation angle ζ equals ζ = arccos([h U − h A ]/L 1 ), while ζ = 0. Furthermore, both azimuth angles θ and θ are RVs with the densities w θ (x) = 1/2π and w θ (x) = 1/2π, respectively.
The only remaining unknown in (32) is the distance (L 1 + L 2 ). Let us denote the projections of these segments L 1 , L 2 , and L 3 onto the xOy axis by L 1,x , L 2,x , and L 3,x , respectively. We see that L 2 = L 2,x , while where L 2,x is the RV with the PDF w L 1,x = 2x/r 2 , 0 < x < r, L 3,x is the RV with the PDF w L 3,x = 2x/R 2 , 0 < x < R, while the distribution of L 1,x is unknown. We can also express ζ as By combining the above facts together, we establish that the received power at the RX in (32) is a function of four RVs, L 2,x , L 3,x , θ, and θ . However, since G cs (ζ, θ, ζ , θ ) is not available in the closed form, a direct application of the RV transformation technique does not lead to a closed-form solution. However, one can estimate the CDF of P RX via numerical integration as where f (x 1 , x 2 , x 3 , x 4 ) is the received signal strength. Accordingly, the probability that the received signal strength at the RX exceeds a certain threshold is given by Finally, if the mutual orientations of the mobile objects are known a-priori (e.g., the considered target and the RX are part of a fleet moving in the same direction), θ 1 and θ 2 become constants. Then, the received power depends only on two RVs, L 2,x and L 3,x , thus leading to a simpler form Note that our analysis can be extended to incorporate variable BS heights as follows. First, one needs to make the said height a random variable. Hence, the variables L 1 , L 3 , and ζ will become functions of more than one RV. Second, solutions for L 1 and L 2 in (33) can be obtained in a closed form for some suitable distributions of BS heights (i.e., those having a single exponential term, such as exponential, Rayleigh, and Weibull).

D. Successful Detection Probability
Having both the target visibility probability and the probability that the received signal strength at the RX is lower than a certain threshold, we are now in a position to determine the probability of successful target detection given that it is located closer than a certain distance r. Recalling that the two derived parameters are implicitly conditioned on the target being closer than r, and by combining (31) and (36), we arrive at where Δ D and Δ P are the Doppler and power thresholds.

VI. SIGNAL RECEPTION AND TARGET DETECTION
In Section V, we presented a RADAR detection probability analysis based on stochastic geometry. While useful for firstorder system-level assessment, that approach has limitations in capturing finer details. In fact, the detection of a target relies on the specific signal waveform, receiver design, and signal processing algorithms. In this section, we elaborate the signal reception and target detection procedure. It is then used in our Monte-Carlo simulations of target detection probability.

A. Matched Filter Reception
The first step in the RX processing is the matched filter, which maximizes the signal-to-noise ratio (SNR) in the presence of additive stochastic noise. The timing of the maximum value of the matched filter outcome characterizes the distance to the detected target, while its magnitude is used for signal detection. Typically, matched filtering is implemented by correlating the signal with the conjugation of a known waveform (or sequence, if digital processing is used) as where s(t) is the transmitted waveform, r(t) = r(t) h(t) is the received signal, and (·) * denotes the conjugation operation. From (39), the waveform s(t) determines the SNR gain of the matched filter. To illustrate the detection capability for the conventional communication signals, we employ the SRS defined in the 3GPP specifications [36], which is typically based on a Zadoff-Chu sequence and brings approximately a 20 dB SNR gain.

B. Range Doppler Profile of Received Signal
The outcome of the matched filter represents a power-range profile of the environment. In our case, however, since the UAV antenna is not perfect, the range alone may not be able to discriminate the targets from the clutter. Hence, we further explore the power-frequency (Doppler) domain profile jointly with the power-range profile to improve the discrimination capability of the detection system. The range-Doppler profile is obtained via cross-ambiguity processing as The output of (40) characterizes the power distribution of the received signal in frequency (Doppler) and time (range) domains as compared to the transmitted waveform. x(τ, f d ) is also known as the ambiguity surface. The time delay and Doppler shift carried by the reflections from the UAVs can be represented by the value of x(τ, f d ) at the specific time instant and frequency point, which is the peak (local maximum) value on the ambiguity surface x(τ, f d ).

C. Detection Rate With Constant False Alarm Rate (CFAR)
The local maximum values of x(τ, f d ) indicate the delay and Doppler shift of a reflected signal. In this work, we utilize a cell-averaging CFAR (CA-CFAR) detector to capture the local maximum (see Fig. 13). The CA-CFAR detector compares the intensity of the cell under test (CUT) with a dynamic threshold that is calculated via the following steps. First, the threshold factor is chosen according to the expected CFAR level P fa and is given by where N is the total number of the selected training cells. Further, the threshold is determined by the threshold factor α and the estimated noise level n est according to where n sum is the sum of all the training cell values. To avoid signal leakage into adjacent cells, this method exploits guard cells to isolate the CUT. Once the threshold is determined, the detection decision is made based on where p CUT is the value (power) of CUT. The CA-CFAR procedure is applied to the ambiguity surface x(τ, f d ) cell by cell for determining the detection status of every cell.
Regarding the computational complexity of our RADAR algorithm, which we use for determining the detection state of every cell on the ambiguity surface, the CFAR procedure averages the values of the selected cell and the training cells around it (see Fig. 13). Since the averaging process repeats for every cell that appeared on the ambiguity surface, the computational complexity of RADAR detection depends on the total size of the cells on the surface. In other words, the CFAR complexity is determined by the multiplication of frequency domain cell numbers (Num F ) and the range domain cell numbers (Num R ). This work employed 1024 and 2048 for Num F and Num R , respectively, and the complexity is thus O(Num F Num R ).

VII. SELECTED NUMERICAL RESULTS
In this section, we aim to verify the scattering model proposed in Section III with the RADAR detection metrics. By utilizing this approach, we can obtain the bistatic RCS of the UAV model in any direction, see Fig. 14.

A. Simulation Setup for Radar Detection
The specific 5G waveform used is the SRS, which features Zadoff-Chu sequences [36]. We consider a system with 2048 subcarriers having 15 kHz subcarrier spacing. The overall spectral span of the signal is 30 MHz. By analyzing a total of 400 SRS symbols, the range and Doppler resolution can reach 9.76 m and 30 Hz, respectively, which is adequate for UAV systems. To reduce sensing times, one may increase the SRS rate at the BSs. After propagation modeling, the range-Doppler profile is obtained via ambiguity processing. The noise floor of the RX is assumed to be −130 dBW. In the CA-CFAR, 5 guard cells and 25 training cells are employed on the range-Doppler surface.
The simulation environment follows the description in Section II with 19 BSs, 33 UAVs (one drone is chosen as the tagged RX), and 100 buildings. We collect all the parameters used by our simulation setup in Table IV. The locations of the UAVs and buildings are randomly generated in 5000 independent Monte-Carlo runs. Hence, the data set of the reflection channel model has 5000 independent instances and 19 × 132 = 2508 signal paths in every run. Fig. 15(a) presents an example range-Doppler profile of the signals reflected from the UAVs and buildings with marks. The strongest area in Fig. 15 corresponds to −110 dBW, while the weakest reflections are below −150 dBW. To select the local maximum values from the range-Doppler plot, CA-CFAR is adopted, which calculates the threshold factor α according to the pre-defined FAR level P fa as per (41). Then, the actual threshold level is determined according to (42).
The detection probability of CA-CFAR at 10 −5 FAR level is demonstrated in Fig. 15(b). The highlighted region embodies the reflections from targets. The FAR is utilized to control the impact of noise, where higher acceptable FAR directly influences the probability of detection. In our case, owing to the matched filter, a signal with the power of ≈ 20 dB below the noise floor of −130 dBW remains detectable; hence, the true sensitivity is −150 dBW.
When using a single nearest BS as the illuminator, we note around 20% UAV detection probability, even with the very high FAR of 10 −3 , primarily due to nulls in the RCS patterns. However, due to diversity in bistatic geometries, RX power from different BSs varies widely. In this situation, a target (UAV or building) that is "invisible" when using one BS may be easily detectable with another BS. Hence, multiple illuminators provide extra robustness for the detection purposes. The detection performance with multiple illuminators may be improved if the detection of the same target with different illuminators is independent, which is generally true, as: where P M is the detection probability by using M -th BS as an illuminator. With fewer BSs, the detection rate of targets is less than 25%. Taking advantage of multiple BSs, the detection rates are improved to over 80% even at low FAR levels.

B. Discrimination of UAVs and Buildings
The detection rate evaluated above does not differentiate between UAVs and buildings, neither is it straightforward with power alone. Further, from a single range-Doppler plot like Fig. 15, one may not discriminate the types of targets. In what follows, the profiles of different types of targets are investigated Fig. 16. Detection rate of one UAV if using multiple reflections from different BSs jointly. Carrier frequency is 800 MHz. and the strategies for target discrimination are discussed. Fig. 16 depicts the detection rate of one UAV if using multiple reflections from different BSs jointly.

1) Target Discrimination With Constant RCS:
For the fixed RCS parameters of UAVs and buildings (as we set the mean values of the actual RCS distributions), the reflected signal power from the UAVs and buildings is distinguishable as confirmed in Fig. 17(a). The mean values of the reflected power for all 19 BSs are collected in Fig. 17(b). It is apparent that the mean signal power from the UAVs is higher than that from the buildings most of the time; hence, one can discriminate based on the signal power. By setting the middle point of the UAV and building reflection power as a threshold, the "Hypothesis Testing" analysis is applied to the CFAR output for each BS. The "Hypothesis Testing" observes the power of all the detected signals and compares these values against the threshold, which is plotted as the yellow line in Fig. 17(b). If the power is above the threshold, we treat the detection outcome as a UAV, otherwise it is a building. The "Hypothesis Testing" results for all the single illuminator cases are demonstrated in Fig. 17(c). As the reflection power gap between the UAVs and the buildings is large for BS 1 to 10, we observe that the percentage of 'TP' remains higher than 60%. When building reflections become stronger and power gap reduces, as shown by the remaining cases in Fig. 17(b), incorrect discrimination occurrences like 'FN' and 'FP' in Fig. 17(c) emerge. To support stable discrimination, we extend to multiple illuminators at the discrimination stage.   With the discrimination based on each individual BS, majority voting is applied. The joint discrimination output is offered in Fig. 18(a), where we see that the correct discrimination rate improves from around 76% to over 90% via combining more signals from different BSs. In comparison, the misclassification rate reduces at the same time from 25% down to 10%.
2) Target Discrimination With Realistic RCS: Fig. 19(a) and (b) report the distribution of the reflection power from buildings and UAVs when using the signals from BS 1 and BS 3, and a realistic RCS. We notice in Fig. 19 that the power from a given type of target tends to reside within a certain range for any given BS. We extend the discrimination strategy proposed above to improve over this. We plot the average reflection power from different BSs in Fig. 19(c) together with the discrimination threshold. With the realistic RCS, one can see that the reflection power from the UAVs and buildings does not show distinguishable features as in the constant RCS cases.
Therefore, we jointly analyze the detection results with multiple BSs, and the output in Fig. 18(b) reaches 65% of the correct discrimination rate if we use 8 BSs as illuminators in total. It should be noted that the varying and weak reflection power in Fig. 19(c) limits the maximum correct discrimination rate to around 65% even if we leverage all 19 BSs. Moreover, this strategy may not be generalized because the power distribution varies as the bistatic geometry changes. Hence, to discriminate different types of targets with realistic RCS, one needs to exploit the time and frequency shift properties of the continuously detected targets from the (time) series of the range-Doppler plots.

VIII. CONCLUSION
In this article, we introduce a novel approach to the design and analysis of an integrated sensing and communication UAV system. This co-design methodology allows for conceptualizing an efficient converged solution that permits UAVs within a city to perform resilient RADAR collision avoidance without significant investments into RADAR hardware on each UAV, and with a possibility to reuse the ground network infrastructure. As a potential future research direction, our system model may be complemented by a detailed BS-to-BS interference model, which then requires a cellular network deployment model.
A key enabler to this work is our novel RCS calculation algorithm, which allows approaching the challenge of bistatic distributed RADAR modeling. This RCS characterization method demonstrates highly accurate results under moderate computational complexity. Our deterministic approach to channel modeling may be extended further by considering the multi-bounce effect for capturing more than one reflection component in multipath environments. For this, a highly efficient voxel cone tracing algorithm can be combined with our introduced point-cloud method.
Moreover, we demonstrate how analytical modeling can be applied to such a complex scenario, all while taking into account the details of RCS patterns and Doppler processing in the RADAR receivers. Finally, we illustrate the system performance and scaling with extensive numerical results, thus conclusively demonstrating the robustness of our method to RCS pattern nulls, which is crucial for safety-centric UAV applications such as collision avoidance.