Optimal Search Strategy of Low-Altitude Target for Airborne Phased Array Radar Using Digital Elevation Model

The search for low-altitude targets becomes a more challenging task for airborne phased array radar, due to the limited system resources and serious ground clutter. In this article, an optimal search strategy using the digital elevation model (DEM) is proposed to allow for a fast and high-probability searches for low-altitude targets. First, a clutter-free detection region is formed by introducing the DEM. Compared with the traditional methods, the detection probability of the low-altitude targets can be improved. Second, the search problem is transformed into a constrained multiobjective optimization problem. The proposed strategy can simultaneously minimize the search time and maximize the target detection probability. Finally, two intelligent evolutionary algorithms are compared to verify the proposed strategy, and the proposed strategy is also testified by utilizing real DEM data of Dujiangyan, Sichuan, China. Simulated results indicate that the proposed optimal search framework improves the detection probability and reduces time cost for low-altitude targets.


I. INTRODUCTION
I N MODERN electronic countermeasures, the advancement and integration of various information technologies have led to the rapid development of airborne phased array radar (APAR) [1], [2], [3], [4], [5], [6], [7]. Compared with traditional mechanical scanning radar systems, the electronically controlled array antenna used in APAR systems provides strong beam agility and beam control capabilities, which is significant for the primary task of APAR: the fast and accurate airspace search [8], [9], [10]. However, in the current APAR system, the conflict in the consumption of system resources between different radar functions severely limits its airspace search performance. That is, inadequate time resource leads to longer target discovery time and lower discovery probability. In addition, for the low-altitude target's searching and detection, the complex ground environment and clutter interference greatly increase the difficulty of target detection. Therefore, the poor detection probability with strong ground clutters needs to be handled. A traditional method of searching for a weak-scattering moving targets against ground clutter is sequential beam scanning with pulse Doppler (PD) technology [3]. PD technology relies on the relative motion between airborne platform and ground, which has been used to separate the target from ground clutter in the range-Doppler domain. However, this technology requires a high pulse repetition frequency to form a clutter-free region by avoiding the inherent Doppler folding problem [3], [11], [12]. Moreover, when the target's radial velocity is low, the target will fall in the sidelobe clutter region, and the proper ground clutter radar cross-section (RCS) fluctuation models and constant false alarm rate (CFAR) techniques should be studied and adopted [13], [14], [15]. The abovementioned methods do not combine the optimal search problem.
The search problems for APAR are usually divided into two categories: one is autonomous search and the other is cued search [3]. Many solutions have been proposed to improve the detection probability by improving the autonomous search strategy. Without any cued information, it is usually necessary to search the wide-range airspace or establish several search screens in different pitch ranges [16]. Relying on a long time cost, a repeated observation search strategy is proposed for space satellites [17]. This strategy provides a reliable detection probability, but it is difficult to extend to low-altitude targets due to the unacceptable observation time requirements. In [18], to minimum the radar search load, while ensuring the detection probability, a scalar optimization strategy is proposed that merges the multiple surveillance parameters. This method improves the computational efficiency of optimization algorithms. However, it considers only the one-off and cumulative detection probability without optimizing the target search time.
Fortunately, some cued information, such as the rough location, the velocity, and the motion trajectory prediction of the target, is provided by the space-based or ground-based early warning systems [3]. Using these cues, the beam dwell time of APAR can be optimally designed to achieve higher detection probability and lower system resources cost. In [19], a search strategy based on the principle of maximum information gain is proposed. The information gain of each beam position is This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ predicted before radar detection, and each beam position is illuminated in sequence according to the sorted information gain. In [20], by considering the target appearance distribution and platform velocity cued information, a resource allocation method is proposed for both airborne and ground-based radars.
Relying on the scenario prior information, the method can enlarge the search surveillance airspace under the same detection probability. In [21], a cued search parameter optimization model is derived under the accumulated detection probability criterion to obtain a theoretical optimal search data rate. In [22], through the introduction of the target threat level prior information, different search airspaces are weighted by manually defined airspace threat vectors. However, the method strongly relies on manual weights, and improper weights may lead to longer target search time. In fact, the cued search strategies can generally provide superior performance over that the autonomous search strategies.
In this article, the digital elevation model (DEM) data are introduced in APAR to improve the detection probability of low-altitude targets. DEM is a 3D computer graphics representation of elevation data to represent terrain. It is widely used in surveying and mapping, hydrology, meteorology, and geomorphology fields [23], [24], [25]. The DEM dataset is quite mature and can be easily acquired by global navigation system and synthetic aperture radar interferometry [26], [27]. To avoid complex clutter RCS fluctuation modeling, in [28], a ground return suppression method using a terrain database has been proposed for weather radar. This system takes advantage of the DEM to provide the accurate elevation information for the low-altitude clouds remote sensing [29], [30]. However, this method is only suitable for cloud imaging, and cannot be used for optimal target search.
Inspired by the abovementioned papers, this article proposes an optimized search strategy relying on the DEM to optimize the time resources and increase the detection probability of APAR. First, by exploiting the geometric relationship between the DEM and the platform position, the optimal range scope of each beam position is calculated. Clutter-free airspace is formed to improve the detection probability of low-altitude targets. Compared with the traditional target detection method in the clutter, the noisedependent detection probability is higher. In addition, there is no CFAR detection loss in the clutter-free region. Second, based on the detection result, the search problem is established as a constrained multiobjective optimization problem. Compared with the traditional method that only considers the minimum search time, this method additionally considers the maximum target detection probability. Multiobjective particle swarm optimization (MOPSO) [31] is used to solve the objective functions to optimize the system resources. Finally, utilizing open-source DEM data of Dujiangyan City, Sichuan, China, the effectiveness of the proposed DEM-aided optimized search strategy is verified by comparing it with the traditional average weighted sequential search method and a method based on the nondominated sorting genetic algorithm II (NSGA-II) [32]. In summary, derived from the proposed DEM-aided optimized search strategy, the maximum average target detection probability is improved, and the minimum average target search time is reduced.
The rest of this article are organized as follows. Section II shows the search geometric models of low-altitude target search in APAR. In Section III, the proposed optimal search strategy is detailed. Section IV gives the simulation results used to verify the performance of the proposed method. Section V discusses the influence of the DEM error and the limitation of our method. Finally, Section VI presents the conclusions.

