A through-wall radar multi-motion target robust tracking method based on multi-domain data association and multi-state trajectory management

Abstract

This invention discloses a robust tracking method for multiple moving targets using through-wall radar based on multi-domain data association and multi-state trajectory management. First, radar measurement results are obtained using a constant false alarm rate (CFAR) detection method. Second, the Euclidean distance between the measurements and the trajectories, the kernel correlation filter response, and the trajectory confidence cost are calculated to construct a minimum-cost maximum flow problem for data association. Subsequently, a 5-state trajectory management mechanism is established. Based on the target state and association results, new trajectories are generated, regular and unmatched trajectories are updated, and cross-trajectories are generated, updated, and separated. Specifically, in the cross-trajectory separation process, Bézier curves are used to fit the cross-trajectories, and the target state is estimated based on the Bézier curves to predict the target's position in the next frame. The re-matching of cross-trajectories and measurements is achieved based on the distance between the predicted position and the measurements.

Description

Technical Field

[0001] This invention belongs to the field of through-wall radar imaging detection, and relates to a robust tracking method for multiple moving targets using through-wall radar based on multi-domain data association and multi-state trajectory management. Background Technology

[0002] Indoor moving target tracking has become a research hotspot in recent years. However, existing indoor detection methods, such as indoor Wi-Fi positioning, are limited in non-cooperative environments. Through-wall radar, by emitting L / S band electromagnetic waves that penetrate walls, can detect indoor moving targets undetected, making it an effective means of tracking indoor moving targets in non-cooperative environments.

[0003] The classic radar tracking method, the Kalman filter, is based on Bayesian filtering theory and achieves tracking by assuming the target's motion state is known. This method and its improvements are widely used in moving target tracking in through-wall radar. One study first achieves two-dimensional positioning through geometric triangulation, then uses the Kalman filter for target tracking. Another study uses both the Kalman filter and particle filter to track the target simultaneously, improving tracking accuracy. The STFT–NGRC method can improve the accuracy of instantaneous target frequency estimation and is used in the target prediction process to improve the tracking performance of targets behind walls. In addition, the Kalman filter is combined with multipath suppression methods to achieve moving target tracking in indoor multipath environments. Another study analyzes the tracking performance of the Kalman filter, extended Kalman filter, and H-infinite a posteriori filter for moving targets behind walls, showing that these filters have limited performance in multi-target tracking. Another study uses interactive multi-model tracking to achieve multi-moving target tracking. While this type of motion model-based tracking method can effectively track moving targets with specific trajectories, human targets are maneuverable, making it difficult to accurately fit them using motion models. Ding et al. applied deep learning methods to through-wall radar and obtained good target tracking results, but such methods are only applicable to Doppler through-wall radar.

[0004] Ultra-wideband MIMO through-wall radar often uses imaging to detect targets behind walls; therefore, image-based tracking methods can be used for target tracking behind walls. Liu et al. were the first to introduce the mean-shift algorithm into through-wall radar for target tracking. The mean-shift algorithm first quantizes the radar image amplitude, then uses the histogram distribution characteristics of the target position to construct a template for matching the target position in the next frame of the radar image. Based on the shape characteristics of the target position in radar imaging, Li et al. proposed a scale-adaptive mean-shift algorithm, which replaces the rectangular region of interest with an ellipse, making it more suitable for radar images and effectively improving the tracking performance of moving targets in through-wall radar. Li et al. further combined the mean-shift algorithm with a Kalman filter to improve the tracking performance of moving targets on specific trajectories. Yao et al. proposed a trajectory association method for target identification (ID) continuity in distributed TWIR. This method achieves association between different trajectories even in the presence of occlusions in the environment. The kernel correlation filtering method uses the directional gradient histogram of the image instead of the traditional amplitude histogram. In our previous work, we introduced kernel correlation filtering into through-wall radar and modified the region of interest (ROI) to vary with the target position based on the target shape in the radar image, achieving good moving target tracking performance. However, due to the non-uniformity of the wall, some frames in the through-wall radar image become defocused or even fail to form a correct image, causing the radar image-based tracking algorithm to fail. Furthermore, when the trajectories of multiple targets intersect or separate, tracking methods relying on radar images will produce erroneous tracking, leading to errors in the fusion of multiple trajectories. Summary of the Invention

[0005] In view of this, this invention provides a robust multi-moving target tracking method for through-wall radar based on multi-domain data association and multi-state trajectory management, applicable to through-wall radar imaging detection. First, radar measurement results are obtained using a constant false alarm rate (CFAR) detection method. Second, the Euclidean distance between the measurement and the trajectory, the kernel correlation filter response, and the trajectory confidence cost are calculated to construct a minimum-cost maximum flow problem for data association. Subsequently, a 5-state trajectory management mechanism is established. Based on the target state and association results, new trajectory generation, regular trajectory and unmatched trajectory updates, and cross trajectory generation, updating, and separation are implemented. Specifically, in the cross trajectory separation process, Bézier curves are used to fit the cross trajectories, and the target state is estimated based on the Bézier curves to predict the target's position in the next frame. Re-matching of the cross and measurement is achieved based on the distance between the predicted position and the measurement. This method can still correctly track multiple targets even in situations such as radar image defocusing and trajectory intersection, making it a general and robust image-domain multi-target tracking method.

[0006] To achieve the above objectives, the robust tracking method for multiple moving targets using through-wall radar based on multi-domain data association and multi-state trajectory management of the present invention includes the following steps: Step 1: Trajectory initialization, generating 4 types of trajectory information and 5 states; The target trajectory is assigned four types of information: identification information, tracking quality information, image information, and status information. There are five trajectory states: "initial" indicates the trajectory is in the initialization state; "UB" indicates the trajectory is in the pending start stage, where the trajectory has begun tracking the target but its lifetime is low; "RUN" indicates the trajectory is in normal tracking state, and when the lifetime of the "UB" segment is greater than 0.5, the trajectory state changes to "RUN" (meaning that if the trajectory correctly matches the focused target for 15 consecutive frames (equivalent to 0.14 seconds of human movement for the radar), it is considered a correct moving target trajectory); "IS" indicates the trajectory intersects with other trajectories, in which case `relation_other` stores the ID of the intersecting trajectory; and "END" indicates the trajectory has stopped.

