Self-adaptive random finite set tracking method based on low-cost spaceborne distributed radar

Abstract

The invention discloses a self-adaptive random finite set tracking method based on a low-cost spaceborne distributed radar, and mainly solves the problems of high cost and poor self-adaptive degree in the prior art. According to the implementation scheme, an equivalent one-transmitting-multiple-receiving model is constructed in a cosatellite forwarding mode, and echo data are recorded at low cost; forming two equivalent-length synthetic apertures according to the echo data; measuring points are obtained through synthetic aperture imaging, clutter suppression and detection; self-adaptively generating new target strength by using the multi-step measuring points; the clutter rate under each step is adaptively estimated by using the multi-step measuring points; and performing random finite set target tracking by using the adaptively generated intensity and clutter rate. The method is low in cost and short in coherent accumulation time, performance reduction caused by parameter-free prior can be avoided on the premise of keeping real-time tracking processing, engineering can be further achieved, and the method can be used for real-time detection tracking of low-cost distributed radar ground moving targets or low-altitude moving targets.

Description

Technical Field

[0001] The invention belongs to the field of radar technology, and in particular relates to an adaptive random finite set tracking method, which can be used for real-time detection and tracking of ground moving targets or low-altitude moving targets of low-cost distributed radars. Background Art

[0002] Spaceborne distributed radars conduct collaborative monitoring of observation scenes through space-based networking. Through the spatial deployment of each satellite and the design of signal transmission and reception modes in combination with specific missions, imaging, detection and other needs can be achieved. It has many advantages such as wide detection range, strong survivability, and no national or regional restrictions. However, spaceborne distributed radars are generally expensive and not conducive to large-scale deployment.

[0003] In order to alleviate the high cost of satellite-borne distributed radar, the equivalent SIMO system can be adopted. The main satellite transmits RF signals to illuminate the detection area, and multiple auxiliary satellites receive echo signals at the same time. In this system, the receiving auxiliary satellite only serves as a relay route. It does not sample or store the received echo signals, but modulates and forwards them to the main satellite, which is sampled after orthogonal reception by the main satellite. Both transmission and reception are completed on the main satellite. Compared with the traditional one-transmit-multiple-receive system, its cost is reduced and the structure is simplified, which is convenient for large-scale configuration of auxiliary satellites. After preprocessing on the main satellite, the echo signal is detected by joint ground targets. After target detection processing, the measurement data of each step is provided for subsequent target tracking processing.

[0004] Existing target tracking methods can be divided into data association processing methods and random finite set processing methods. Among them, the random finite set RFS target tracking method can achieve efficient target tracking in complex environments where the number of targets is time-varying and unknown. However, the traditional RFS target tracking method is implemented under the assumption that the probability density of the new target is known a priori and the clutter is given a priori, which greatly limits the further application and engineering implementation of the algorithm.

[0005] At present, there are mainly the following methods for tracking random finite sets of probability hypothesis density of adaptive new targets:

[0006] B. Ristic et al. proposed a measurement-driven adaptive random finite set tracking method for newborn targets in the TAES paper "Adaptive Target Birth Intensity for PHD and CPHD Filters". The core idea is to introduce labels to distinguish newborn targets from surviving targets, construct intensity prediction update formulas for newborn targets and surviving targets respectively, and model the PHD intensity of newborn targets in a measurement-driven way, and then use the first-order Bayesian moment approximation formula for prediction update.

[0007] In the TII paper “GM-PHD-based multi-target visual tracking using entropy distribution and game theory”, Xiaolong Zhou et al. proposed a nascent intensity estimation algorithm based on entropy distribution and coverage to automatically and accurately track nascent targets in noisy videos.

[0008] The above method only considers the problem of adaptive new target intensity modeling under single-step observation, and does not consider the comprehensive use of multi-step observation information, so it does not make full use of time dimension information. In addition, the above method only considers the problem of adaptive new generation, and does not consider the problem of adaptive estimation of clutter rate at the same time, so the degree of adaptability is not strong enough. Summary of the invention

[0009] The purpose of the present invention is to address the deficiencies of the above-mentioned prior art and propose a low-cost spaceborne distributed radar adaptive random finite set tracking method, so as to reduce the performance degradation caused by the absence of parameter priors and improve the degree of tracking adaptability by iteratively performing adaptive new intensity modeling and parameter estimation using multi-step measurements.