II. PROBLEM FORMULATION
The geometric model of the APAR platform and search airspace is shown in Fig. 1. In the Oxyz Cartesian coordinates, the APAR platform moves toward the y-axis with a constant velocity v. The initial position is (0, 0, H), where H represents the platform height. In the radar spherical coordinates, the azimuth 0 • direction is the x-axis, and counterclockwise rotation of the angle denotes the positive direction. The pitch 0 • direction is the y-axis, with the positive angle oriented toward the ground, as seen in Fig. 1. The azimuth and pitch ranges of the search airspace are [θ min , θ max ] and [ϕ min , ϕ max ], respectively. Following the triangular beam packing lattice [33], the search airspace is divided into N Ω cells. Each cell, Ω i , i = 1, 2, . . ., N Ω , represents the radar main lobe beamwidth coverage. [R min , R max ] means the detected range scope.
The foremost purpose of this article is to solve the low-altitude target search problem for APAR. First, within the detection range scope, the detection probability of the target is influenced by ground clutters generated by the complex environment. Actually, only the aerial targets are of interest. Thus, the detected range scope of different beam positions should be truncated properly according to the distance between the APAR platform and the true ground elevation. Second, the traditional uniform sequential search strategy cannot offer an optimized scheme for time resources, resulting in the illumination redundancy of each beam position. Hence, how to reasonably allocate the illumination time of each beam position to achieve the optimal search are essential point of this article.

