Successive rejection track-before-detect method based on EM algorithm

Abstract

The invention discloses a successive elimination track-before-detect method based on an EM algorithm, and belongs to the field of radar weak target detection, and the method comprises the following steps: 1, carrying out the initialization setting; step 2, current frame data pre-detection: CFAR low-threshold pre-detection is carried out on the data of the current time frame of the beam in the same direction; step 3, current frame trace point condensation: performing centroid estimation on trace points continuously distributed in a distance direction by using amplitude weighting to obtain a condensed trace point zi; 4, sliding window multi-frame trace point: forming a multi-frame trace point z = [zi-F + 1, zi-F + 2,..., zi] T according to the length of the sliding window, wherein F is the frame number of the sliding window; step 5, motion parameter initialization based on a halton sequence; step 6, performing weighted iteration single-target motion parameter estimation based on an EM algorithm; step 7, multi-target parameter estimation of used trace points is successively eliminated; and 8, repeating the steps 2-7 until the data frame is ended, and the method can effectively suppress the false alarm and stably detect the flicker target at the same time, so that the detection probability is remarkably improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of radar weak target detection, and particularly relates to a sequential rejection detection before tracking method based on the EM algorithm. Background Art

[0002] Traditional radar target detection detects targets by comparing the amplitude information of a single frame with the statistical background under the assumption of a Gaussian white noise background. When the clutter conditions and target scintillation occur, the false alarm rate is high and the target detection performance is unstable. The detection before tracking technology uses the continuous information of the target state before detection to find possible target points in multiple frames of data, so as to achieve the purpose of detecting weak targets and suppressing false alarms.

[0003] The current detection before tracking technologies mainly adopt three methods: the Hough transform method [1], the particle filter method [2], and the dynamic programming method [3].

[0004] The Hough transform detection before tracking maps points in space to the parameter space for accumulation and then detection. The detection performance depends on the granularity of the parameter space division and the setting of the threshold. When the parameters are not set properly, the detection performance is poor. The particle filter detection before tracking uses a weighted sampling particle to approximate the target probability distribution, which depends on the local approximation of the current state and is a local method. The detection performance for scintillating targets is unstable.

[0005] The dynamic programming detection before tracking constructs the global optimal solution by solving the optimal solution in stages through constructing a cost function of the target state over time. The method based on stages is prone to be inconsistent with the global state of the target, resulting in the phenomenon of batch switching under the condition of multiple crossing targets.

[0006] References: [1] Xu Luping, Xu Na, Yang Sheng, Yan Bo, etc., Detection before tracking method based on three-dimensional Hough transform [P], CN109557532A.

[0007] [2] Huang Yong, Pei Jiazheng, Liu Ningbo, Chen Baoxin, etc., Detection before tracking method and system based on LMB particle filter [P], CN110865343A.

[0008] [3] Yang Xiaobo, Ku Feilong, Xu Longxiao, Fu Lingzhi, etc., Dynamic programming detection before tracking method for moving platform radar with motion compensation search, CN109917373A. Summary of the Invention

[0009] Aiming at the above problems, the purpose of the present invention is to overcome the disadvantages of the poor detection performance of the existing technology for scintillating weak targets, and provide a new detection before tracking method based on the EM algorithm, which can effectively and stably detect scintillating weak targets.