[0010] The technical solution for achieving the purpose of the present invention comprises the following steps:

[0011] (1) A satellite-borne distributed radar cooperative detection model is constructed, which includes a transmitting primary satellite and N receiving auxiliary satellites, where each auxiliary satellite has two x The receiving channel of the satellite uses the sliding beam mode to continuously observe the target, and each satellite receives the scattered echo signal.

[0012] (2) Design the method of relaying by the auxiliary satellite and orthogonal reception by the primary satellite of the satellite-borne distributed radar system, perform coherent mixing and AD sampling on the echo signal, and obtain the dual-path echo signal Ω of N auxiliary satellites on the primary satellite;

[0013] (3) In the kth observation, two equivalent long synthetic apertures that meet the DPCA conditions of the offset phase center antenna are constructed by sequentially connecting the dual-path echo signals of N auxiliary satellites. and The symbol || indicates serial splicing and the aperture spatial positions are overlapped;

[0014] (4) Based on two long synthetic apertures and After SAR imaging, DPCA clutter suppression and constant false alarm detection, the moving target measurement position set Z of the kth observation is obtained. k ;

[0015] (5) Extracting the target state through sliding window correlation of multi-step measurements and the target state covariance matrix And design a weight mapping function to adaptively generate the new target intensity v of the kth observation nb,k (x);

[0016] (6) Record the total number of measurements in the multi-step measurement sliding window and adaptively estimate the clutter rate κ of the kth observation k ;

[0017] (7) Using Γ obtained in step (5) k and κ obtained in step (6) k , for the moving target measured position set Z of the kth observation k Perform adaptive random finite set filtering to obtain the posterior strength v of the kth observation k|k (x), extract and output the multi-target tracking result of the kth observation in real time from the intensity;

[0018] (8) Repeat steps (2) to (7) to achieve real-time continuous detection and tracking of multiple targets in the observation scene.

[0019] Compared with the prior art, the present invention has the following advantages:

[0020] 1) The present invention designs a new system of satellite-borne distributed radar signal transceiver system, adopts the equivalent one-transmit-multiple-receive mode of primary satellite transmission-secondary satellite forwarding-primary satellite reception, and can achieve low-cost data acquisition; by designing a satellite formation, two equivalent synthetic apertures are formed, which can realize ground moving target indication, shorten the coherent accumulation time, and realize rapid moving target detection;

[0021] 2) Since the present invention adopts adaptive random finite set tracking in multi-target tracking processing, it can jointly utilize multi-step measurement information to iteratively estimate the intensity and clutter rate of new targets in real time, reduce the performance degradation caused by parameter-free priors, and improve adaptive tracking performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is the overall flow chart of the implementation of the present invention;

[0023] Figure 2 It is a sub-flow chart of iterative multi-step tracking and processing of measurement data in the present invention;

[0024] Figure 3 It is a comparison chart of OSPA indicators of tracking simulation experiments performed using the present invention and the existing method;

[0025] Figure 4 It is a comparison chart of potential estimation indexes of tracking simulation experiments performed using the present invention and the existing method. DETAILED DESCRIPTION

[0026] In order to enable those skilled in the art to better understand the scheme of the present invention, the technical scheme in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, other embodiments obtained by ordinary technicians in this field without creative work should all fall within the scope of protection of the present invention.

[0027] It should be noted that the step numbers in the specification and claims of the present invention are only for a clear description of the implementation scheme of the present invention to facilitate understanding, and the order of the step numbers is not limited.

[0028] Reference Figure 1 and Figure 2 , the implementation steps of the present invention are as follows:

[0029] Step 1: Build a satellite-borne distributed radar collaborative detection model.

[0030] The satellite-borne distributed radar cooperative detection model constructed in this step includes one transmitting main satellite and N receiving auxiliary satellites;

[0031] The primary satellite continuously observes the target in a sliding beam mode, with the starting position at (x, y, z) and the ending position at (x, y + L, z), where L is the short aperture length;

[0032] The N receiving auxiliary satellites all have two receiving channels with a spacing of half a wavelength for receiving scattered echo signals. These auxiliary satellites are arranged at equal intervals. The starting position of the i-th auxiliary satellite is (x, y+(2i-1)L, z), and the ending position is expressed as (x, y+2iL, z), i=1, 2,…, N.