[0007] Step 2: Generate radar imaging results using multi-channel radar echoes; The compression result of the received echo pulse from the nth channel of the radar in the p-th period can be expressed as: In the formula, l represents the two-way distance at which the radar detects the target. The wall echo pulse compression result of the nth channel is represented, which does not change with the detection period and can be removed by mean cancellation. Tpr represents the pulse duration of a single period, B represents the signal bandwidth, f0 represents the initial frequency of the signal, and c represents the speed of electromagnetic wave propagation in air. This represents the distance between the j-th human target and the n-th radar channel during the p-th period.

[0008] Back projection algorithms are commonly used to generate radar images. They divide the imaging region into different grids and utilize phase compensation to achieve coherent energy accumulation at the target location. For the imaging grid points... The radar imaging results can be expressed as . In the formula, This represents the nth channel and the imaging grid points. The distance between them. Furthermore, high-quality radar images can be achieved using CF-type algorithms.

[0009] Step 3: Divide the training units and obtain the measurement points using the constant false alarm rate (CFAR) detection method; For the p-th frame of the radar image Its training units can be represented as a set in, , and These represent the major axis, minor axis, and rotation angle of the ellipse corresponding to position C of the imaging grid. It is a constant. Therefore, the image after CFAR detection can be represented as In the formula, α is the threshold factor for CFAR detection, which is related to the false alarm rate setting and the number of training units. By performing connected component detection on the detected image and determining its centroid position, the measurement result of the p-th frame can be obtained, which can be expressed as: .

[0010] Step 4: Calculate the multi-domain cost between the existing trajectory and the measurement points, and then establish the minimum cost maximum flow problem to achieve data association by solving this problem; Assume that the total number of existing trajectories in the p-th radar cycle is This includes trajectories in all trajectory states. This module is used to implement... Matching one trajectory with Q measurements.

[0011] First, the cost between the q-th measurement and the t-th trajectory can be calculated as follows: In the formula, , and There are three weighting coefficients, and they satisfy... . This represents the distance cost between the q-th measurement and the t-th trajectory, used to measure the distance of a human target's movement within a single radar cycle. Its calculation method is as follows: In the formula, This represents the tracking result of the t-th trajectory in the (p-1)-th radar cycle.

[0012] This represents the kernel correlation filter (KCF) response cost between the q-th measurement and the t-th trajectory, used for the target range after extraction and rotation of the q-th measurement position. The correlation between the t-th trajectory and the radar image can be calculated as follows: In the formula, The last_area represents the t-th trajectory. The filter_para represents the t-th trajectory.

[0013] Let represent the trajectory lifetime of the t-th trajectory, which represents the confidence level for that trajectory, and can be expressed as: The data association problem is modeled as a minimum-cost maximum-flow problem. Trajectory nodes and measurement nodes are added between the source and sink. The trajectory nodes represent existing... A trajectory. Measurement nodes represent the existing Q measurements. The capacity and cost between the source point and trajectory nodes depend on the current state of the trajectory. For a trajectory in state 'UB', the capacity of the path between the source point and this trajectory node is 1, and the cost is 1e-6. For a trajectory in state 'RUN', the capacity of the path between the source point and this trajectory node is 1, and the cost is 2e-6. For a trajectory in state 'IS', the capacity of the path between the source point and this trajectory node is the number of trajectories intersecting with this trajectory, and the cost is 2e-6. For a trajectory in state 'END', the capacity of the path between the source point and this trajectory node is 0. The capacity between the source point and trajectory nodes represents the maximum number of measurements that each trajectory can be associated with. Cost allocation ensures that measurement values ​​are preferentially allocated to trajectories in motion rather than those tentatively starting.

[0014] The capacity between trajectory nodes and measurement nodes is 1, meaning that each trajectory can be matched with each measurement at most once. The cost is the price between the corresponding trajectory and measurement.

[0015] The cost between the measurement node and the sink is 0. Capacity represents the maximum number of trajectories a target can be associated with, and this value affects how many trajectories will intersect. Considering radar resolution and human body size, this value is usually taken as 3.

[0016] Data association can be achieved by solving the minimum-cost maximum-flow problem described above, ultimately yielding the data association matrix DA. When the q-th measurement is associated with the t-th trajectory... When the q-th measurement is associated with the t-th trajectory, .

[0017] Step 5: Establish a new tracking trajectory using measurements that are not associated with any trajectory, and initialize the trajectory parameters; When the qth measurement of the p-th frame satisfies the following formula Then, a new trajectory needs to be generated starting from the q-th measurement. That is, the q-th measurement... The initial value of the trajectory is .

[0018] The target region corresponding to the target can be extracted as follows: Among them, set Let the elliptical range corresponding to the new trajectory in frame p be represented as follows: In the formula, , and These represent the major axis, minor axis, and rotation angle of the ellipse corresponding to the measurement position, respectively. Given... In this case, the target region can be rotated using bilinear interpolation, expressed as: A correlation filter can be constructed, where the filter parameters can be expressed as... In the formula, and These represent the forward and inverse Fourier transforms, respectively. This represents a two-dimensional Gaussian function, whose dimensions are... Maintain consistency This indicates the conjugate operation, and σ represents the regularization parameter. The operator is defined as being computable as In the formula, β is a constant representing the bandwidth of the Gaussian space.

[0019] Therefore, the last_area corresponding to this trajectory is obtained. and filter_para In addition, the ID of this trajectory is... Both begin and end are p, total is 1, SAN is 1, CLN is 0, and now_state is 'UB'.

[0020] Step 6: For the trajectory corresponding to the associated measurement, update all trajectory information based on the associated measurement; For the t-th trajectory, when its now_state is not 'IS', and the trajectory and its matched q-th measurement satisfy the following formula, This trajectory is a normal trajectory. For normal trajectories, the results obtained using the KCF method can be used as prediction results.

