3.1 Algorithm flow

The first problem in multi-target tracking is to start the target automatically. After target initiation, the data association algorithm based on a particle filter associates the measurements at the current time with the possible events one by one, including a new generation, continuation and extinction of a target, while simultaneously eliminating the influence of clutter. After the measurements have been correlated with the estimated state of all known targets and the initial state of new targets, the updated values of the current state for all existing targets are obtained using a Kalman state method. Finally, all the Kalman filtering results are smoothed to obtain a smooth track for each existing target. The algorithm flow is shown in Figure 2.

Figure 2. Multi-target automatic initiation and tracking algorithm.

Every target has a lifecycle process from birth to death, so the number of targets is always changing. The data association algorithm can be applied to judge whether a measurement belongs to an existing target, a new target or clutter, thereby simplifying the complex radar multi-target tracking problem to the common single-target tracking problem. At the same time, the lifecycle of each target is recorded to obtain real-time statistics on the total number of existing targets.

3.2 Data association

Multi-target data association based on a particle filter algorithm can be regarded as a generalisation of the Multiple Hypothesis Tracking (MHT) algorithm. Each particle represents a different data association hypothesis for target initiation, known target or target death, which reduces the complexity of the MHT algorithm. For each particle, the algorithm associates the measurement with several possible events, calculates the prior probability P _p of each event and the event correlation similarity P _l and then estimates the correlation probability P of each event. The correlation probability models of all the possible events are listed in Section 3.3.

The algorithm gives each of the N particles the same initial weight $ \{ w_{k}^{(i)}=1/N, $ $ i=1,{\ldots},N \} $ , while the correlation probability P is calculated as

(1) \begin{equation} P = {P_{\mathrm{p}}} \cdot {P_{\mathrm{l}}}\end{equation}

where P _p is the prior probability and P _l is the event correlation similarity. To estimate the measurement correlation, the Monte Carlo sampling method is used to determine the event category.

The weight of each sample is calculated as

(2) \begin{equation}w_k^{(i)} = w_{k - 1}^{(i)} \cdot P_{j,{\mathrm{p}}}^{} \cdot {P_{j,{\mathrm{l}}}}/P_j^{}{\mathrm{ }}\end{equation}

where P _j,p , P _j,l and P _j respectively represent the normalised prior probability, event correlation similarity and normalised correlation probability of event j. Note that the index j represents the event category while the index k represents the particle number.

Then, the weight of each sample is normalised as

(3) \begin{equation}w_k^{*(i)} = \frac{{w_k^{(i)}}}{{\sum\nolimits_{i=1}^N {w_k^{(i)}} }}\end{equation}

The number of valid samples n is estimated as

(4) \begin{equation}n = \frac{1}{{\sum\nolimits_{i=1}^N {{{\left( {w_k^{*(i)}} \right)}^2}} }}\end{equation}

If n is too small, resampling is performed.

3.3 Probability estimation of possible event correlation

In the absence of any prior knowledge, the correlation results of measurements depend on the calculation of the probability of occurrence of all possible events, such as clutter, new targets, known targets and target extinction. Firstly, p _b is set as the probability of a new target, which means that the current measurement is the starting state of a new target. In this paper, an active hotspot of bird targets is taken as the reference point for target track initiation, and track initiation is achieved by estimating the degree of correlation between the measurements and the hotspot, which reduces the complexity of the tracking algorithm and improves its efficiency and engineering applicability. In addition, p _d is taken as the target death probability, which is estimated by a gamma distribution function as

(5) \begin{equation}gamma(\unicode[Times]{x03B8} |\unicode[Times]{x03B1} ,\unicode[Times]{x03B2} ) = {\unicode[Times]{x03B8}^{\unicode[Times]{x03B1} - 1}}\frac{{{\unicode[Times]{x03B2}^{\unicode[Times]{x03B1}} }{e^{\frac{{ - \unicode[Times]{x03B8} }}{\unicode[Times]{x03B2} }}}}}{{\Gamma (\unicode[Times]{x03B1} )}}\end{equation}

where $\unicode[Times]{x03B8} $ is the interruption time of the association between the target track and measurements. The target death probability means the probability of the end of the target trajectory. The probability of the target disappearing increases with the extension of the interruption time. Note that the slope of the function changes with different values of $\unicode[Times]{x03B1} $ and $\unicode[Times]{x03B2} $ .