[0033] Step 2: Design a satellite-borne distributed radar system solution with auxiliary satellite relay forwarding and orthogonal reception by the primary satellite, and obtain the dual-path echo signal Ω of N auxiliary satellites on the primary satellite.

[0034] 2.1) The primary satellite sends a direct reference signal f to the nth secondary satellite rn ;

[0035] 2.2) The echo received by the ith receiving channel of the nth auxiliary satellite With the direct reference signal f rn Superposition is performed to obtain a mixed signal

[0036] 2.3) Use a coherent mixer to mix the signals The frequency shift is The up-mixed signal is obtained by frequency shift and forwarded to the primary satellite, i = 1, 2;

[0037] 2.4) After receiving the forwarded signal, the primary satellite uses a coherent mixer to perform a frequency shift of The frequency shift is to remove the orthogonal frequency modulation of the nth satellite, and AD sampling is performed using the coherent I / Q demodulator driven by the local oscillator, i = 1, 2;

[0038] 2.5) Perform range compression and azimuth Fourier transform on the digital signal after AD sampling, and separate the zero-frequency direct reference signal in the range Doppler domain. and echo signal

[0039] 2.6) Estimating the two reference signals f in the primary satellite rn and The amplitude and phase errors of the echo signal Compensation is performed to obtain the echo signal after correcting the up-down mixing error.

[0040] 2.7) Repeat steps 2.1) to 2.6) to obtain the dual-path echo signal Ω of N auxiliary satellites on the primary satellite:

[0041]

[0042] Where n=1,...,N.

[0043] Step 3: Obtain two equivalent long apertures based on the echo signals of N auxiliary satellites and

[0044] In the kth observation, two equivalent long synthetic apertures satisfying the DPCA conditions are constructed by sequentially connecting the echo signals of N auxiliary satellites. and

[0045]

[0046] The symbol || indicates serial splicing and overlapping positions.

[0047] Step 4: Based on two long synthetic apertures and Get the moving target measurement position set Z of the kth observation k .

[0048] 4.1) Use polar coordinate formatting algorithm to and The echo data in the image is focused to obtain the corresponding SAR image and

[0049] 4.2) Obtaining a difference image based on complex image cancellation Among them, i x ,i y is the two-dimensional pixel index of the difference image, i x =1,...,N x ,i y =1,...,N y , N x and N y Respectively represent the number of resolution units in the azimuth and range directions of the difference SAR image;

[0050] 4.3) Yes and Perform interference processing to obtain interference images in(·) * It means taking the conjugate of the image;

[0051] 4.4) Let the moving target measurement position set be Z k , set the threshold V T , using constant false alarm detection to make a decision:

[0052] If ΔI(i x ,i y )≥V T , the target exists and the measurement point is recorded Execute step 4.5);

[0053] If ΔI(i x ,i y )<V T , then the target does not exist;

[0054] in is the difference image ΔI(i x ,i y ) in i x ,i y The corresponding two-dimensional position under the index, the superscript q is the measurement point index;

[0055] 4.5) From the measurement point Extract the interference phase at the corresponding interference image position

[0056] 4.6) Calculate the measurement points The azimuth offset Where y0 is the slant distance from the center of the scene;

[0057] 4.7) Compensate the offset of the measuring point to obtain the relocated measuring point

[0058]

[0059] in To relocate the azimuth position of the moving target, To determine the distance to the moving target after repositioning;

[0060] 4.8) Repeat steps 4.4) to 4.7) until all the two-dimensional pixel indexes of the difference image are traversed to obtain the moving target measurement position set m k is the total number of measurements for the kth observation.

[0061] Step 5: Extract the target state through sliding window correlation of multi-step measurements and the target state covariance matrix And design a weight mapping function to adaptively generate the new target intensity v of the kth observation nb,k (x).

[0062] 5.1) Let L ≥ 3 be the sliding window length and compare it with the number of observation steps k:

[0063] If the number of observation steps k<L, then execute step 5.2),

[0064] If the number of observation steps k ≥ L, then execute step 5.3),

[0065] 5.2) Store the measured position set Z of the current observation step k , execute step eight to process the next observation step;

[0066] 5.3) Obtain the index set of the association pairs from the k-L+1th observation to the kth observation through the nearest neighbor measurement association Where M is the index set Π (k-L+1):k The number of associated pairs in , represents the measurement index in the measurement location set with observation step j for association pair i, Represents Π (k-L+1):k The index of the i-th associated pair in , j = k-L+1,…,k-1,k;