III. PROPOSED OPTIMAL SEARCH STRATEGY
In this section, the DEM is first introduced to reject the ground clutter. A ground clutter-free region is calculated to improve the detection probability. Then, the cued information of the rough target location is used to reconstruct the prior distribution weight of the target in the search airspace. Finally, a constrained multiobjective optimization problem and its solution are given to optimize the illuminations for the each beam.
A. DEM-Aided Ground Clutter Rejection 1) Principle: In [34] and [28], a terrain dataset is employed to reject ground clutter for airborne weather radar. Inspired by these papers, in this section, the DEM is introduced to form a clutter-free region for the APAR target search. Fig. 2 portrays the radial projection of radar sweeps in the cross-section. As seen in Fig. 2(a), the radar receiver suffers the interference from strong ground clutter, which will dramatically deteriorate the search difficulty. The transmission along radial  ) three produces at least two distinct echoes, as the transmitted beam will illuminate at least two distinct terrain surfaces [28]. As neither of these reflections represents aerial targets, the information is not merely unnecessary but generates confusion as a result of its presence. Because the aerial target search is the principal purpose of the APAR, the elimination of both of the reflections will not degrade the effectiveness of the APAR. Using the DEM and the position of the APAR platform, for each search beam, an optimized range scope can be calculated. According to this range scope, a high detection probability can be obtained within the clutter-free region, as shown at the bottom of Fig. 2(c).
2) Implementation: From Fig. 1, the position of the radar platform at time t is expressed as (0, vt, H). Suppose the beam direction of the ith beam position is (θ i , ϕ i ), and the azimuth and pitch half-power beamwidth are, respectively, Δθ × Δϕ. Then, the available beam coverage footprint can be written as During the flight, the DEM information of the real-time search airspace is extracted part by part to calculate the clutter-free region. This prior information is stored in the system before the radar platform takes off. Suppose that η(x, y) is extracted at t. Considering a large pitch angle, the beam at ϕ i may be affected by the terrain of the higher mountains at ϕ < ϕ i ; thus, the position of the closest point from the radar platform to η(x, y) To find R i max , the ground elevation η(x, y) in x-y coordinates is first mapped to the θ-r polar coordinates as η polar (θ, r) with an origin point of (0, vt). The mapping relationship can be expressed as where and denote the row ordinal number and column ordinal number of the respectively. x interp ∈ R M , and y interp ∈ R N are the range vector, and azimuth vector in the x-y coordinates after the nearest interpolation. m = 1, 2, . . ., M and n = 1, 2, . . ., N are the x interp and y interp coordinates when (2) and (3) are minimized, respectively. M and N are the interpolated range and azimuth dimensions. θ n θ and r n R denote the azimuth angle and ground range after mapping transformation, respectively. For the ∀θ k ∈ (θ i − Δθ, θ i + Δθ) cross-section, the equation of the radial with the pitch angle ϕ l on the plane ROZ is expressed as where ϕ l ∈ (ϕ i − Δϕ, ϕ i + Δϕ), r = [r 1 , r 2 , . . ., r n R , . . ., r N R ] ∈ R N R denotes the ground range vector after the mapping, and z = [z 1 , z 2 , . . ., z n R , . . ., z N R ] ∈ R N R denotes the height vector corresponding to r on the radial (4). Then, the ground range from the focal point of the radial z and the ground elevation η polar to the platform can be expressed as where the row order of η polar (θ k , r n R ) is calculated as argmin n θ |θ − θ k |. For beam direction (θ k , ϕ l ), the distance from the radar platform to the ground can be expressed as where the column order of η polar (θ k ,r) is similar to (5), which is calculated by argmin n R |r −r| is regarded as the collection of the ground range from the radar platform to η polar (θ, r) in the ith beam footprint. k = 1, 2, . . ., K i and l = 1, 2, . . ., L i represent the azimuth and pitch variables within a beam footprint. K i and L i denote the azimuth and pitch dimensions of this footprint. Finally, we can obtain the optimal range gate for the ith beam position as where min θ and min, ϕ respectively, represent the minimum element of the row and column of R i max . This means that R i max is optimized as the longest range within the footprint of the ith beam position without touching any ground reflector. That is, at the current time t, for the beam position (θ i , ϕ i ), within the range gate [R min , R i max ], the detection of a low-altitude target will not be interfered with by ground clutter.

B. Target Detection
The search radar equation is first given in this section. According to the obtained signal-to-noise ratio, the disadvantages of traditional target detection methods with clutter are analyzed. Then, the target detection method and detection probability of the clutter-free region are refined.