Therefore, the probability models of all possible measurement-related events can be established and summarised as follows:

1. (1) Clutter, no target death

(6) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = (1 - {p_{\mathrm{b}}}) \cdot {(1 - {p_{\mathrm{d}}})^n} \cdot cp\\{P_{\mathrm{l}}} = cd\end{array}\end{equation}

where cp is the prior probability of clutter, cd is the probability of clutter distribution and the current number of targets is n.

1. (2) Clutter, one target death

(7) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = \left\{ \begin{array}{l}(1 - {p_{\mathrm{b}}}) \cdot {p_{\mathrm{d}}} \cdot {(1 - {p_{\mathrm{d}}})^{n - 1}} \cdot cp{\mathrm{, }}\;n\,>\,1\\(1 - {p_{\mathrm{b}}}) \cdot {p_{\mathrm{d}}} \cdot cp{\mathrm{, }}\;n=1\end{array} \right.\\{P_{\mathrm{l}}} = cd\end{array}\end{equation}

1. (3) Known target j, no target death

(8) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = \frac{{(1 - {p_{\mathrm{b}}}) \cdot {{(1 - {p_{\mathrm{d}}})}^n} \cdot (1 - cp)}}{n}\\{P_{\mathrm{l}}} = {p_{{\mathrm{tlh}}}}\end{array}\end{equation}

where p _tlh is the target association similarity, which satisfies a Gaussian distribution with T _j and S _j as parameters

(9) \begin{equation}{p_{{\mathrm{tlh}}}} = N({{\textbf{y}}_k}|{{\textbf{T}}_j},{{\textbf{S}}_j})\end{equation}

where T _j and S _j are the estimated state and covariance matrix of target j, and ${{\textbf{y}}_k}$ is the measurement.

1. (4) Known target j, one target death

(10) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = \frac{{(1 - {p_{\mathrm{b}}}) \cdot {p_{\mathrm{d}}} \cdot {{(1 - {p_{\mathrm{d}}})}^{n - 1}} \cdot (1 - cp)}}{{n - 1}}{\mathrm{, }}n\,>\,1\\{P_{\mathrm{l}}} = {p_{{\mathrm{tlh}}}}\end{array}\end{equation}

1. (5) New target initiation, no target death

(11) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = {p_{\mathrm{b}}} \cdot {(1 - {p_{\mathrm{d}}})^n}\\{P_{\mathrm{l}}} = {p_{{\mathrm{ilh}}}}\end{array}\end{equation}

where ${p_{{\mathrm{ilh}}}}$ is the target initial association similarity.

1. (6) New target initiation, one target death

(12) \begin{equation}\begin{array}{l}{P_{\mathrm{p}}} = \left\{ \begin{array}{l}{p_{\mathrm{b}}} \cdot {p_{\mathrm{d}}} \cdot {(1 - {p_{\mathrm{d}}})^{n - 1}}{\mathrm{ }}\!\quad\qquad n>1\\{p_{\mathrm{b}}} \cdot {p_{\mathrm{d}}}{\mathrm{ }}\qquad\qquad\qquad\qquad n=1\end{array} \right.\\{P_{\mathrm{l}}} = {p_{{\mathrm{ilh}}}}\end{array}\end{equation}

Based on the estimation results of the correlation probability of these six kinds of possible event, management of the whole lifecycle of each target from birth to death can be realised.

3.4 Kalman filtering and smoothing

The essence of the Kalman filter is to estimate the future state of a target according to its current state, and to modify it based on the measured value. The Kalman filter can be divided into two steps. The prediction part estimates the next state of the system according to the measurement value in the previous step, while the update part estimates the current state of the system according to the measurement value at the current time. The prediction part is completed using the equation

(13) \begin{equation}\begin{array}{l}{\textbf{m}}_k^ - = {{\textbf{A}}_{k - 1}}{{\textbf{m}}_{k - 1}}\\{\textbf{P}}_k^ - = {{\textbf{A}}_{k - 1}}{{\textbf{P}}_{k - 1}}{\textbf{A}}_{k - 1}^T + {{\textbf{Q}}_{k - 1}}\\{{\textbf{T}}_k} = {{\textbf{H}}_k}{\textbf{m}}_k^ - \\{{\textbf{S}}_k} = {{\textbf{H}}_k}{\textbf{P}}_k^ - {\textbf{H}}_k^T + {{\textbf{R}}_k}\end{array}\end{equation}

where ${\textbf{m}}_k^ - $ and ${\textbf{P}}_k^ - $ are the mean and variance of the estimated state before measurement at time k. A, Q, R and H are the transfer matrix, process noise matrix, measurement noise matrix and measurement model matrix of the dynamic model, respectively. T and S are the mean and variance of the estimated measurement state. The update part is completed using the equation

(14) \begin{equation}\begin{array}{l}{{\textbf{v}}_k} = {{\textbf{y}}_k} - {{\textbf{T}}_k}\\{{\textbf{K}}_k} = {\textbf{P}}_k^ - {\textbf{H}}_k^T{\textbf{S}}_k^{ - 1}\\{{\textbf{m}}_k} = {\textbf{m}}_k^ - + {{\textbf{K}}_k}{{\textbf{v}}_k}\\{{\textbf{P}}_k} = {\textbf{P}}_k^ - - {{\textbf{K}}_k}{{\textbf{S}}_k}{\textbf{K}}_k^T\end{array}\end{equation}

where ${{\textbf{m}}_k}$ and ${{\textbf{P}}_k}$ are the mean and variance of the estimated state after the measurement at time k, v is the measurement modification and K is the filtering gain, which defines the degree of the modification of the estimated value.

The process of Kalman smoothing is similar to that of Kalman filtering, the difference being that the recursion of filtering is forward, while the recursion of smoothing is backward. Higher precision can be obtained after Kalman smoothing compared with only Kalman filtering.

4.1 Simulation of target track initiation

The state of the moving target is expressed as

(15) \begin{equation}{{\textbf{x}}_k} = {\left( {\begin{array}{*{20}{c}}{{x_k}}\quad{{y_k}}\quad{{{\dot x}_k}}\quad{{{\dot y}_k}}\end{array}} \right)^{\mathrm{T}}}\end{equation}

where $({x_k},{y_k})$ and $({\dot x_k},{\dot y_k})$ are respectively the position and speed of the target in the two-dimensional rectangular coordinate system at time k.

Figure 3. Simulation data of target tracking initiation.

The discrete dynamic model of a target can be expressed as a linear time-invariant construction equation

(16) \begin{equation}{{\textbf{x}}_k} = \left( {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}1& 0&{\Delta t}&0\\0&1&0&{\Delta t}\\0&0&1&0\\0&0&0&1\end{array}} \right){{\textbf{x}}_{k - 1}} + {{\textbf{q}}_{k - 1}}\end{equation}

where ${{\textbf{q}}_{k - 1}}$ is the white noise of the discrete Gaussian process, and the moment characteristic is

(17) \begin{equation}\begin{array}{l}E[{{\textbf{q}}_{k - 1}}] = {\textbf{0}}\\E[{{\textbf{q}}_{k - 1}}{\textbf{q}}_{k - 1}^T] = \left( {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\textstyle{1 \over 3}}\Delta {t^3}}&0&{{\textstyle{1 \over 2}}\Delta {t^2}}&0\\0&{{\textstyle{1 \over 3}}\Delta {t^3}}&0&{{\textstyle{1 \over 2}}\Delta {t^2}}\\{{\textstyle{1 \over 2}}\Delta {t^2}}&0&{\Delta t}&0\\0&{{\textstyle{1 \over 2}}\Delta {t^2}}&0&{\Delta t}\end{array}} \right)q\end{array}\end{equation}