[0067] 5.4) According to the index set Π of the association pairs from the k-L+1th observation to the kth observation (k-L+1):k The associated pair index p i , extract the corresponding measurements in each observation step and form the correlation pair measurement matrix and the corresponding set of measurement error covariance matrices for the associated pairs in is the measurement corresponding to pair i in the measurement position set under the j-th observation step, is the measurement error covariance matrix corresponding to the measurement position set association pair i under the j-th observation step;

[0068] 5.5) Estimate the position, velocity, and acceleration of the target at the kth observation based on the measurement matrix. Let T be the time interval between two adjacent observation steps, and let Calculate the estimated target velocity for the kth observation Estimated acceleration and the target estimated speed at the k-1th observation step

[0069]

[0070] in are the vectors corresponding to the kth, k-1th and k-2th columns in the correlation pair measurement matrix, are the components of the target estimated velocity in the x and y directions, are the components of the target estimated acceleration in the x and y directions, q=1,2,…,J nb,k , J nb,k is the number of new target Gaussian components in the kth observation step;

[0071] 5.6) The measured value in step 5.4) and the target estimated speed in step 5.5) Estimated acceleration Write the estimated new target state in the form of the existing CA model state vector

[0072]

[0073] in[·] T Indicates transposition processing;

[0074] 5.7) Let H k is the observation matrix of the kth observation, according to the qth measurement error covariance matrix of the kth observation Get the new target state covariance matrix Among them, β1 and β2 are constant factors;

[0075] 5.8) The prediction strength of the k-th observation is expressed in Gaussian mixture form as Calculates Gaussian distribution probability density function

[0076]

[0077] Among them, J k|k-1 is the number of Gaussian components for the predicted intensity, are the weight, mean and covariance matrix of the jth Gaussian component of the predicted intensity, j = 1, 2, ..., J k|k-1 , is the new target state obtained in step 5.6), n z is the measurement dimension;

[0078] 5.9) Calculation The accumulated value of q :

[0079]

[0080] 5.10) Define the mapping function as f(x) = e -βx , where β is a scaling factor that weights the Gaussian mixture component of the new density Set to Λ q The mapping function is

[0081] 5.11) Take the maximum value of the newborn density weight

[0082] 5.12) Set the upper limit of the new student target weight

[0083] 5.13) Define limiting factors

[0084] 5.14) Scale the Gaussian mixture component weights of the new density to obtain the scaled Gaussian mixture component weights of the new density:

[0085] 5.15) According to the new target state obtained in step 5.6) Step 5.7) The new target state covariance matrix obtained The Gaussian mixture component weights of the new density obtained in step 5.14) The new target intensity v observed at the kth time is converted into nb,k (x) is represented as follows:

[0086]

[0087] in, is the Gaussian distribution symbol.

[0088] Step 6: Record the total number of measurements in the multi-step measurement sliding window and adaptively estimate the clutter rate κ of the kth observation k .

[0089] 6.1) Record the number of measurements m in each step within the sliding window step length L k′ , and measure the number m of each step k′ Sum and accumulate to get the total number of measurements m within the sliding window step sum :

[0090]

[0091] Where k′=k-L+1,…,k-1,k;

[0092] 6.2) Record the potential estimate number n of the first L-1 steps within the sliding window step size L k″ , get the total number of potential estimates within the sliding window step n sum :

[0093]

[0094] Where, k″=k-L+1,…,k-1;

[0095] 6.3) Let p c is the probability density function of the clutter, according to the total amount m sum Total number of potential estimates n sum , get the estimated clutter rate κ of the kth observation k :

[0096]

[0097] Step 7: Use the Γ obtained in step 5 k and κ obtained in step 6 k , for the k-th observation Z k Perform adaptive random finite set filtering to obtain the posterior strength v of the kth observation k|k (x), the multi-target tracking result of the k-th observation is extracted and output in real time from the intensity.

[0098] 7.1) Obtain the posterior strength v expressed in the form of a Gaussian mixture under the k-1th observation k-1|k-1 (x):

[0099]

[0100] Among them, J k-1|k-1 is the posterior strength v of the k-1th observation k-1|k-1 The number of Gaussian components of (x), are the posterior strength v of the k-1th observation respectively. k-1|k-1 The weight, mean and covariance matrix of the i-th Gaussian component of (x), i=1,2,…,J k-1|k-1 ;