1) Search Radar Equation:
To search the target efficiently, the APAR illumination resources are reasonably allocated for different beam positions Ω i . T d denotes the dwell time of a single illumination. ΔΩ represents the solid angle of the antenna beamwidth. According to the search radar equation, the signalto-noise ratio of the ith beam position is expressed as [10], [35] where P av denotes the average transmit power, A e = Gλ 2 /4π is the effective aperture area, G is the antenna gain, σ is the RCS, k B is Boltzmann's constant, and T e is the radar system temperature in kelvins. L s denotes the total search loss, which includes the atmosphere attenuation loss, beam shape loss, and integration loss [35]. It is shown that for a given search range and solid angle, the signal-to-noise ratio increases with the beam dwell time.
2) Target Detection Without R i max : Since the distribution of average clutter power in each range segment is diverse, to ensure a CFAR, it is necessary to estimate the average power of the clutter and set the threshold according to this power. In this case, the cell averaging CFAR (CA-CFAR) detection method is commonly used [36]. Assuming the noise of the radar's I/Q channel obeys the Gaussian distribution, using a square-law detector, Swerling case 1 target, the detection performance of CA-CFAR is expressed as [36] where SNR denotes the mean signal to noise ratio of N samples, is the threshold product factor, andP f denotes the mean false alarm rate. When there is a strong interfering target in the reference range bins, the detection threshold of CA-CFAR is increased, which makes it difficult to detect weak targets. In this case, it is suitable to use the smallest-of CFAR method [37]. But when detecting low-altitude targets, the range distribution of clutter is seriously uneven. The detection performance will decrease sharply, and the greatest-of CFAR detector should be used in this case [14]. However, the abovementioned detection methods suffer from varying CFAR losses, which means that they require higher SNRs [36]. To guarantee the desiredP f , a higher threshold is necessary to compensate for the imperfectly known interference power.
3) Target Detection Using R i max : According to the previous section, we use the height information of the DEM to obtain the clutter-free detection range of each beam position. During the process of target detection, the square-law detector is used for the clutter-free range. The detection probability for N pulses under the Swerling case I model is written by [38] where T is the detection threshold, and I(·) denotes the Pearson's form of the incomplete gamma function [36]. When N = 1, the false alarm rate satisfies P f = e −T . In this case, the detection probability of a single pulse for the ith beam position is written as [20], [36], [39], [40] where SNR 0 = P av A e σ/4πk B T e L s denotes a constant related to the target RCS and radar system parameters.

C. Target Search
The target search algorithm can be divided into two parts. One is the transformation of the relationship between cued information and target weight distribution. The other is the establishment and solution of the search optimization function.