where the time step is set to $\Delta t = 0.01$ and the power spectral density of the process noise is set to $q = 0.001$ . Note that the time unit in the simulation is set to seconds.

Figure 3 shows the simulation tracking process of the target in cluttered environments with different parameters. There are 50 simulation steps, with the target measurements shown by open circles and the clutter by crosses. The amount of clutter in each scanning period follows a Poisson distribution with a mean value of $\lambda$ . Taking (0, 0) as the starting reference point for the target, the distance between the starting point of the target track and the reference point is d, the absolute value of the starting speed of the target is 1 and the starting direction is randomly distributed in the range from 0° to 360°. Figure 3(a) shows a global view and Figure 3(b) shows a local enlarged view, where the target initiation point is labelled with an open circle. In this experiment, the parameters of the tracking algorithm are set to d = 0.2 and $\lambda$ = 0.5, and the particle number is N = 50. Figure 3(c) shows a global view of the target track in another simulation experiment, while Figure 3(d) shows a local view of the target trajectory. In this experiment, the parameters are set as d = 3 and $\lambda$ = 1 with a particle number of N = 50. Comparing the two target tracking results shown in Figure 3, it is clear that the target track starting point in Figure 3(a) is close to the reference point (d = 0.2), leading to no delay in track starting, while the target track starting point in Figure 3(c) is relatively far from the reference point (d = 3), leading to three cycles of delay in track starting.