[0021] Based on the target position in the previous frame The target region can be extracted. Then, the rotated target area can be obtained. To prevent discrepancies in target region size due to varying target locations, the extracted target region needs to be reconstructed. In the formula, Indicates the target area in the first frame The dimensionality. Using the reconstructed target region, the relevant response can be calculated, represented as... Maximum value in the relevant response The position represents the estimated target position in frame p. Finally, using this location, the positioning result can be obtained. and the target area corresponding to this positioning result. and updated filter parameters Based on the results, the regular trajectory can be updated, where the total number of the trajectory is increased by 1, the end is set to p, the SAN number is increased by 1, the CLN is set to 0, and the new last_area can be represented as... The new filter_para can be represented as The radar positioning result of this frame can be represented as follows: .

[0022] Step 7: For trajectories not associated with any measurements, update the trajectories based on trajectory information and radar images; For the t-th trajectory, when its now_state is not 'IS' and it satisfies the following formula: Then the trajectory did not match any measurement points. The trajectory is considered a match if now_state is 'UB' or CLN is greater than or equal to... If this happens, its now_state is changed to 'END'. In addition, the total number of tracks is reduced. The tracking results suggest that The result was discarded because it was invalid.

[0023] If the trajectory now_state is not 'UB' and CLN is less than 'UB' If the trajectory is not found, KCF is used to continue tracking it, and the KCF method localization result is obtained. and the maximum value of the KCF response Based on the results, unassociated trajectories can be updated, where the end of the trajectory is set to p, CLN is incremented by 1, last_area and filter_para remain unchanged, and now_state is updated. The radar tracking result for this frame can be represented as... .

[0024] Step 8: Generate cross trajectories based on the association results and update the state of the generated cross trajectories; update or separate the cross trajectories based on the measurements; in the cross trajectory separation, first use Bézier curves to fit each cross trajectory, then use the motion state to predict the trajectory, use the prediction results to re-associate with the associated measurements, and finally use the re-association results to separate the cross trajectories. The now_state of an intersecting trajectory can be represented as 'IS'.

[0025] One of the conditions for trajectory intersection is that a certain measurement q matches more than one trajectory with now_state 'RUN'. When a certain measurement matches multiple trajectories, the trajectory with now_state 'UB' needs to be changed to 'END'.

[0026] Assume the number of intersecting trajectories is When a trajectory intersection occurs, the intersection... The now_state of all trajectories is converted to "IS", and the IDs of the remaining trajectories are placed in the relation_other of that trajectory. Following the normal trajectory update steps, the measurement q is considered as the measurement corresponding to the intersecting trajectory and updated accordingly. A trajectory.

[0027] For trajectories that have already intersected, the number of measurements that need to be matched to these trajectories needs to be determined. Make a judgment. Perform corresponding operations based on different measurement quantities. When the measurement quantity... At that time, different measurement trajectories are updated based on the criteria for unmatched measurements. When the number of measurements is... At that time, the measurement was considered as follows: For each trajectory, a measurement is taken, and the regular trajectory update steps are performed. The difference is that only the trajectory with the longest lifespan and the largest total number of SANs is incremented by 1, while the CLN of all other trajectories is incremented by 1. Then, if the CLN of a certain trajectory is greater than the gateCLN, the now_state of that trajectory is changed to "END".

[0028] When the cross trajectory matches the measurement quantity In such cases, it is necessary to separate multiple trajectories. The method presented in this paper utilizes the target's motion state to achieve target separation. However, since human targets are maneuvering targets with highly variable motion states, they cannot be described by simple motion models, thus motion model-based tracking methods are not applicable. Literature has demonstrated that, under non-external force conditions, the acceleration curve of a human target's motion is smooth; therefore, the target's acceleration at the next moment can be estimated using existing trajectories, thereby obtaining the target's position estimation result.

[0029] Bézier curves can be fitted to smooth curves using control points. This paper uses Bézier curve fitting to fit the trajectory of a human target, estimates the human acceleration based on the fitted Bézier curve, and finally obtains the target tracking prediction result.

[0030] For a Bézier curve of order Z, it can be represented as: In the formula, Let Z represent the z-th control point. For a Bézier curve of order Z, the number of control points is Z+1. The Bernstein function corresponding to the control point can be expressed as: For the There are intersecting trajectories, and their total trajectory length is The overall trajectory can be represented as Therefore, we can obtain By solving (31) using the least squares method, the control points corresponding to the intersection trajectory can be obtained. Furthermore, the first Bézier curve fitting can be obtained. By analyzing the intersecting trajectories and then using trajectory difference, the corresponding velocity and acceleration curves can be obtained. Since the acceleration curve of a human target's motion is smooth, the position prediction result of the trajectory can be obtained. .

[0031] By constructing a new cost function, cross trajectories can be achieved. Matched measurements The cost of a quadratic matching between the q'-th measurement and the t'-th intersection trajectory can be expressed as: Similarly, the minimum-cost maximum flow problem can be used to rematch the intersecting trajectories with the measurements, such as... Figure 4As shown. The capacity between the source node and the intersection trajectory node is 1, and the cost is 0, indicating that each intersection trajectory requires a re-matching process; the capacity between the intersection trajectory node and the already matched target node is 1, and the cost is... The capacity between the matched target node and the sink node is... The cost is 0. Similarly, when the q'-th measurement is associated with the t'-th trajectory, When the q'th measurement is not associated with the 'tth trajectory, .

[0032] Based on the reassociation results, when a trajectory meets the condition that all targets reassociated with it are associated only with this trajectory, its now_state becomes "RUN," and it is removed from the "relation_other" of other intersecting trajectories. Otherwise, the trajectory's now_state remains "IS." Subsequently, the average of all measurements reassociated with this trajectory is used as its new measurement, and the regular trajectory update steps are performed to achieve the separation of intersecting trajectories.