[0101] 7.2) The posterior intensity v expressed in the form of Gaussian mixture for the k-1th observation k-1|k-1 The mean of each Gaussian component of (x) Perform Kalman filter prediction to obtain the predicted mean of each Gaussian component

[0102]

[0103] Among them, F k,k-1 is the state transition matrix from the k-1th observation to the kth observation;

[0104] 7.3) The posterior intensity v expressed in the form of Gaussian mixture for the k-1th observation k-1|k-1 The covariance matrix of each Gaussian component of (x) Perform Kalman filter prediction to obtain the prediction covariance matrix of each Gaussian component

[0105]

[0106] Among them, Q k,k-1 is the error covariance matrix of the state transfer from the k-1th observation to the kth observation;

[0107] 7.4) Using the predicted mean and the prediction covariance matrix and the new target strength v obtained in step 5 nb,k (x), obtain the predicted strength v expressed in the form of a Gaussian mixture from the k-1th observation to the kth observation k|k-1 (x):

[0108]

[0109] Among them, p S,k is the survival probability of the target under the kth observation, is the Gaussian distribution symbol, J k|k-1 is the predicted strength v of the kth observation k|k-1 The number of Gaussian components of (x), are the predicted strength v of the kth observation respectively. k|k-1 The weight, mean and covariance matrix of the i-th Gaussian component of (x), i=1,2,…,J k|k-1 ;

[0110] 7.5) Using the predicted strength v expressed in the form of a Gaussian mixture from the k-1th observation to the kth observation k|k-1 The covariance matrix of each Gaussian component of (x) Get the innovation covariance matrix of each Gaussian component and Kalman gain

[0111]

[0112] Among them, H k is the observation matrix under the kth observation, R k is the measurement error covariance matrix;

[0113] 7.6) The predicted strength v from the k-1th observation to the kth observation is expressed in the form of a Gaussian mixture k|k-1 The mean of each Gaussian component of (x) Perform Kalman filter update to obtain the updated mean of each Gaussian component

[0114]

[0115] Where z is the measurement position set Z k The vector elements in ;

[0116] 7.7) The predicted strength v from the k-1th observation to the kth observation is expressed in the form of a Gaussian mixture k|k-1 The covariance matrix of each Gaussian component of (x) Perform Kalman filter update to obtain the updated covariance matrix of each Gaussian component

[0117]

[0118] 7.8) Using the predicted strength v expressed in the form of a Gaussian mixture from the k-1th observation to the kth observation k|k-1 The covariance matrix of each Gaussian component of (x) And the clutter rate κ estimated in step 6 k , get the updated weights for each Gaussian component

[0119]

[0120] Among them, p D,k is the detection probability of the target under the k-th observation step;

[0121] 7.9) Using the predicted strength v given in step 7.4) k|k-1 (x) and steps 7.5) to 7.8) for the measurement position set Z k Updated intensity of element z Get the posterior strength v of the kth observation k|k (x):

[0122]

[0123] Among them, J k|k is the posterior strength v of the kth observation k|k (x) the number of Gaussian components;

[0124] 7.10) From v k|k (x) extracts and outputs the multi-target tracking result of the kth observation.

[0125] Step 8: Repeat steps 2 to 7 to achieve real-time continuous detection and tracking of multiple targets.

[0126] The effects of the present invention are further described below in conjunction with simulation.

[0127] 1. Simulation parameters:

[0128] The number of Monte Carlo experiments is set to 1000, the total observation time is 50s, the time interval between two adjacent observation steps is 1s, the target detection probability is 0.99, the target survival probability is 0.98, the maximum number of Gaussian constraints is 100, and the pruning threshold is 10. -5 , the merging threshold is 4.

[0129] 2. Simulation content:

[0130] Simulation 1: Under the above simulation parameters, the adaptive random finite set tracking is performed using the present invention and the existing 4×4 grid adaptive newborn processing algorithm, 8×8 grid adaptive newborn processing algorithm, and Ristic adaptive newborn processing algorithm to obtain the optimal sub-mode allocation OSPA distance curves for each tracking, which are used to evaluate the accuracy of target tracking state estimation. The results are shown in FIG. Figure 3 shown.