Figure 4. Target initiation delay under conditions with different amounts of clutter.

Figure 5. Target initiation delay under different initiation conditions.

Setting different clutter environments (Λ value) and target track starting conditions (d value), the average target tracking results obtained from 100 Monte Carlo simulation experiments are shown in Figures 4 and 5, and the stability of the target tracking and the timeliness of track initiation under different parameter settings are investigated.

Figure 4 shows the change of the initial delay period of the target track with the cluttered environment when the distance between the initial position of the target track and the reference point is varied (d = 0.5, 1.5, 2.5 and 3). It is clear that, as the amount of clutter is increased (from Λ = 0 to 1.5), the fluctuation of the delay period of target track initiation will also increase, and the increase of the clutter density will affect the timeliness of target correlation, resulting in a delay in track initiation. The fluctuation of the delay period curve may be related to the number of Monte Carlo experiments, and its fluctuation will decrease with increasing number of experiments. The experimental results shown in Figure 4 reveal that, when d < 3 (d = 0.5, 1.5 or 2.5), the range of the delay period of track initiation does not change greatly, and it is generally less than 1, which is obviously better than the logic method for track initiation in three scan cycles^{(Reference He, Xiu, Zhang and Guan12)}. Even when d = 3, the value of the initiation delay period curve is still less than 3 on the whole.

The logical initiation method is the most widely used engineering method for target tracking, usually taking the nearest measurement from the gate centre as the effective measurement and completing target initiation in three or four scanning cycles^{(Reference He, Xiu, Zhang and Guan12)}. The logic initiation method in this paper refers to the three-cycle logic initiation method, which can realise the start of a target trajectory based on target measurements in three scanning cycles under the ideal condition with no clutter. In fact, in a cluttered environment or in the case of target measurement loss, the number of cycles required to achieve target initiation by the logic initiation method may be further increased.

Table 1 Comparison of delay periods of target initiation

Figure 5 shows the variation of the delay period of the target track initiation with the distance between the track initiation point and the reference point in environments with different amounts of clutter (Λ = 0.1, 0.5, 1.5 and 2). Similar to the conclusion related to Figure 4, with increasing relative distance (from d = 0 to 2.5), the delay period curve of each track initiation shows some fluctuation, but the overall value range changes little. In addition, when the amount of clutter is small (Λ = 0.1), the delay period of track initiation is close to zero, with almost no delay. When the value of Λ is gradually increased, the delay period of track initiation increases, but the overall delay period does not exceed 3.

Table 1 presents the average value of the delay period (T _D) of track initiation for different values of Λ, all of which are less than 2, being shorter than the logic method, which requires three scanning cycles to complete track initiation. When the value of Λ is set low enough, there is barely any delay in target initiation. When Λ = 2, the average value of the initiation delay period is 1.8583, which is caused by the clutter in the simulation data. In the simulation model established in this paper, part of the clutter data is used to replace the target measurement data, which leads to a certain degree of initiation delay of the target track.

To further verify the advantages of the proposed algorithm, two additional indexes are considered: the average number of true track initiation and the average number of false track initiation, represented by N _In ^T and N _In ^F. In 100 experiments of target tracking in environments with different amounts of clutter, only one real track is generated in each simulation. For example, in 100 simulations, if the track is started correctly 99 times while a false track is generated 3 times, then N _In ^T = 0.99 and N _In ^F = 0.03. The average experimental results for environments with different amounts of clutter are presented in Table 1, where all the performance indexes of the proposed algorithm are better than those of the three-cycle logic method. In an environment with more clutter (Λ = 2), the current algorithm can still maintain the correct target initiation (N _In ^T = 1) with fewer false tracks (N _In ^F = 0.02), but the average number of correct target initiations using the logic method will decrease (N _In ^T = 0.95) while the false alarm tracks will increase (N _In ^F = 0.06) significantly.

4.2 Simulation analysis of multi-target lifecycle