[0010] In order to achieve the above object, the technical solution adopted by the present invention is: a tracking method before detection based on the EM algorithm, comprising the following steps: Step 1, Initialization: Set the target speed range, accumulated frame number, frame detection ratio threshold, average error threshold, and distance threshold; Step 2: Pre-detection of current frame data: Perform CFAR low threshold pre-detection on the data of the current time frame of the same azimuth beam; Step 3: Current frame trace condensation: Use amplitude weighting to estimate the centroid of the traces continuously distributed in the distance direction to obtain the condensed traces zi; Step 4, sliding window multi-frame traces: generate multi-frame traces z=[zi-F+1, zi-F+2,…,zi]T according to the sliding window length, where F is the number of sliding window frames; Step 5: Initialize motion parameters based on the halton sequence: Use the halton sequence to generate multiple sets of linear motion parameters X, calculate the error and average error between the multi-frame traces and the theoretical point of each set of parameters, and take the motion parameter with the smallest average error as the initial target motion parameter a corresponding to the current multi-frame traces; Step 6: Weighted iterative single target motion parameter estimation based on EM algorithm: According to the error and average error between the trace and the theoretical point, the weight of the point is calculated according to the Gaussian distribution, and the weight matrix W=diag(w) is generated. Recalculate the motion parameters until the error between the two parameters is less than the threshold, and output the corresponding target motion parameter a; Step 7: Multi-target parameter estimation by removing used points one by one: According to parameter a, the points with the closest distance to the theoretical point trace and less than the error threshold are selected in each frame time. If the number of points is greater than the proportional threshold, new target motion parameters are formed; otherwise, no new target motion parameters are formed; the selected points are removed from the point trace, and the remaining number of points is compared with the total number of frames. If it is greater than the proportional threshold, step 5 is repeated; otherwise, the final multi-target parameters of the current time frame are output; Step 8: Repeat steps 2 to 7 until the data frame ends.

[0011] The step 2 is specifically as follows: for each distance unit, a protection unit of g units is used, and the noise background of each of the M points on the left and right is averaged. If the amplitude of the current unit is higher than the threshold than the background, it is judged that there is a target and the position and amplitude are output. Otherwise, it is judged that there is no target.

[0012] The step 3 is specifically as follows: for each frame of detected points, condense them along the distance direction: if there are detection points in continuous distance units, weight them according to the amplitude to form a point at the center position; if the points are not continuous, directly output their positions.

[0013] Step 5 specifically includes the following: Step 5.1: Generate K sets of Halton sequences with a dimension of 2; and generate theoretical points of multiple sets of parameters according to Y = CX, where Y is multiple sets of target theoretical points formed by multiple sets of motion parameters, C is a matrix constructed by the vector t formed by the time of the tracks and the all-ones sequence 1, C = [t, 1], t is the frame time sequence of the tracks, and 1 is an all-ones sequence of the same length as t, and the length is the number of tracks N; Step 5.2: Generate multiple sets of target motion parameters X = [v, r]T according to the target distance, speed range, and Halton sequence, where v = v0 + vlmt * halton(k, 1), r = r0 + rlmt * halton(k, 2), k = 1, 2,..., K; Step 5.3: For each set of parameters Xi, calculate the error dj = Yj - z between multiple frames of tracks and the theoretical track CXi and the average error , where Yj is the j-th column of Y; Step 5.4: Take the parameters with the minimum average error as the initial target motion parameters corresponding to the current track , , is the k-th column of the K sets of linear motion parameters X.

[0014] Step 6 specifically includes the following: Step 6.1: Calculate the error d = z - Ca between the track z and the theoretical track Ca of the parameter a and the average error , where N is the total number of tracks z; Step 6.2: Generate the weights of the tracks according to the Gaussian distribution , n = 1, 2,..., N, is the n-th component of d, and the weight matrix W = diag(w), w = [w1, w2,..., wN]T; Step 6.3: Let , and recalculate the motion parameters ; Step 6.4: Calculate the theoretical target motion point zt = Ca under the new parameters, and recalculate the error d, average error and the weight ; Step 6.5: Calculate the difference e between the parameters before and after ; Step 6.6: If e is greater than the threshold, continue with Step 6.1; otherwise, exit and output a.

[0015] Step 7 specifically includes the following: Step 7.1: Select the points that are closest to the theoretical plot points and less than the error threshold within each frame time according to the parameter a output in Step 6; Step 7.2: If the number of selected points is greater than the ratio threshold, form new target motion parameters; otherwise, do not form new target motion parameters; Step 7.3: Remove the selected points from the plot points and replace z with the remaining plot points; Step 7.4: Compare the number of remaining plot points z with the total number of frames. If it is greater than the ratio threshold, go to Step 6 for weighted iterative single-target motion parameter estimation based on the EM algorithm; otherwise, output the final multi-target parameters.