[0131] from Figure 3 It can be seen that the distance error curve of the present invention is the lowest, which is better than the three existing methods, indicating that the target tracking state estimation of the present invention is more accurate.

[0132] Simulation 2: Under the above simulation parameters, the adaptive random finite set tracking is performed using the present invention and the existing 4×4 grid adaptive newborn processing algorithm, 8×8 grid adaptive newborn processing algorithm, and Ristic adaptive newborn processing algorithm. The target number estimation performance of each algorithm is evaluated, and the potential estimation curves of each comparative analysis tracking are obtained. The results are shown in FIG. Figure 4 shown.

[0133] from Figure 4 It can be seen that the potential estimation curve of the present invention is closest to the true value curve, which is better than the above three methods, and the target number estimation is more accurate.

Claims

1. An adaptive random finite set tracking method based on low-cost spaceborne distributed radar, characterized in that: The steps include: (1) A satellite-borne distributed radar cooperative detection model is constructed, which includes a transmitting primary satellite and N receiving auxiliary satellites, where each auxiliary satellite has two x The receiving channel of the satellite uses the sliding beam mode to continuously observe the target, and each satellite receives the scattered echo signal. (2) Design the method of relaying by the auxiliary satellite and orthogonal reception by the primary satellite of the satellite-borne distributed radar system, perform coherent mixing and AD sampling on the echo signal, and obtain the dual-path echo signal Ω of N auxiliary satellites on the primary satellite; (3) In the kth observation, two equivalent long synthetic apertures that meet the DPCA conditions of the offset phase center antenna are constructed by sequentially connecting the dual-path echo signals of N auxiliary satellites. and The symbol || indicates serial splicing and the aperture spatial positions are overlapped; (4) Based on two long synthetic apertures and After SAR imaging, DPCA clutter suppression and constant false alarm detection, the moving target measurement position set Z of the kth observation is obtained. k ; (5) Extracting the target state through sliding window correlation of multi-step measurements and the target state covariance matrix And design a weight mapping function to adaptively generate the new target intensity v of the kth observation nb,k (x); (6) Record the total number of measurements in the multi-step measurement sliding window and adaptively estimate the clutter rate κ of the kth observation k ; (7) Using Γ obtained in step (5) k and κ obtained in step (6) k , for the moving target measured position set Z of the kth observation k Perform adaptive random finite set filtering to obtain the posterior strength v of the kth observation k|k (x), extract and output the multi-target tracking result of the kth observation in real time from the intensity; (8) Repeat steps (2) to (7) to achieve real-time continuous detection and tracking of multiple targets in the observation scene.

2. The method according to claim 1, characterized in that Step (1) constitutes the transmitting main satellite and N receiving auxiliary satellites in the satellite-borne distributed radar cooperative detection model, and their distribution positions are as follows: The starting position of the primary star is (x, y, z) and the ending position is (x, y + L, z), where L is the short aperture length; N auxiliary stars are arranged at equal intervals. The starting position of the i-th auxiliary star is (x, y+(2i-1)L, z), and the ending position is expressed as (x, y+2iL, z), i=1,2,…,N.

3. The method according to claim 1, characterized in that In step (2), the dual-path echo signals Ω of the N auxiliary satellites are obtained from the primary satellite, and the implementation steps include the following: 2a) The primary satellite sends a direct reference signal f to the nth secondary satellite rn ; 2b) The echo received by the ith receiving channel of the nth auxiliary satellite With the direct reference signal f rn Superposition is performed to obtain a mixed signal 2c) Use a coherent mixer to mix the signals The frequency shift is The up-mixed signal is obtained by frequency shift and forwarded to the primary satellite, i = 1, 2; 2d) After receiving the forwarded signal, the primary satellite uses a coherent mixer to perform a frequency shift of The frequency shift is used to remove the orthogonal frequency modulation of the nth auxiliary satellite, and AD sampling is performed using a coherent I / Q demodulator driven by a local oscillator, i = 1, 2; 2g) Perform range compression and azimuth Fourier transform on the digital signal after AD sampling, and separate the zero-frequency direct reference signal in the range Doppler domain and echo signal 2f) Estimating the two reference signals f in the primary satellite rn and The amplitude and phase errors of the echo signal Compensation is performed to obtain the echo signal after correcting the up-down mixing error. 2h) Repeat steps 2a) to 2f) to obtain dual-path echo signals of N auxiliary satellites on the primary satellite.