As shown in Figure 6, the total number of simulation steps designed in this experiment was 500 with a time step of 0.01, when three target tracks starting in sequence are simulated. The lifecycle of target 1 is t = 0 to 5. Starting from (–0.25, –1.5), it moves at a constant speed of (0, 1) from t = 0 to 2, makes a left turn from t = 2 to 4 and moves at a constant speed (–1, 0) from t = 4 to 5 until the end of its track. The lifecycle of target 2 is from t = 1 to 4.5, starting from (0, –1.5) and moving at a constant speed of (0, 1) throughout its lifecycle until the end of its track. The lifecycle of target 3 is from t = 0.5 to 4.5. Starting from (0.25, –1.5), it moves at a uniform speed of (0, 1) from t = 0.5 to 3, makes a right turn from t = 3 to 4, and moves at a uniform speed of (1, 0) from t = 4 to 4.5 until the end of the track. The real moving track of the target is shown by a solid line in Figure 6, where an open circle represents a target measurement with certain Gaussian noise added and a cross represents clutter with a parameter value of Λ = 0.1, being randomly distributed in the area of [–3, 2] × [–2, 3]. In the data association, the number of particles is set to N = 50, p _b is set to 0.01, cp is set to 0.02, cd is set to1/16 and the starting reference point of the target track is set to (0, 0).

Figure 6. Multi-target tracking simulation experiment.

Figures 6 and 7 compare the tracking results and target number estimation results for three targets with different parameters set in the target extinction model, both achieving timely initiation and stable tracking of all targets in the cluttered environment. In the estimation of the tracking end time, when the parameters in the target death model are set to $\unicode[Times]{x03B1} = 2$ and $\unicode[Times]{x03B2} = 0.5$ , the judgement of the target track initiation time and death time is basically consistent with the real situation, as shown in Figures 6(a) and 7(a). When the parameters in the target death model are set to $\unicode[Times]{x03B1} = 2$ and $\unicode[Times]{x03B2} = 2.5$ , the judgement of the target death time is delayed, and the tracks of targets 2 and 3 fail to finish at a time at t = 4.5, resulting in deviation of the estimation results from the target data, as shown in Figures 6(b) and 7(b). It is demonstrated that the parameter $\unicode[Times]{x03B2} $ controls the target death probability by changing the correlation interruption time between the target track and the measurements, and increasing $\unicode[Times]{x03B2} $ leads to a decrease of the target death probability.

Figure 7. Target number estimation.

5.1 System platform

Figure 8 shows the experimental avian radar system developed by China Academy of Civil Aviation Science and Technology (CAST). The system uses a solid-state transmitter to replace the original magnetron of the S-band navigation radar, and introduces the signal processing technologies of pulse compression, pulse Doppler and constant false alarm rate. The algorithm for bird target detection and tracking is also improved continuously. So far, the system has collected and accumulated a large amount of bird radar target information over more than two years at several airports. In this paper, two bird observation cases from different locations are described, where the analysis of each example is based on observation data taken over 10 days. Due to the unsatisfactory performance of the avian radar in bad weather conditions, all the bird data processed in this paper were collected under conditions without precipitation.

Figure 8. CAST experimental avian radar system.

5.2 Case 1

The resident birds around the airport usually follow the rule of leaving their nests in the morning and returning home in the evening. By using the avian radar to obtain the bird information around an airport in southwest China, some statistical analysis could be carried out on the activity rhythm of the airport resident birds to determine the number and distribution of birds around the airport. In this section, based on the bird target data of an airport obtained by the CAST experimental avian radar system, the basic steps and examples of statistical analysis of bird target number around the airport are presented. Figure 9 shows a flowchart of the statistics and analysis of the bird situation data to summarise and clarify the innovation presented herein, revealing that the proposed tracking algorithm is especially suitable for the case of automatic initiation of multiple targets to determine the starting reference point of a flock of targets, as in this case of a flock of birds leaving their nests

Figure 9. Analysis flowchart of bird distribution situation.

5.2.2 Grid distribution statistics of bird activities

As shown in Figure 10, the surveillance area of L × L of the avian radar was divided into a series of rectangular grids, each having an area of l × l, both in units of m². In this example, the radius of coverage of the avian radar is 3km with L = 6000m and l = 100m, and the total number of grids in the monitoring range is 60 × 60. In the figure, the upper-left corner is set as the origin of the coordinate system (0, 0), while the X-axis is horizontally to the right and the Y-axis is vertically downwards.