1) Target Weight Distribution Calculation:
We assume there is a cued location at (x 1 , y 1 , z 1 ) in the search airspace. According to the geometric relationship between the radar platform and the cued location, the real-time (time t) beam position corresponding to the airport can be expressed as where θ 1 and ϕ 1 represent the corresponding azimuth and elevation beam positions, respectively. Since the precision of the cued location is limited, referring to [19], the true target position can be expressed as where (δ θ , δ ϕ ) is the error in the azimuth and the pitch direction, respectively. It is assumed that the error satisfies a two-dimensional Gaussian distribution model with a mean value of (μ θ , μ ϕ ) and variance of (σ 2 θ , σ 2 ϕ ). Therefore, the probability density function of the position of the target in the azimuth and pitch directions can be expressed as f 1 (θ 0 , ϕ 0 ).
If there are multicues in the search airspace at the same time, the weights should be assigned by their prior information, such as target type, distance, and velocity. The weight of the mth cued location is written as [22] , q = 1, 2, . . ., M (15) where S type represents a [0,1] threat level defined by the target type from the actual demand and the expert experience. M denotes the number of the cued location within the search range at t. S r denotes a [0,1] threat level related to the target range where [r min , r max ] denotes a critical distance. S r denotes a [0,1] threat level related to target velocity where [v min , v max ] denotes a critical velocity. According to (15)- (17), the target appearance weights of the ith beam is written as where f m (θ 0 , ϕ 0 ) is the probability density function of the target appearance distribution of the mth cued location. Due to the beam broadening effect of APAR, the calculation of the integral in (18) is quite complicated. Considering the invariable characteristics of the beam shape of the sine space coordinates, the beam position arrangement and solution of (18) are usually calculated in sine space coordinates. The target position probability density function f m (θ 0 , ϕ 0 ) is converted to sine space coordinates as [41] Then, (18) can be rewritten as To simplify the calculation, it is assumed that the weights of target appearance at any position within the same beam coverage are equal to the weights of the center of the beam. That is, ignoring the difference within the beam, w i can be further written as where g m (α i , β i ) denotes the target appearance probability density at the center of the ith beam position for the mth cued location. D represents the beamwidth in the sine space coordinates, which satisfies D = sin( 0 ). 0 is the PAR beamwidth in the radar spherical coordinates. In this section, the prior weights of the target appearance for each beam position from the cued location and the threat level. It is seen that the prior weights are related to the beam position (α, β) and the weight ρ m of the cued location. In this way, for each beam Ω i , we can optimally allocate the illuminations of each beam position according to w i . For beams with high target appearance weights, the number of illuminations is increased, and for beams with low target appearance weights, the number of illuminations is decreased. Thus, the limited time resources can be maximized for the search procedure by a certain optimization.
2) Construction and Solution of Multiobjective Optimization Problem for Searching: Wu et al. [22], simply considered the minimum target discovery time consumption without considering the maximum target detection probability. In this section, both target detection probability and target discovery time are taken into consideration based on w i . Within a period, the ith beam position is illuminated k i times. The probability of the target being detected with a certain RCS for each illumination is calculated by (2) as P d i ,j . Then, the target detection probability for the ith beam position during this period is expressed as [33] Suppose that in a short period, the probability of each detection for the ith beam position remains a constant P d i ; then, (22) can be simplified as which satisfies P d i > 0, lim k i →∞ P i = 1. Using the prior weights from the cued information, the average target detection probability of all beam positions for the entire airspace is written asP where w i represents the probability distribution of the target's appearance at each beam position, which was calculated in the previous section. If a target is located at the ith beam position, the probability that it is found in the nth search frame can be denoted as Let the search frame of the ith beam position be t f i and the target appearance obeys a uniform distribution. Then, the time consumed to find the target in the nth search frame can be written as [22] T n = (n − 0.5) t f i .
Combine (25) and (26), the average target discovery time consumption of the ith beam position is expressed as The search frame for each beam position can be written as where T g is the total time resources of the APAR. Substituting (25), (26), and (28) into (27), the average target search time consumption is For (29), Wu et al. [22] uses the Lagrangian multiplier method to find the analytical solution of the objective function. However, this analytical solution is solved based on the search in a small indicated airspace, and P d i is assumed to be a constant.
In this article, due to the large search airspace, the changes in P d i at different beam positions cannot be ignored. In addition, for the target optimal search, we always hope to find the target as soon as possible while increasing its detection probability as much as possible, that is, the maximum detection probability and the minimum search time consumption should be satisfied at the same time. Consequently, under limited search time resources, a constrained multiobjective optimization problem that maximizes the average target detection probability and minimizes the average target search time consumption is constructed as where (32a) denotes the limitation of search time resources and (32b) expresses the boundary of k i for the ith beam position. L s is the ratio of the search time to the APAR's total time resources. In this case, it is difficult to solve the analytical expression of the optimal illuminations k i,opt through some mathematical approaches [22]. Solutions for multiobjective optimization problems are usually divided into two categories: 1) traditional optimization algorithms and 2) intelligent optimization algorithms [42]. After establishing the optimization model of the search problem, we analyze the type of the optimization problem, which belongs to the dual-objective optimization problem. Considering that methods for solving this optimization problem have been extensively studied, stable and accurate algorithms have been proposed. The solution algorithms adopted in this article are NSGA-II and MOPSO [31], [32], which are widely used in applications [43], [44].
In [31], a new particle swarm optimization algorithm is proposed to handle the problem with several objective functions, which originated from the study of bird predation behavior. The MOPSO algorithm first initializes a group of particles in the feasible solution space. Each particle represents a potential optimal solution to the optimization problem. The position, speed, and fitness value are used to represent the characteristics of the particles. The particle moves in the solution space, and the individual position is updated by tracking the individual extremum and the group extremum. In [32], a genetic algorithm is proposed to solve multiobjective optimization problems. This method is improved based on a genetic algorithm [45]. There are the following three key steps: 1) fast nondominated sorting operator design; 2) individual crowding distance operator design; and 3) elite strategy selection operator design. For the application of the abovementioned two algorithms and parameter adjustment, refer to [31], [32], as it will not be discussed in detail in this article.