4. The method according to claim 1, characterized in that In step (4), the moving target measurement position set Z of the kth observation is obtained k , the implementation steps include the following: 4a) Use polar coordinate formatting algorithm to and The echo data in the image is focused to obtain the corresponding SAR image and 4b) Obtaining a difference image based on complex image cancellation Among them, i x ,i y is the two-dimensional pixel index of the difference image, i x =1,...,N x ,i y =1,...,N y , N x and N y Respectively represent the number of resolution units in the azimuth and range directions of the difference SAR image; 4c) Yes and Perform interference processing to obtain interference images in(·) * It means taking the conjugate of the image; 4d) Let the moving target measurement position set be Z k , set the threshold V T , using constant false alarm detection to make a decision: If ΔI(i x ,i y )≥V T , the target exists and the measurement point is recorded Execute step 4e); If ΔI(i x ,i y )<V T , then the target does not exist; in is the difference image ΔI(i x ,i y ) in i x ,i y The corresponding two-dimensional position under the index, the superscript q is the measurement point index; 4e) From the measurement point Extract the interference phase at the corresponding interference image position 4f) Calculate the measurement points The azimuth offset Where y0 is the slant distance from the center of the scene; 4g) Compensate the offset of the measuring point to obtain the relocated measuring point in To relocate the azimuth position of the moving target, To determine the distance to the moving target after repositioning; 4h) Repeat steps 4d) to 4g) until the two-dimensional pixel index of the difference image is traversed to obtain the moving target measurement position set m k is the total number of measurements for the kth observation.

5. The method according to claim 1, characterized in that In step (5), the target state is extracted by sliding window correlation of multi-step measurements and the target state covariance matrix The implementation steps include the following: 5a) Let L ≥ 3 be the sliding window length and compare it with the number of observation steps k: If the number of observation steps k<L, then execute step 5b), If the number of observation steps k ≥ L, execute step 5c); 5b) Store the measured position set Z of the current observation step k , implement step (8) in claim 1, and proceed to the next observation step; 5c) Obtain the index set of the measurement association pairs from the k-L+1th observation to the kth observation through the nearest neighbor measurement association: Where M is the index set Π k-L+1:k The number of associated pairs in , represents the index of the i-th associated pair, represents the measurement index in the measurement position set with observation step j for association pair i, j=k-L+1,…,k-1,k; 5d) Measure the index set Π of the association pairs from the k-L+1th observation to the kth observation k-L+1:k The associated pair index p i , extract the corresponding measurements in each observation step and form the correlation pair measurement matrix and the corresponding set of measurement error covariance matrices for the associated pairs in is the measurement corresponding to pair i in the measurement position set under the j-th observation step, is the measurement error covariance matrix corresponding to the measurement position set association pair i under the j-th observation step; 5e) Estimate the position, velocity, and acceleration of the target at the kth observation based on the measurement matrix. Let T be the time interval between two adjacent observation steps, and let Calculate the target estimated velocity for the kth observation Estimated acceleration and the target estimated speed at the k-1th observation step in are the vectors corresponding to the kth, k-1th and k-2th columns in the correlation pair measurement matrix, are the components of the target estimated velocity in the x and y directions, Estimate the components of the target acceleration in the x and y directions respectively; 5f) The measured value in step 5d) and the target estimated speed in step 5e) Estimated acceleration Write the estimated new target state in the form of the existing CA model state vector 5g) Set H k is the observation matrix of the kth observation, according to the qth measurement error covariance matrix of the kth observation Get the new target state covariance matrix Among them, β1 and β2 are constant factors.

6. The method according to claim 1, characterized in that Step (5) Design a weight mapping function to adaptively generate the new target intensity v of the kth observation nb,k (x), the implementation steps include the following: 5h) The prediction strength of the k-th observation is expressed in Gaussian mixture form as Calculates Gaussian distribution probability density function Among them, J k|k-1 is the number of Gaussian components of the predicted intensity, are the weight, mean and covariance matrix of the jth Gaussian component of the predicted intensity, j = 1, 2, ..., J k|k-1 , is the new target state obtained in step 5f), n z is the measurement dimension; 5i) Calculation The accumulated value of q : Among them J nb,k is the number of new target Gaussian components in the kth observation step; 5j) Define the mapping function as f(x) = e -βx , where β is a scaling factor that weights the Gaussian mixture component of the new density Set to Λ q The mapping function is 5k) Take the maximum value of the newborn density weight 5l) Set the upper limit of the weight of the new student target 5m) Define limiting factors 5n) Scale the Gaussian mixture component weights of the new density to obtain the scaled Gaussian mixture component weights of the new density: 5o) The new target state obtained according to step 5f) Step 5g) The new target state covariance matrix obtained Step 5n) Gaussian mixture component weights of the new density obtained The new target intensity v observed at the kth time is converted into nb,k (x) is represented as follows: in, is the Gaussian distribution symbol.