[0016] Compared with the prior art, the technical solution adopted by the present invention has the following beneficial effects: The present invention calculates the probabilities that each pre-detected point belongs to each moving target in the E-step (expectation step) of the EM algorithm, uses these probabilities as weights, and then calculates the motion parameters of each target in the M-step (expectation maximization step). Such iteration is performed until the multi-target parameters converge.

[0017] Since this method performs parameter estimation by jointly using multi-frame data, it is a global algorithm. Also, during the iteration process, it automatically and accurately estimates the target parameters using the attributed points without searching in the parameter space. Therefore, it can stably detect flickering targets while effectively suppressing false alarms, and the detection probability is significantly improved. Description of the Drawings

[0018] Figure 1 It is a schematic diagram of the pre-tracking method based on the EM algorithm for successive elimination detection in this embodiment.

[0019] Figure 2 It is a schematic diagram of the motion trajectories of two targets set in this embodiment.

[0020] Figure 3 It is a schematic diagram of the echo amplitudes of two targets set in this embodiment.

[0021] Figure 4 It is a schematic diagram of the plot points condensed after echo detection in this embodiment.

[0022] Figure 5 It is a schematic diagram of multi-target detection for successive elimination of used plot points in this embodiment. Detailed Embodiment

[0023] The following further explains the technical solution of the present invention in conjunction with the drawings of the specification and specific embodiments.

[0024] As Figure 1 shown, the schematic diagram of the pre-tracking method based on the EM algorithm for successive elimination detection in this embodiment includes the following steps: Step 1. Initialization: Set the target speed range, accumulate the number of frames, frame detection ratio threshold, average error threshold, and distance threshold; Step 2. Pre-detection of current frame data: Perform CFAR low-threshold pre-detection on the data of the current time frame of the same azimuth beam; Step 3. Current frame trace condensation: Use amplitude weighting to estimate the centroid of the traces continuously distributed in the range direction to obtain the condensed trace zi; Step 4. Sliding window multi-frame traces: Form multi-frame traces z = [zi - F + 1, zi - F + 2, …, zi]T, where F is the number of frames in the sliding window; Step 5. Initialization of motion parameters based on Halton sequence Use the Halton sequence to generate multiple sets of linear motion parameters X, calculate the error and average error between the multi-frame traces and the theoretical points of each set of parameters, and take the motion parameters with the smallest average error as the initial target motion parameters a corresponding to the current multi-frame traces; Step 6. Weighted iterative single-target motion parameter estimation based on EM algorithm According to the error and average error between the trace and the theoretical point Ca, calculate the weight of this point according to the Gaussian distribution, generate the weight matrix W = diag(w), and then use Recalculate the motion parameters until the parameter error between the previous and the current time is less than the threshold, and output the corresponding target motion parameters a; Step 7. Multi-target parameter estimation by successively excluding used traces Select the point with the closest distance to the theoretical trace and less than the error threshold within each frame time according to the parameter a. If the number of points is greater than the ratio threshold, form new target motion parameters; otherwise, do not form new target motion parameters; Exclude the selected points from the traces. Compare the remaining number of points with the total number of frames. If it is greater than the ratio threshold, repeat Step 5; otherwise, output the final multi-target parameters of the current time frame; Step 8. Repeat Step 2 - Step 7 until the data frame ends.

[0025] Example 1: 1. Data scenario Set 2 target uniform linear motion targets, and the motion equations are 1) r1 = 200 + 200 * t, 2) r2 = 1000 - 140 * t, that is, the initial distances are 200m and 1000m respectively, and the initial radial velocities are 200m / s and -140m / s respectively. After 6 frames, the frame interval is 1 second, that is, t = 0, 1, 2, 3, 4, 5 seconds, see Figure 2 . The signal-to-noise ratio of the echo of each target per frame is 9dB, and the echo amplitudes of each frame with a total of 20 range cells (a total of 1500m, range resolution is 75m) are shown in the following table:

[0026]

[0027] The corresponding echo amplitude is shown in Figure 3 。

[0028] 2. Initialization: Set the CFAR detection threshold to 2.5119 (8 dB), the protection unit to 2, the unilateral window length to 12; the frame ratio threshold to 0.55; the target speed range to [-600, 600] m / s, and the distance range to [-1500, 1500] 3. Frame data detection: Perform bilateral CFAR detection on each frame of data, and the background is obtained as shown in the following table:

[0029]

[0030] And the obtained plot distribution is shown in the following table (1 indicates that there is a target point at this position, and 0 indicates that there is no target point at this position):

[0031]

[0032] 4. Plot condensation: For the plots detected in each frame, perform condensation along the distance direction: If there are detection points in consecutive distance units, weight them according to the amplitude to form a point at the central position; if the plots are not consecutive, directly output their positions.

[0033] The plots after condensation of each frame of data are shown in Figure 4 。As can be seen from the figure, in the 3rd frame (time t = 2 s), the original number of plots is 3, and after condensation, there are only 2 plots. Among them, the two consecutive plots at distances of 525 m and 600 m are condensed into one plot, and the position is 561.1. Comparing with the set target moving points Figure 2 and the condensed points after detection Figure 4 , it can be seen that there are a total of 12 plots in 6 frames of data, 5 plots are detected for target 1, and 4 plots are detected for target 2; in addition, there are 3 false alarm points. The obtained plots are z = [225, 375, 525, 825, 561.1273, 750, 600, 825, 1350, 975, 300, 1500 ]T.

[0034] 5. Initialization of motion parameters based on the Halton sequence Take the trace time \(t\) and a sequence of all 1s with the same length as the time to form \(C = [t, l]=[0\ 1\ 1\ 1\ 2\ 2\ 3\ 3\ 3\ 4\ 5\ 5; 1\ 1\ 1\ 1\ 1\ 1\ 1\ 1\ 1\ 1\ 1\ 1]^T\). Then use the function \(P = haltonset(2)\) provided by Matlab 2019 to generate a 2D Halton object, and use \(xx = net(P, 100)\) to generate 100 groups of 2D sequences. Use \(ii = randperm(100, 10)\) to take the indices of the first 10 groups of sequences, and obtain the parameters of the first 10 groups \(x = xx(ii, :)\). See the following table:

[0035] Then use \(v = v0+vlmt*x(:, 1); r = r0+rlmt*x(:, 2)\) to generate the motion parameters of velocity and initial distance for 10 groups as the motion parameter \(X\). See the following table:

[0036] Form 10 groups of theoretical traces through \(Y = CX\), compare them with 12 point condensation traces, calculate the average error, and obtain that the group of parameters with the smallest average error is the 6th group, and the parameters are \(a = [150, 833.3333]^T\), and the average error = 332921.549.

[0037] 6. Weighted Iterative Single Target Motion Parameter Estimation Based on EM Algorithm According to the parameter \(a\), the weight \(w\) obtained by calculating the distance difference between the theoretical trace and the actual trace is