D. Summary of the Proposed Strategy
A flowchart of the optimal search method proposed in this article is shown in Fig. 3. As seen from Fig. 3, the proposed search strategy is divided into two parts, which are the DEM-based target detection in the clutter-free region and search resources optimization.
1) First, according to the APAR location, the Cartesian coordinates DEM is projected to polar coordinates by (1)-(3). 2) Then, the platform-to-ground distances of all grids within the coverage of the current search beam position are calculated by (4)-(6). 3) Clutter-free detection region of the current beam position is obtained by searching for the smallest element in R i . 4) Detection threshold is set according to the noise power in the range to be detected, and then the target is detected by a square-law detector, whose performance is shown in (11). 5) From the cued information of the target, the weight w i of the current beam position is obtained by (12)

7) Dual-objective optimization function (30)-(31) is estab-
lished that maximizes the average target detection probability and minimizes the average target discovery time. 8) Multiobjective optimization problem in step 7 is solved by the MOPSO method to obtain the optimal beam dwell for each beam position. 9) Search airspace is scanned using the optimal beam dwell parameter.

IV. SIMULATIONS
In this section, we present our simulation results, including the solutions mentioned in previously sections.
The radar system parameters of the simulation in this article are set as follows. The search beamwidth is set to 2 • × 2 • . The azimuth and pitch range of the solid angle covered by the search airspace are [60 • , 120 • ] and [15 • , 60 • ], respectively. The normal azimuth and pitch direction of the radar array is [90 • , 25 • ]. The dwell time of a single illumination at each beam position T d = 5 ms, the radar transmit power P t = 10 KW, the antenna gain G = 40 dB, the wavelength λ = 0.1 m, Boltzmann's constant k B = 1.38 × 10 −23 J/K, the radar system temperature T e = 627.1 K, the total system loss G Loss = 7 dB, the bandwidth is B = 10 MHz, and the total radar system time resource T g = 3.5 s.

A. Formation of a Clutter-Free Region
Part of the DEM for the search airspace of interest at the current moment can be selected from the total stored ground elevation data according to the location of the airborne platform (0, vt, H), where H = 9000 m. The DEM of the airspace of interest at t is η t (x, y), as shown in Fig. 4. These DEM data are extracted from the open-source Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) dataset and cover the region near Dujiangyan City, Sichuan, China (103.0 • E, 31.0 • N). ASTER GDEM is acquired through stereo images, and the elevation accuracy can reach 20 m (95% confidence level). Its advantage is that it covers a wide area, reaching 95% of the land area [46]. ASTER GDEM is selected because of its satisfactory spatial resolution (30 m) and accuracy, as well as its world terrain coverage.
The purple-red sector indicates the airspace coverage of interest, and the triangle points A1 and A2 are the cued locations assumed in our simulations. Using (1)-(3), the DEM matrix η polar,t (R, θ) in polar coordinates is calculated, as shown in Fig. 5.
Using the DEM and (4)- (7), the optimal detection range gate, R i max , of each search beam position is calculated. Fig. 7 shows the results of the optimal range gates for all beam positions in sine space coordinates. The coordinate transformation relationship is deduced by [41] α = cos ϕ sin θ β = sin ϕ cos θ T − cos ϕ cos θ sin θ T (33) where θ T = arctan − cos ϕ min − cos ϕ max sin ϕ min − sin ϕ max cos θ max (34) denotes the inclination of the array. The distribution of circles in Fig. 7 represents the arrangement of the beam positions [33], and the inner coverage of each circle denotes the solid angle covered by the main lobe of each beam position. A total of 617 beams cover the entire search airspace. The intensity of the fill color of each circle indicates the optimal clutter-free range gate of each beam, and the color bar shows the relationship between the color and the range. According to the DEM data and optimal range gates from Figs. 5 and 7, a 3D simulation of the clutter-free region is shown in Fig. 6. It is seen that the search airspace is divided into a clutter-free region and a ground clutter region. To better distinguish the two regions visually, a separatrix is formulated by the MATLAB function "grid data" with the blue dots. These blue dots on the separatrix calculated by (4)-(7) denote the position of the optimal range gates for all 617 beams in Oxyz coordinates. Consequently, using the calculated optimal range gates R i max , the ith beam's SNR i and P d i with noise only is calculated by (8) and (11).