Figure 10. Hotspot statistics of bird activities around the airport in case 1.

During the period of bird departures (07:00 to 08:00) around the airport, the grids with bird target measurements exceeding a certain threshold value are marked in colour. In this case, after ten scanning cycles of radar data accumulation, grids with more than ten bird target measurements are determined to be the cells representing the hotspot. Using the marked grid distribution map, binary image processes such as erosion and dilation are carried out to eliminate isolated points, thereby forming a connected area called the hotspot, as shown in Figure 10. In combination with the bird survey results around the airport, the bird habitat can be determined. The centre point of the hotspot area is taken as the reference point of the target initiation. In this case, the reference point coordinate is (4750m, 2650m).

Figure 11. Flight path distribution of birds in case 1.

Figure 12. Bird population estimation of case 1.

5.2.3 Statistical analysis of bird population

Figure 11 shows the flight trajectory distribution of birds leaving the habitat in the hotspot area in the early morning of a certain day, including a global view within the radar coverage and a local enlarged view of the bird flight paths. The bird habitat is located to the northeast, outside the flight area. The process of birds leaving their nests, as shown in Figure 11, occurs between 07:00 and 08:00, lasting for about 1min. The population number is about 20, and the flight direction is generally southwest, when the birds usually enter the flight area to feed.

Figure 12 shows the estimated results for the number of birds leaving the nest in the habitat between 07:00 and 08:00 on different days, when the target tracking time lasts about 70s. In the data association, the number of particles is set to N = 50, p _b is set to 0.01, cp is set to 0.02, cd is set to 1/16 and the reference point of target track initiation is set to (4750, 2650). The parameters in the target death model are set to $\unicode[Times]{x03B1} = 2$ and $\unicode[Times]{x03B2} = 0.5$ . The maximum value of each bird target number estimation process is taken as the estimated value of the population number, thus the first estimate is 23 and the second estimate is 24.

Table 2 presents the bird population statistics of the habitat for 10 consecutive days. It is clear that the bird population number automatically counted by the algorithm is almost consistent with the results of artificial observations completed by professional bird watchers with the help of telescopes.

Table 2 Estimation results of bird population in case 1

5.3 Case 2

The processing and analysis results of the bird situation in this example are based on radar data obtained during the application of the CAST experimental avian radar system at a general airport in northeast China, trying to reveal the general number of birds around the airport and any regularity in their activities.

Similar to case 1, the surveillance area of the avian radar is again divided into a series of rectangular grids with coverage radius of 2km. As shown in Figure 13, the upper-left corner is also set as the origin of the coordinate system (0, 0), while the X-axis is horizontally to the right and the Y-axis is vertically downwards. Based on a large amount of bird data accumulated by the avian radar and combined with artificial observation results, it is confirmed that two bird activity hotspots in the northern area outside the airport are bird habitats, as marked in Figure 13. Moreover, in this case, the centre points of the two hotspots are respectively taken as the reference points of the target initiation: (630m, 835m) and (3430m, 410m).

Figure 13. Hotspot statistics of bird activities around the airport in case 2.

Figure 14. Flight path distribution of birds in case 2.

Figure 15. Bird population estimation of case 2.

Figure 14 shows the flight path distribution of the birds leaving the habitat in this case and a local enlarged view of the two hotspots. The two bird flocks usually left the habitat between 07:00 and 08:00 in the morning, when the flight direction of birds from hotspot 1 was scattered in the easterly direction, while that of hotspot 2 was mainly concentrated in the southwest direction. Figure 15 shows the estimated results for the number of birds in the two flocks when leaving their nest in the habitat on a certain day, lasting about 60s. The parameter settings used for target initiation and tracking are the same as in case 1. The maximum values in the process of the two target number estimations are taken as the number of bird targets in the two flocks, viz. 24 and 15, respectively. Table 3 presents a comparison of the bird populations in the two hotspots over ten consecutive days, including the results estimated using the avian radar and artificial observations. According to the results of the artificial bird survey and the avian radar observation, there are about 25 birds in hotspot 1 on the northwest side of the airport, and their activity range is mainly in the flat area on the north side of the airport. There are about 15 birds in hotspot 2 located in the northeast side of the airport. Note that the birds in hotspot 2 often fly over the airport runway and enter the mountainous area to the south side of the airport, which causes a high risk of bird strikes.

Table 3 Estimation results of bird population in case 2