[0033] Beneficial effects: 1. This invention proposes a multi-domain data association method: First, the association cost between measurements and trajectories is calculated using Euclidean distance, kernel correlation filter response, and trajectory information; then, a minimum cost maximum flow model is constructed to realize the association between multiple measurements and multiple trajectories.

[0034] 2. This invention introduces a multi-state trajectory management strategy: First, a five-state trajectory management mechanism is established; then, new trajectories are generated, regular trajectories are updated, and unassociated trajectories are processed based on the association results. Specifically, when an intersecting trajectory is associated with multiple measurements, Bézier curves are used to fit the intersecting trajectory for prediction, and then the trajectory is reassociated with the measurements to achieve trajectory separation.

[0035] Compared with existing methods, this method can still correctly track multiple targets even when radar images are defocused or have intersecting trajectories. Attached Figure Description

[0036] Figure 1 Method flowchart; Figure 2 1. A scene diagram of through-wall radar detection; Figure 3 Solving the minimum cost maximum flow problem for data association; Figure 4 Solving the minimum cost maximum flow problem and re-associating data; Figure 5 Simulation scenarios; Figure 6 Simulation results; Figure 7 1. Schematic diagram of the measured radar; Figure 8 Actual measurement results. Detailed Implementation

[0037] The present invention discloses a robust tracking method for multiple moving targets using through-wall radar based on multi-domain data association and multi-state trajectory management. The processing flowchart is as follows: Figure 1 As shown, the specific steps are as follows: Step 1: Trajectory initialization, generating 4 types of trajectory information and 5 states; The target trajectory is assigned four types of information: identification information, tracking quality information, image information, and status information. The specific contents of each type of information are shown in Table 1.

[0038] Table 1 Tracking trajectory information