B. Target Prior Weight Distribution
Based on the ground prior information, we can easily obtain all known airport locations along the route. Suppose there are two cued locations in the search airspace at the same time. For cases with more cues, the weights are calculated in the same way. In this simulation, it is assumed that the cued locations A 1 and A 2 , as shown in Fig. 4. A 1 and A 2 are located at (−3.21, 11.75, 2.91 km) and (23.21, 74.01, 3.47 km), respectively. Their velocity are S A ≈ 0.53 km 2 and S B ≈ 0.24 km 2 , respectively. Using (12) and (13), the azimuth-pitch coordinates of A 1 and A 2 are (72.59 • , 5.95 • ) and (105.30 • , 31.64 • ), respectively. Then, the target appearance weights are derived using (12)- (17). From Fig. 8, the filled color intensity of each beam position indicates the target appearance weight. It is seen that the radius of the high-weight area corresponding to A 2 is larger than that of A 1 due to the higher velocity of A 2 . However, because A 2 is located at the edge of the search airspace, the high-weight area of A2 does not appear completely in Fig. 8.

C. Multiobjective Optimization and Algorithm Comparison
In this section, a detailed search performance comparison is conducted for the traditional average weighted sequential search method and the proposed optimization method with MOPSO and NSGA-II. First, under limited search time resources, the average target detection probability and the average target discovery time consumption are used to evaluate and analyze the performances of these three methods. In addition, a search performance factor is introduced to further analyze the search performances for varying SNR 0 .
The parameter configuration of the intelligent algorithms (MOPSO and NSGA-II) is given as follows. NSGA-II [32] can be implemented using the function named gamultiobj in the MATLAB [43]. The optimization parameters are set as follows: the number of variables nvar = 617, the variable lower boundary is 0, the upper boundary is 2L s T g /N Ω T d , the Pareto fraction is 0.3, the population size is 300, the number of generations is 200, the crossover fraction is 0.7, and the migration fraction is 0.2.
For the implementation of the MOPSO algorithm, refer to [31]. The optimization parameters are set as follows: the number of variables, bound of variables, and population size are the same as for NSGA-II. The repository size is 100, the number of iterations is 200, the inertia weight is 0.73, the self-learning factor is 1.5, the population-learning factor is 1.5, and the grid expansion parameter is 0.1.   Fig. 9. shows the average target detection probability using the traditional average weighted sequential search method and the proposed optimized search methods under varying search resource constraints L s . It can be seen that the proposed method provides a higher average target detection probability than that of the traditional one. In addition, when the search resources are not sufficient, it is shown that the MOPSO-based method achieves the highest target detection probability of the three methods. When the search resources increase, the average detection probability of the three methods gradually approaches 1. NSGA-II and MOPSO are still superior to traditional methods, but the gap gradually decreases. Fig. 10 shows the average target search time consumption using the traditional average weighted sequential search method and the proposed optimized search methods with varying search resource constraints L s . It is shown that the proposed method provides a shorter average target search time compared with the traditional method. When the search resources are not sufficient, the improvement in the search time consumption achieved by MOPSO is more obvious than that achieved by NSGA-II.
To further evaluate the performance of the method proposed with different SNRs, a search performance improvement factor where the improvement factor has a positive correlation with the average discovery time consumption and a negative correlation with the average detection probability. Fig. 11 shows the search performance of the three methods as a function of the signalto-noise ratio. The signal-to-noise ratio here refers to SNR 0 in Section I. It reflects the global signal-to-noise ratio of each beam position for the entire search airspace. From Fig. 11, when the SNR 0 increases, the performance factor of the three methods decreases, which means the search performance is improved. It can be seen that the proposed method can maintain its optimal performance at different SNR 0 , which verifies the robustness of the proposed search strategy. When the SNR 0 increases, the κ of the NSGA-II algorithm gradually approaches that of the traditional method, but the improvement of the MOPSO algorithm still performs well.

V. DISCUSSION
Considering the inherent measurement errors of DEM, this section includes the analysis of the influence of DEM error. In addition, the limitations of this method are discussed in the following.