[0038] Then press \(a = Recalculate, and get \(a = [204.037, 472.5921]^T\). The modulus of the difference between \(a\) before and after the first time is 364.766; continue to calculate the theoretical trace and the actual weight obtained, and recalculate \(a = [221.373, 307.4429]^T\). The modulus of the difference between \(a\) before and after the second time is 166.0566; the third time obtained \(a = [224.0392, 251.8065]^T\), and the modulus of the difference between \(a\) before and after is 55.7003; the fourth time obtained \(a = [226.0059, 219.3082]^T\), and the modulus of the difference between \(a\) before and after is 32.4779; it can be seen that the parameters gradually approach the motion parameters \([200, 200]^T\) of the first target. Until the 25th time \(a = [200, 175]^T\), the modulus of the difference between \(a\) before and after is \(1.536e - 5\), which is lower than the modulus difference threshold of 0.01, and the loop is exited and \(a\) is output. Compared with the set parameters, the velocity component is the same, and the distance component only differs by 1 / 3 of a distance resolution unit; output this parameter as the parameter of the first target found.

[0039] 7. Multi-target parameter estimation based on successive elimination of used traces According to the first target parameter, find the point closest to the theoretical point in each frame, and remove the points with an error greater than the resolution cell of 75 m, obtaining 5 traces [225, 375, 561.1273, 825, 975] corresponding to the target of this parameter, and the corresponding times are [0, 1, 2, 3, 4]. Since 5 / 6 = 0.8333 is greater than the set frame ratio threshold of 0.55, the found traces are valid and are included in the new target parameter. After removing the points found for the first target from the original traces, the remaining traces are z = [525, 825, 750, 600, 350, 300, 1500]T, and the corresponding frame times are t = [1, 1, 2, 3, 3, 5, 5]T. Use the remaining traces to continue with step 6, obtaining the parameters of another target as [-150, 1050], which are close to the set parameters [-140, 1000]T of the second target. The corresponding traces found are [825, 750, 600, 300]T, and the corresponding frame times are [1, 2, 3, 5]. After removing the traces of the second target, only 3 traces remain. Since 3 / 6 = 0.5 is less than the frame detection threshold, the target parameters have been found, the loop is exited, and the final 2 target parameters and their respective corresponding traces are output.

[0040] Although the present invention has been disclosed above in a preferred embodiment, the embodiments and the drawings are not used to limit the present invention. Any person skilled in this art can make various changes or modifications without departing from the spirit and scope of the present invention, but still within the protection scope of the present invention. Therefore, the protection scope of the present invention should be defined by the protection scope of the claims of this application.

Claims

1. A sequential elimination detection-before-tracking method based on the EM algorithm, characterized in that, It includes the following steps: Step 1, Initialization: Set the target speed range, accumulate the number of frames, frame detection ratio threshold, average error threshold, and distance threshold. Step 2, Pre-detection of current frame data: Perform CFAR low-threshold pre-detection on the data of the current time frame of the same azimuth beam. Step 3, Current frame trace condensation: Use amplitude weighting to estimate the centroid of the traces continuously distributed in the range direction to obtain the condensed trace zi. Step 4, Sliding window multi-frame traces: Form multi-frame traces z = [zi - F + 1, zi - F + 2, …, zi]T, where F is the number of frames in the sliding window. Step 5, Initialization of motion parameters based on the Halton sequence: Use the Halton sequence to generate multiple groups of linear motion parameters X, calculate the error and average error between the multi-frame traces and the theoretical points of each group of parameters, and take the motion parameters with the smallest average error as the initial target motion parameters a corresponding to the current multi-frame traces. Step 6, Weighted iterative single-target motion parameter estimation based on the EM algorithm: According to the error and average error between the track and the theoretical point, calculate the weight of this point according to the Gaussian distribution, generate the weight matrix W = diag(w), and then use Recalculate the motion parameters until the parameter error between the previous and the next time is less than the threshold, and output the corresponding target motion parameter a; Step 7, Multi-target parameter estimation by sequentially excluding used traces: According to the parameter a, pick out the points closest to the theoretical trace and less than the error threshold within each frame time. If the number of points is greater than the ratio threshold, form new target motion parameters, otherwise do not form new target motion parameters; exclude the picked points from the traces. Compare the remaining number of points with the total number of frames. If it is greater than the ratio threshold, repeat Step 5, otherwise output the final multi-target parameters of the current time frame. Step 8, Repeat Step 2 - Step 7 until the data frame ends.

2. The sequential elimination detection-before-tracking method based on the EM algorithm according to claim 1, wherein, The specific content of Step 2 is as follows: For each range cell, use g cells as guard cells, average the statistical noise background of M points on both the left and right. If the amplitude of the current cell is higher than the background by the threshold, it is judged that there is a target, and the position and amplitude are output, otherwise it is judged that there is no target.

3. A sequential rejection tracking-before-detection method based on the EM algorithm according to claim 1, characterized in that, The specific content of Step 3 is as follows: For the traces detected in each frame, perform condensation along the range direction: If there are detected points in consecutive range cells, weight them by amplitude to form a point at the center position; if the traces are not continuous, directly output their positions.

4. The sequential rejection detection-before-tracking method based on the EM algorithm according to claim 1, wherein The specific content of Step 5 includes: Step 5.1: Generate K groups of Halton sequences with dimension 2; and generate the theoretical points of multiple groups of parameters according to Y = CX. Y is the multiple groups of target theoretical points formed by multiple groups of motion parameters, C is a matrix constructed by the vector t formed by the frame time of the traces and the all-ones sequence 1, C = [t, 1], t is the frame time sequence of the traces, and 1 is an all-ones sequence of the same length as t, with the length being the number of traces N. Step 5.2: Generate multiple groups of target motion parameters X = [v, r]T according to the target distance, speed range, and Halton sequence, where v = v0 + vlmt * halton(k, 1), r = r0 + rlmt * halton(k, 2), k = 1, 2, …, K. Step 5.3: For each set of parameters Xi, calculate the error d between the multi-frame traces and the theoretical trace CXi j =Y j -z and the average error ,Y j is the j-th column of Y; Step 5.4: Take the parameters with the minimum average error as the initial target motion parameters corresponding to the current trace , , is the k-th column of the K sets of linear motion parameters X.

5. A sequential rejection detection-before-tracking method based on the EM algorithm according to claim 5, characterized in that: The specific content of Step 7 includes: Step 7.1: According to the parameter a output in Step 6, pick out the points closest to the theoretical trace and less than the error threshold within each frame time. Step 7.2: If the number of picked points is greater than the ratio threshold, form new target motion parameters, otherwise do not form new target motion parameters. Step 7.3: Remove the selected points from the point tracks and replace z with the remaining point tracks; Step 7.4: Compare the number of remaining point tracks z with the total number of frames. If it is greater than the ratio threshold, go to Step 6 for weighted iterative single-target motion parameter estimation based on the EM algorithm; otherwise, output the final multi-target parameters. According to a successive elimination pre-detection tracking method based on the EM algorithm as claimed in claim 1, wherein said Step 6 specifically includes: Step 6.1: Calculate the error d = z - Ca between the trace z and the theoretical trace Ca of the parameter a and the average error , where N is the total number of traces z; Step 6.2: Generate the weights of the traces according to the Gaussian distribution , where n = 1, 2, …, N, is the n-th component of d, the weight matrix W = diag(w), and w = [w1, w2, …, wN]T; Step 6.3: Let and recalculate the motion parameters . Step 6.4: Calculate the theoretical target moving point \(z_t = Ca\) under the new parameters, and recalculate the error \(d\) between the trace \(z\) and the theoretical moving point \(z_t\), the average error and the weight ; Step 6.5: Calculate the difference e between the parameters of the two times before and after as follows ; Step 6.6: If e is greater than the threshold, continue with Step 6.1; otherwise, exit and output a.

6. The sequential rejection detection-before-tracking method based on the EM algorithm according to claim 5, characterized in that: Said Step 7 specifically includes: Step 7.1: Select the points closest to the theoretical point tracks and less than the error threshold within each frame of time according to the parameter a output by Step 6; Step 7.2: If the number of selected points is greater than the ratio threshold, form new target motion parameters; otherwise, do not form new target motion parameters; Step 7.3: Remove the selected points from the point tracks and replace z with the remaining point tracks; Step 7.4: Compare the number of remaining point tracks z with the total number of frames. If it is greater than the ratio threshold, go to Step 6 for weighted iterative single-target motion parameter estimation based on the EM algorithm; otherwise, output the final multi-target parameters.

Images