[0039] In Table 1, the method for calculating the trajectory lifetime is as follows: There are five track states. "initial" indicates the track is in the initialization state; "UB" indicates the track is in the pending start stage, where it has begun tracking the target but has a short lifespan; "RUN" indicates the track is in normal tracking state. When the lifespan of the "UB" segment is greater than 0.5, the track state changes to "RUN" (meaning that when the track's lifespan is greater than 0.5 for 15 consecutive frames...). For radar, if a human movement (within 0.14 seconds) correctly matches the focused target, then the trajectory is considered a correct moving target trajectory. "IS" indicates that the trajectory intersects with other trajectories; in this case, `relation_other` stores the ID of the intersecting trajectory. "END" indicates that the trajectory has stopped.

[0040] Step 2: Generate radar imaging results using multi-channel radar echoes; An ultrawideband MIMO through-wall radar is placed close to a wall to detect multiple moving targets behind the wall, as shown in the scene diagram. Figure 2 As shown. An ultra-wideband radar consists of N channels. It detects the scene behind a wall by emitting ultra-wideband signals (such as linear frequency modulated continuous wave signals, stepped frequency signals, etc.). Assume the radar emits P cycles of signals to detect J human targets behind the wall. Due to human movement, the human's position changes in different radar cycles. Therefore, the pulse compression result of the echo received by the nth channel of the radar in the p-th cycle can be expressed as... In the formula, l represents the two-way distance at which the radar detects the target. The wall echo pulse compression result of the nth channel is represented, which does not change with the detection period and can be removed by mean cancellation. Tpr represents the pulse duration of a single period, B represents the signal bandwidth, f0 represents the initial frequency of the signal, and c represents the speed of electromagnetic wave propagation in air. This represents the distance between the j-th human target and the n-th radar channel during the p-th period.

[0041] Back projection algorithms are commonly used to generate radar images. They divide the imaging region into different grids and utilize phase compensation to achieve coherent energy accumulation at the target location. For the imaging grid points... The radar imaging results can be expressed as . In the formula, This represents the nth channel and the imaging grid points. The distance between them. Furthermore, high-quality radar images can be achieved using CF-type algorithms.

[0042] Step 3: Divide the training units and obtain the measurement points using the constant false alarm rate (CFAR) detection method; CFAR detection is a commonly used detector in radar signal processing. Its principle is to estimate the noise amplitude by measuring the amplitude within the training unit and then use this as a threshold to achieve target detection. Common two-dimensional CFAR detection in radar images often uses rectangular calibration of the guard and training units. In the algorithm proposed in this paper, the calibration range of the guard and training units is replaced with an ellipse based on the shape of the target in the radar image.

[0043] For the p-th frame of the radar image Its training units can be represented as a set in, , and These represent the major axis, minor axis, and rotation angle of the ellipse corresponding to position C of the imaging grid. It is a constant. Therefore, the image after CFAR detection can be represented as In the formula, α is the threshold factor for CFAR detection, which is related to the false alarm rate setting and the number of training units. By performing connected component detection on the detected image and determining its centroid position, the measurement result of the p-th frame can be obtained, which can be expressed as: .

[0044] Step 4: Calculate the multi-domain cost between the existing trajectory and the measurement points, and then establish the minimum cost maximum flow problem to achieve data association by solving this problem; Assume that the total number of existing trajectories in the p-th radar cycle is This includes trajectories in all trajectory states. This module is used to implement... Matching one trajectory with Q measurements.

[0045] First, the cost between the q-th measurement and the t-th trajectory can be calculated as follows: In the formula, , and There are three weighting coefficients, and they satisfy... . This represents the distance cost between the q-th measurement and the t-th trajectory, used to measure the distance of a human target's movement within a single radar cycle. Its calculation method is as follows: In the formula, This represents the tracking result of the t-th trajectory in the (p-1)-th radar cycle.

[0046] This represents the kernel correlation filter (KCF) response cost between the q-th measurement and the t-th trajectory, used for the target range after extraction and rotation of the q-th measurement position. The correlation between the t-th trajectory and the radar image can be calculated as follows: In the formula, The last_area represents the t-th trajectory. The filter_para represents the t-th trajectory.

[0047] Let represent the trajectory lifetime of the t-th trajectory, which represents the confidence level for that trajectory, and can be expressed as: The minimum-cost maximum flow problem is a classic problem in graph theory, often applied in deep learning and network traffic optimization. Essentially, it involves finding the maximum flow between the source and sink nodes, while also considering the cost between different nodes. For example... Figure 3 As shown, the data association problem can be modeled as a minimum-cost maximum-flow problem. Trajectory nodes and measurement nodes are added between the source and sink. The trajectory nodes represent existing... A trajectory. Measurement nodes represent the existing Q measurements. The capacity and cost between the source point and trajectory nodes depend on the current state of the trajectory. For a trajectory in state 'UB', the capacity of the path between the source point and this trajectory node is 1, and the cost is 1e-6. For a trajectory in state 'RUN', the capacity of the path between the source point and this trajectory node is 1, and the cost is 2e-6. For a trajectory in state 'IS', the capacity of the path between the source point and this trajectory node is the number of trajectories intersecting with this trajectory, and the cost is 2e-6. For a trajectory in state 'END', the capacity of the path between the source point and this trajectory node is 0. The capacity between the source point and trajectory nodes represents the maximum number of measurements that each trajectory can be associated with. Cost allocation ensures that measurement values ​​are preferentially allocated to trajectories in motion rather than those tentatively starting.

[0048] The capacity between trajectory nodes and measurement nodes is 1, meaning that each trajectory can be matched with each measurement at most once. The cost is the price between the corresponding trajectory and measurement.

[0049] The cost between the measurement node and the sink is 0. Capacity represents the maximum number of trajectories a target can be associated with, and this value affects how many trajectories will intersect. Considering radar resolution and human body size, this value is usually taken as 3.

[0050] There are many solutions to the minimum-cost maximum flow problem, the simplest of which is the improved Dijkstra's method. This method is a greedy algorithm that finds a global optimum by finding a local optimum. Data association can be achieved by solving the minimum-cost maximum flow problem, ultimately yielding the data association matrix DA. When the q-th measurement is associated with the t-th trajectory, When the q-th measurement is associated with the t-th trajectory, .

[0051] Step 5: Establish a new tracking trajectory using measurements that are not associated with any trajectory, and initialize the trajectory parameters; When the qth measurement of the p-th frame satisfies the following formula Then, a new trajectory needs to be generated starting from the q-th measurement. That is, the q-th measurement... The initial value of the trajectory is .

[0052] The target region corresponding to the target can be extracted as follows: Among them, set Let the elliptical range corresponding to the new trajectory in frame p be represented as follows: In the formula, , and These represent the major axis, minor axis, and rotation angle of the ellipse corresponding to the measurement position, respectively. Given... In this case, the target region can be rotated using bilinear interpolation, expressed as: A correlation filter can be constructed, where the filter parameters can be expressed as... In the formula, and These represent the forward and inverse Fourier transforms, respectively. This represents a two-dimensional Gaussian function, whose dimensions are... Maintain consistency This indicates the conjugate operation, and σ represents the regularization parameter. The operator is defined as being computable as In the formula, β is a constant representing the bandwidth of the Gaussian space.

[0053] Therefore, the last_area corresponding to this trajectory is obtained. and filter_para In addition, the ID of this trajectory is... Both begin and end are p, total is 1, SAN is 1, CLN is 0, and now_state is 'UB'.

[0054] Step 6: For the trajectory corresponding to the associated measurement, update all trajectory information based on the associated measurement; For the t-th trajectory, when its now_state is not 'IS', and the trajectory and its matched q-th measurement satisfy the following formula, This trajectory is a normal trajectory. For normal trajectories, the results obtained using the KCF method can be used as prediction results.

[0055] Based on the target position in the previous frame The target region can be extracted. Then, the rotated target area can be obtained. To prevent discrepancies in target region size due to varying target locations, the extracted target region needs to be reconstructed. In the formula, Indicates the target area in the first frame The dimensionality. Using the reconstructed target region, the relevant response can be calculated, represented as... Maximum value in the relevant response The position represents the estimated target position in frame p. Finally, using this location, the positioning result can be obtained. and the target area corresponding to this positioning result. and updated filter parameters Based on the results, the regular trajectory can be updated, where the total number of the trajectory is increased by 1, the end is set to p, the SAN number is increased by 1, the CLN is set to 0, and the new last_area can be represented as... The new filter_para can be represented as The radar positioning result of this frame can be represented as follows: .

[0056] Step 7: For trajectories not associated with any measurements, update the trajectories based on trajectory information and radar images; For the t-th trajectory, when its now_state is not 'IS' and it satisfies the following formula: Then the trajectory did not match any measurement points. The trajectory is considered a match if now_state is 'UB' or CLN is greater than or equal to... If this happens, its now_state is changed to 'END'. In addition, the total number of tracks is reduced. The tracking results suggest that The result was discarded because it was invalid.

[0057] If the trajectory now_state is not 'UB' and CLN is less than 'UB' If the trajectory is not found, KCF is used to continue tracking it, and the KCF method localization result is obtained. and the maximum value of the KCF response Based on the results, unassociated trajectories can be updated, where the end of the trajectory is set to p, CLN is incremented by 1, last_area and filter_para remain unchanged, and now_state is updated. The radar tracking result for this frame can be represented as... .

[0058] Step 8: Generate cross trajectories based on the association results and update the state of the generated cross trajectories; update or separate the cross trajectories based on the measurements; in the cross trajectory separation, first use Bézier curves to fit each cross trajectory, then use the motion state to predict the trajectory, use the prediction results to re-associate with the associated measurements, and finally use the re-association results to separate the cross trajectories. The now_state of an intersecting trajectory can be represented as 'IS'.

[0059] One of the conditions for trajectory intersection is that a certain measurement q matches more than one trajectory with now_state 'RUN'. When a certain measurement matches multiple trajectories, the trajectory with now_state 'UB' needs to be changed to 'END'.

[0060] Assume the number of intersecting trajectories is When a trajectory intersection occurs, the intersection... The now_state of all trajectories is converted to "IS", and the IDs of the remaining trajectories are placed in the relation_other of that trajectory. Following the normal trajectory update steps, the measurement q is considered as the measurement corresponding to the intersecting trajectory and updated accordingly. A trajectory.

[0061] For trajectories that have already intersected, the number of measurements that need to be matched to these trajectories needs to be determined. Make a judgment. Perform corresponding operations based on different measurement quantities. When the measurement quantity... At that time, different measurement trajectories are updated based on the criteria for unmatched measurements. When the number of measurements is... At that time, the measurement was considered as follows: For each trajectory, a measurement is taken, and the regular trajectory update steps are performed. The difference is that only the trajectory with the longest lifespan and the largest total number of SANs is incremented by 1, while the CLN of all other trajectories is incremented by 1. Then, if the CLN of a certain trajectory is greater than the gateCLN, the now_state of that trajectory is changed to "END".

[0062] When the cross trajectory matches the measurement quantity In such cases, it is necessary to separate multiple trajectories. The method presented in this paper utilizes the target's motion state to achieve target separation. However, since human targets are maneuvering targets with highly variable motion states, they cannot be described by simple motion models, thus motion model-based tracking methods are not applicable. Literature has demonstrated that, under non-external force conditions, the acceleration curve of a human target's motion is smooth; therefore, the target's acceleration at the next moment can be estimated using existing trajectories, thereby obtaining the target's position estimation result.

[0063] Bézier curves can be fitted to smooth curves using control points. This paper uses Bézier curve fitting to fit the trajectory of a human target, estimates the human acceleration based on the fitted Bézier curve, and finally obtains the target tracking prediction result.

[0064] For a Bézier curve of order Z, it can be represented as: In the formula, Let Z represent the z-th control point. For a Bézier curve of order Z, the number of control points is Z+1. The Bernstein function corresponding to the control point can be expressed as: For the There are intersecting trajectories, and their total trajectory length is The overall trajectory can be represented as Therefore, we can obtain By solving (31) using the least squares method, the control points corresponding to the intersection trajectory can be obtained. Furthermore, the first Bézier curve fitting can be obtained. By analyzing the intersecting trajectories and then using trajectory difference, the corresponding velocity and acceleration curves can be obtained. Since the acceleration curve of a human target's motion is smooth, the position prediction result of the trajectory can be obtained. .

[0065] By constructing a new cost function, cross trajectories can be achieved. Matched measurements The cost of a quadratic matching between the q'-th measurement and the t'-th intersection trajectory can be expressed as: Similarly, the minimum-cost maximum flow problem can be used to rematch the intersecting trajectories with the measurements, such as... Figure 4 As shown. The capacity between the source node and the intersection trajectory node is 1, and the cost is 0, indicating that each intersection trajectory requires a re-matching process; the capacity between the intersection trajectory node and the already matched target node is 1, and the cost is... The capacity between the matched target node and the sink node is... The cost is 0. Similarly, when the q'-th measurement is associated with the t'-th trajectory, When the q'th measurement is not associated with the 'tth trajectory, .

[0066] Based on the reassociation results, when a trajectory meets the condition that all targets reassociated with it are associated only with this trajectory, its now_state becomes "RUN," and it is removed from the "relation_other" of other intersecting trajectories. Otherwise, the trajectory's now_state remains "IS." Subsequently, the average of all measurements reassociated with this trajectory is used as its new measurement, and the regular trajectory update steps are performed to achieve the separation of intersecting trajectories.

[0067] Example 1 A two-dimensional MIMO ultra-wideband radar was placed against a wall to verify the effectiveness of the algorithm. The radar transmitted an ultra-wideband signal with a start frequency of 2.5 GHz, a bandwidth of 1 GHz, and a pulse duration of 1 ms to detect multiple moving targets behind the wall. The azimuth array aperture of the two-dimensional MIMO radar was 1.544 m. The wall was a 0.24 m thick air-brick wall. The simulation parameters are shown in Table 2.

[0068] Table 2 Simulation Parameter Table

[0069] The proposed method is compared with existing tracking methods, including the mean shift method based on radar images, the scale-adaptive rotating kernel correlation filter method, and the motion state-based interactive multi-model-Kalman filter tracking method (which uses the JPDA method to achieve data association).

[0070] like Figure 5 As shown, there are two moving targets behind the wall. Target 1 starts at (0, 7) and moves clockwise on a circle with center (0, 5) and radius 2m at a speed of 1.5m / s. Target 2 starts at (4, 5.5) and moves in the negative x-axis direction along the line y=5.5 at a speed of 2m / s. Both targets reach the intersection of their trajectories simultaneously. During the simulation, a random positional jitter of 0.2m was added to simulate the random swaying that occurs during human movement. Simultaneously, to simulate defocusing in the through-wall radar image due to the non-uniformity of the wall, a 90° phase shift was randomly added to all channels of certain frames. A total of 200 frames of radar data were collected in this simulation scenario.

[0071] Tracking results as follows Figure 6 As shown. The tracking results of the comparison method are also shown in... Figure 6As shown in the figure, the proposed algorithm successfully tracked two moving targets. However, while the KCF-based and mean-drift-based methods only have two trajectories, they cannot effectively determine the target's motion state when targets intersect, leading to tracking errors. The proposed method solves this problem by using Bézier curve fitting. The interactive multi-model method tracks multiple discontinuous trajectories. This problem arises because radar images can defocus, leading to false alarms and missed detections. Furthermore, the JPDA data association method can cause target association errors, resulting in trajectory breaks.

[0072] Example 2 A well-designed MIMO radar is placed against a wall to detect moving people behind the wall. The MIMO radar array layout and scene diagram are shown below. Figure 7 As shown, a MIMO with an azimuth array aperture of 0.71m transmits an ultra-wideband signal with a bandwidth of 0.51GHz, penetrating a 24cm thick air brick wall to track moving targets behind the wall.

[0073] Actual test photos and scene settings are as follows Figure 8 As shown in (a) and (b), the two targets intersect during their movement. The tracking results of different methods are as follows... Figure 8 As shown in (c)-(f), the proposed method correctly tracked the trajectories of both targets. However, the KCF and MS methods only tracked the trajectory of target 1 because target 2 only started moving at frame 80. The IMM method is limited by the image defocusing problem, thus producing multiple discontinuous trajectories. The tracking error in the actual test scene 1 is shown below. Figure 8 As shown in (a), in the actual measurement, the tracking error is calculated as the shortest distance from each tracking trajectory point to the moving straight line. It can be seen that the tracking error of the proposed method is consistently low. As shown in Table 2, the tracking accuracy of the proposed method is 90.93%, and the mean square error of tracking is 0.27m, making it the most robust tracking algorithm.

[0074] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. This invention proposes a robust tracking method for multiple moving targets using through-wall radar based on multi-domain data association and multi-state trajectory management, characterized in that... include: Step 1: Trajectory initialization, generating 4 types of trajectory information and 5 states; Step 2: Generate radar imaging results using multi-channel radar echoes; Step 3: Divide the training units and obtain the measurement points using the constant false alarm rate (CFAR) detection method; Step 4: Calculate the multi-domain cost between the existing trajectory and the measurement points, and then establish the minimum cost maximum flow problem to achieve data association by solving this problem; Step 5: Establish a new tracking trajectory using measurements that are not associated with any trajectory, and initialize the trajectory parameters; Step 6: For the trajectory corresponding to the associated measurement, update all trajectory information based on the associated measurement; Step 7: For trajectories not associated with any measurements, update the trajectories based on trajectory information and radar images; Step 8: Generate cross trajectories based on the association results and update the state of the generated cross trajectories; update or separate the cross trajectories based on the measurements; in the cross trajectory separation, first use Bézier curves to fit each cross trajectory, then use the motion state to predict the trajectory, use the prediction results to re-associate with the associated measurements, and finally use the re-association results to separate the cross trajectories.

2. The method as described in claim 1, characterized in that, Step 1, which initializes the tracking trajectory, is performed as follows: The target motion trajectory is assigned four types of information: identification information, tracking quality information, image information, and status information. There are five trajectory states: "initial" indicates the trajectory is in the initialization state; "UB" indicates the trajectory is in the pending start stage, where the trajectory has begun tracking the target but its lifetime is low; "RUN" indicates the trajectory is in normal tracking state, and when the lifetime of the "UB" segment is greater than 0.5, the trajectory state changes to "RUN" (meaning that if the trajectory correctly matches the focused target for 15 consecutive frames (equivalent to 0.14s of human movement for the radar), it is considered a correct moving target trajectory); "IS" indicates the trajectory intersects with other trajectories, and in this case, `relation_other` stores the ID of the intersecting trajectory; "END" indicates the trajectory has stopped.

3. The method as described in claim 1, characterized in that, Step 4, which involves multi-domain data association, is performed as follows: Assume that the total number of existing trajectories in the p-th radar cycle is This includes trajectories in all trajectory states; this module is used to implement... Matching one trajectory with Q measurements; First, calculate the cost between the q-th measurement and the t-th trajectory, denoted as: ； In the formula, , and There are three weighting coefficients, and they satisfy... ; This represents the distance cost between the q-th measurement and the t-th trajectory, used to measure the distance of a human target's movement within a single radar cycle. Its calculation method is as follows: ； In the formula, This represents the tracking result of the t-th trajectory in the (p-1)-th radar cycle; This represents the kernel correlation filter (KCF) response cost between the q-th measurement and the t-th trajectory, used for the target range after extraction and rotation of the q-th measurement position. The correlation between the t-th trajectory and the radar image is calculated as follows: ； In the formula, The last_area represents the t-th trajectory. The filter_para represents the t-th trajectory; Let represent the trajectory lifetime of the t-th trajectory, which represents the confidence level for that trajectory, and is expressed as . ； The data association problem is modeled as a minimum-cost maximum flow problem; trajectory nodes and measurement nodes are added between the source and sink nodes; where trajectory nodes represent existing... A trajectory; measurement nodes represent the existing Q measurements. The capacity and cost between the source point and trajectory nodes are related to the current state of the trajectory. For a trajectory in state 'UB', the capacity of the path between the source point and the trajectory node is 1, and the cost is 1e-6. For a trajectory in state 'RUN', the capacity of the path between the source point and the trajectory node is 1, and the cost is 2e-6. For a trajectory in state 'IS', the capacity of the path between the source point and the trajectory node is the number of trajectories intersecting with this trajectory, and the cost is 2e-6. For a trajectory in state 'END', the capacity of the path between the source point and the trajectory node is 0. The capacity between the source point and the trajectory node represents the maximum number of measurements that each trajectory can be associated with. The cost allocation ensures that measurement values ​​are preferentially allocated to trajectories in motion rather than those tentatively starting. The capacity between trajectory nodes and measurement nodes is 1, indicating that each trajectory and each measurement are matched at most once; the cost is the cost between the corresponding trajectory and measurement. The cost between the measurement node and the sink is 0. The capacity represents the maximum number of trajectories that a target can be associated with, and this value affects how many trajectories will intersect. Considering radar resolution and human body size, this value is usually taken as 3. Data association is achieved by solving the aforementioned minimum-cost maximum-flow problem, ultimately yielding the data association matrix DA. When the q-th measurement is associated with the t-th trajectory... When the q-th measurement is associated with the t-th trajectory, .

4. The method as described in claim 1, characterized in that, Step 5, generating the new trajectory, is performed as follows: When the qth measurement of the p-th frame satisfies the following formula ； Then, a new trajectory needs to be generated starting from the q-th measurement; that is, the q-th measurement. The initial value of the trajectory is ; Extract the target region corresponding to the target. ； Among them, set Let the elliptical range corresponding to the new trajectory in frame p be represented as follows: ； In the formula, , and These represent the major axis, minor axis, and rotation angle of the ellipse corresponding to the measurement position, respectively; given... In this case, the target region is rotated using bilinear interpolation, represented as... ； Construct a correlation filter, where the filter parameters are expressed as follows: ； In the formula, and These represent the forward and inverse Fourier transforms, respectively. This represents a two-dimensional Gaussian function, whose dimensions are... Maintain consistency This indicates the conjugate operation, and σ represents the regularization parameter. Represents an operator that is computed as ； In the formula, β is a constant representing the bandwidth of the Gaussian space; Therefore, the last_area corresponding to this trajectory is obtained. and filter_para In addition, the ID of this trajectory is Both begin and end are p, total is 1, SAN is 1, CLN is 0, and now_state is 'UB'.

5. The method as described in claim 1, characterized in that, Step 6, updating the regular trajectory, is performed as follows: For the t-th trajectory, when its now_state is not 'IS', and the trajectory and its matched q-th measurement satisfy the following formula, ； Therefore, the trajectory is a normal trajectory; For regular trajectories, the results obtained using the KCF method are used as the prediction results; Based on the target position in the previous frame Extract the target region ； Then, the rotated target area is used. To prevent discrepancies in target region size due to varying target locations, the extracted target region needs to be reconstructed. ； In the formula, Indicates the target area in the first frame The dimensionality; using the reconstructed target region, the relevant response is calculated and represented as... ； Maximum value in the relevant response The position represents the estimated target position in frame p. Finally, using this location, the positioning result is obtained. and the target area corresponding to this positioning result. and updated filter parameters Based on the results, the regular trajectory is updated, where the total number of the trajectory is increased by 1, the end is set to p, the SAN number is increased by 1, the CLN is set to 0, and the new last_area is represented as... ； The new filter_para is represented as ； The radar positioning result for this frame is represented as follows: 。 6. The method as described in claim 1, characterized in that, Step 7, updating the unmatched trajectory, is performed as follows: For the t-th trajectory, when its now_state is not 'IS' and it satisfies the following formula: ； Then the trajectory did not match any measurement points; when the trajectory satisfies now_state is 'UB' or CLN is greater than or equal to If this happens, its now_state is changed to 'END'; otherwise, the total number of the trajectory is reduced. The tracking results suggest that The result was discarded because it was invalid; If the trajectory now_state is not 'UB' and CLN is less than 'UB' If the trajectory is not found, KCF is used to continue tracking it, and the KCF method localization result is obtained. and the maximum value of the KCF response Based on the results, unassociated trajectories are updated, where the end of the trajectory is set to p, CLN is incremented by 1, last_area and filter_para remain unchanged, and now_state is updated; the radar tracking result for this frame is represented as follows. 。 7. The method as described in claim 1, characterized in that, Step 8, which involves generating and updating the intersection trajectory, is performed as follows: The now_state of the intersecting trajectory is denoted as 'IS'; One of the conditions for trajectory intersection is that a certain measurement q matches more than one trajectory with now_state 'RUN'; when a certain measurement matches multiple trajectories, the trajectory with now_state 'UB' needs to be changed to 'END'. Assume the number of intersecting trajectories is When a trajectory intersection occurs, the intersection... The current_state of all trajectories is converted to "IS", and the IDs of the remaining trajectories are placed in the relation_other of that trajectory. Following the normal trajectory update steps, the measurement q is considered to be the measurement corresponding to the intersecting trajectory and updated accordingly. A trajectory; For trajectories that have already intersected, the number of measurements that need to be matched to these trajectories needs to be determined. Make a judgment; perform corresponding operations for different measurement quantities; when the measurement quantity When the number of measurements is [number missing], update different measurement trajectories according to the criteria for unmatched measurements; when the number of measurements is [number missing], [number missing] measurements are performed. At that time, the measurement was considered as follows: The measurement corresponding to each trajectory is performed, and the regular trajectory update steps are executed. The difference is that only the trajectory with the longest trajectory life and the largest total number increases the SAN number by 1, while the CLN of the other trajectories increases by 1. Then, if the CLN of a certain trajectory is greater than the gateCLN, the now_state of that trajectory is changed to "END". When the cross trajectory matches the measurement quantity When multiple trajectories need to be separated, the target separation is achieved by utilizing the target's motion state. However, since human targets are mobile targets with variable motion states, they cannot be described by simple motion models, so tracking methods based on motion models are not applicable. In the literature, it has been proven that the acceleration curve of human target motion is smooth under non-external force. Therefore, the target acceleration at the next moment is estimated by using the existing trajectory, and then the target position estimation result is obtained. For a Bézier curve of order Z, it is represented as: ； In the formula, Let z represent the z-th control point. For a Bézier curve of order Z, the number of control points is Z+1. The Bernstein function corresponding to the control point is expressed as: ； For the There are intersecting trajectories, and their total trajectory length is The overall trajectory is represented as Therefore, we obtain ； Solving using the least squares method yields the control points corresponding to the intersection trajectory. Furthermore, the first control point after Bézier curve fitting is obtained. The system generates several intersecting trajectories. Trajectory difference is then used to obtain the velocity and acceleration curves corresponding to each trajectory. Since the acceleration curve of a human target's motion is smooth, the position prediction result of the trajectory is obtained. ; By constructing a new cost function, cross trajectories are achieved. Matched measurements The cost between the q'-th measurement and the t'-th intersection trajectory in a quadratic matching is expressed as: ； Similarly, the minimum-cost maximum flow problem is used to rematch the intersecting trajectories with the measurements; the capacity between the source node and the intersecting trajectory node is 1, and the cost is 0, indicating that each intersecting trajectory requires a rematch process; the capacity between the intersecting trajectory node and the already matched target node is 1, and the cost is... The capacity between the matched target node and the sink node is... The cost is 0; similarly, when the q'th measurement is associated with the t'th trajectory, ; When the q'th measurement is not associated with the 'tth trajectory ; Based on the reassociation results, when a trajectory meets the condition that all targets reassociated with it are associated only with this trajectory, the current_state of this trajectory becomes "RUN" and it is removed from the "relation_other" of other intersecting trajectories. Otherwise, the current_state of this trajectory remains "IS". Subsequently, the average value of all measurements reassociated with this trajectory is used as its new measurement, and the regular trajectory update steps are performed to achieve the separation of intersecting trajectories.

Images