A. Influence of the DEM Elevation Error
The inherently embedded measurement errors of the DEM, namely, the elevation error, can decrease the performance of the proposed method. The error directly decreases the target detection probability, which in turn affects the search performance.
Suppose an average elevation error of the used DEM is ΔH. When calculating the clutter-free detection range for a search beam position with a pitch angle of ϕ i , a large elevation estimate will lead to a narrowing of the clutter-free detection range. This makes the low-altitude target out of the detection range, resulting in missed detection. In addition, a low elevation estimate will lead to a widening of the clutter-free detection range, which makes the ground clutter enter the detection range and causes false alarms. If the elevation error is known, the projection of the error on R i max needs to be subtracted, then the corrected clutterfree range can be expressed as When ΔH > 0, the minimum searchable target range is narrowed, namely, the target whose height is less than ΔH from the real DEM cannot be detected. From [46], the average elevation accuracy of DEM data (e.g., ASTER DEM) is 20 m. If the bandwidth of the transmitted signal during target search and detection is B = 10 MHz [3], then the range resolution is Δr = c/2B = 15 m. In fact, there are differences in DEM errors in different regions. For ASTER DEM, the elevation error is less than 15 m in most areas [47]. Therefore, in most cases, the clutter detection range is reduced by no more than 1-2 range bins, which has little impact on the detection of low-altitude targets.
When ΔH < 0, the clutter-free detection range is widened. In this case, the unpredictable strong ground clutter enters the end of the detection range. If it is still considered that the detection range is free of clutter, the echo signal-to-noise ratio and detection probability will decrease. A worse single detection probability means that more accumulation time is required to search a target. Suppose that the average noise power is 20 dB, and the clutter power per range bin is 1.5 dB, the relationship of the target detection performance varying with the DEM errors for different pitch angles ϕ is shown in Fig. 12. It can be seen that when ΔH < 0, the detection performance decreases severely for each ϕ. To avoid this detection probability decrease, we set a protection threshold according to the maximum elevation error (ΔH) max of the selected DEM. The modified clutter-free region is expressed as [R i min , R i max − (ΔH) max cos ϕ i ].

B. Limitations of the Method
The low-altitude search capability decreases as the detection range increases. Since the radar scanning beam width is constant, the beam coverage footprint increases as the pitch angle decreases. When calculating the clutter-free region, the excluded airspace where the target may exist is increased, as shown in Fig. 13. When the pitch angle is large, the excluded low-altitude airspace is small, as shown in the red areas 1 and 2 in Fig. 13. In this case, the minimum target detection height is slightly higher than the DEM, and the low-altitude target search performance is good. When the pitch angle is small, the excluded low-altitude airspace is large, as shown in the red areas 3, 4, and 5 in Fig. 13. In this case, the distribution of the minimum detection heights of targets at different ranges in the same search beam is quite different. There may be targets that are much higher than the ground but are excluded and cannot be detected, so the search performance of low-altitude targets will be affected. The detectable range of low-altitude targets is a sawtooth shape that gradually changes with range on an azimuth profile, such as the red areas 1, 2,..., 3, 4, 5, etc. in Fig. 13. To quantitatively analyze the trend of the detectable range of low-altitude targets, the average ground height of a certain area is set as H DEM , and the flight height of the aircraft is H. The maximum height loss H for different pitch angles is defined as The maximum height loss represents the maximum height of airspace that is removed by the proposed method with a clutterfree detection range but where there may exist low-altitude targets. The same ASTER DEM data were used as the simulation with an average ground height of 1754 m. Other parameters are consistent with the simulation. The maximum height loss for different pitch angles is shown in Fig. 14. It can be seen that when the pitch angle is small, the maximum height loss is too large and unacceptable. When the pitch angle is greater than 20 • , the maximum height loss is less than 500 m. And as the pitch angle increases, the maximum height loss becomes smaller. In summary, our method works well for large down-view angles (pitch angles greater than 20 • ), but the performance degrades drastically for smaller down-view angles (pitch angles less than 20 • ).

VI. CONCLUSION
This article proposes an optimal search strategy using the DEM-aided approach. First, clutter-free airspace is calculated by introducing the DEM to improve the target detection probability. Then, based on the detection result, the search problem is transformed into a constrained multiobjective optimization problem that can simultaneously minimize the search time and maximize the target detection probability. Compared with the traditional method, the proposed search strategy improves the detection probability and reduces the time cost for low-altitude targets.