7. The method according to claim 1, characterized in that In step (6), the clutter rate κ of the kth observation is adaptively estimated k , the implementation steps include the following: 6a) Record the number of measurements m in each step within the sliding window step length L k′ , and sum and accumulate to get the total number of measurements m within the sliding window step sum : Where k′=k-L+1,…,k-1,k; 6b) Record the number of potential estimates n in the previous L-1 steps within the sliding window step size L k″ , get the total number of potential estimates within the sliding window step n sum : Where, k″=k-L+1,…,k-1; 6c) Let p c is the probability density function of the clutter, according to the total amount m sum Total number of potential estimates n sum , get the estimated clutter rate κ of the kth observation k :

8. The method according to claim 1, characterized in that In step (7), the moving target position set Z of the kth observation is k Perform adaptive random finite set filtering to obtain the posterior strength v of the kth observation k|k (x), the implementation steps include the following: 7a) Obtain the posterior intensity v expressed in the form of a Gaussian mixture under the k-1th observation k-1|k-1 (x): Among them, J k-1|k-1 is the posterior strength v of the k-1th observation k-1|k-1 The number of Gaussian components of (x), are the posterior strength v of the k-1th observation respectively. k-1|k-1 The weight, mean and covariance matrix of the i-th Gaussian component of (x), i=1,2,…,J k-1|k-1 ; 7b) The posterior intensity v expressed in the form of Gaussian mixture for the k-1th observation k-1|k-1 The mean of each Gaussian component of (x) Perform Kalman filter prediction to obtain the predicted mean of each Gaussian component Among them, F k,k-1 is the state transition matrix from the k-1th observation to the kth observation; 7c) The posterior intensity v expressed in the form of Gaussian mixture for the k-1th observation k-1|k-1 The covariance matrix of each Gaussian component of (x) Perform Kalman filter prediction to obtain the prediction covariance matrix of each Gaussian component Among them, Q k,k-1 is the error covariance matrix of the state transfer from the k-1th observation to the kth observation; 7d) Using the predicted mean and the prediction covariance matrix and the new target strength v obtained in step (5) nb,k (x), obtain the predicted strength v expressed in the form of a Gaussian mixture from the k-1th observation to the kth observation k|k-1 (x): Among them, p S,k is the survival probability of the target under the kth observation, is the Gaussian distribution symbol, J k|k-1 is the predicted strength v of the kth observation k|k-1 The number of Gaussian components of (x), are the predicted strength v of the kth observation respectively. k|k-1 The weight, mean and covariance matrix of the i-th Gaussian component of (x), i=1,2,…,J k|k-1 ; 7e) Using the clutter rate κ estimated in step (6) k , the predicted strength v expressed in the form of a Gaussian mixture from the k-1th observation to the kth observation k|k-1 The mean of each Gaussian component of (x) and the covariance matrix Perform Kalman filter update and weight calculation to obtain the new information covariance matrix of each Gaussian component Kalman Gain Update mean Update the covariance matrix And update the weights 7f) Using the predicted strength v given in step 7d) k|k-1 (x) and step 7e) for the measurement position set Z k Updated intensity of element z Get the posterior strength v of the kth observation k|k (x): Among them, J k|k is the posterior strength v of the kth observation k|k The number of Gaussian components of (x).

9. The method according to claim 1, characterized in that: Step 7e) Obtain the innovation covariance matrix of each Gaussian component Kalman Gain Update mean Update the covariance matrix and update weights The formula is as follows: Where i = 1, 2, ..., J k|k-1 , H k is the observation matrix under the kth observation, R k is the measurement error covariance matrix, z is the measurement position set Z k The vector elements in p D,k is the detection probability of the target under the k-th observation step.

Images