A high-precision imaging method of distributed synthetic aperture radar based on trajectory reconstruction

Abstract

The application discloses a kind of distributed synthetic aperture radar high-precision imaging methods based on trajectory reconstruction, first solve each platform high-order error in distributed system, ensure the focusing degree of each platform imaging result, then solve the linear error of each platform, utilize the quantitative relationship between each platform linear trajectory error parameter and image shift, linear error parameter is inverted, finally reconstruct distributed platform high-precision motion trajectory, compensate motion error, realize distributed SAR focusing imaging.The method of the application solves the motion error compensation problem in distributed synthetic aperture radar imaging, can high-precisionly realize distributed SAR multi-platform motion error reconstruction, provides premise for multi-platform cooperative focusing imaging, has the advantages of high imaging precision, wide scene applicability.

Description

Technical Field

[0001] This invention belongs to the field of radar imaging technology, specifically relating to a high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction. Background Technology

[0002] Synthetic Aperture Radar (SAR) provides all-weather, all-time remote sensing capabilities and has wide applications in the civilian sector. Resolution and frame rate are two key parameters of SAR. Higher resolution requires a larger azimuth synthetic aperture, which in turn requires a longer synthetic aperture time, leading to a decrease in the imaging frame rate. Therefore, improving both parameters simultaneously for existing monostatic SAR is a challenge. Distributed SAR can overcome the limitations of existing monostatic SAR systems. In distributed SAR imaging, multiple radars work together to quickly synthesize an equivalent long aperture. Specifically, distributed SAR spatially decomposes the aperture synthesis process of existing SAR, dividing the entire synthetic aperture into multiple short aperture segments, each synthesized by a corresponding radar. Therefore, distributed SAR can reduce aperture synthesis time while maintaining aperture length, overcoming the contradiction between high-resolution imaging and high frame rate imaging in microwave band SAR, and has significant research value.

[0003] Currently, the trajectory accuracy acquired by navigation equipment is limited and does not meet the requirements of distributed SAR imaging, introducing range migration error and phase error into the echo. Unlike existing monostatic SAR, distributed SAR has a special geometric configuration, and its range migration error and phase error are discontinuous and segmented. This greatly deteriorates the imaging results and increases the difficulty of error inversion and compensation.

[0004] To achieve motion error inversion and compensation in synthetic aperture radar (SAR), the paper "An Estimation and Compensation Method for Motion Trajectory Error in Bistatic SAR. RemoteSensing. 2022, 14, 5522" proposes a method for estimating and compensating motion trajectory errors in bistatic SAR. This method can extract Doppler parameters of strong scatterers and invert bistatic SAR trajectories, but it cannot reconstruct linear trajectory errors. However, the linear trajectory errors of each platform in distributed SAR lead to higher-order trajectory errors in the overall system, which is significant. Furthermore, this method requires solving for trajectory errors on a per-site basis, which is very time-consuming. The paper "Algorithm on the Estimation of Residual Motion Errors in Airborne SAR Images, in IEEE Transactions on Geoscience and Remote Sensing, vol.52, no.2, pp.1311-1323, Feb.2014" proposes an algorithm for estimating residual motion errors in airborne synthetic aperture radar images. This algorithm solves for the actual trajectory based on the relationship between the phase difference of sub-aperture images and trajectory error parameters. However, it requires a series of isolated strong points distributed along the azimuth, which places relatively strict requirements on the imaging scene. The paper "Microwave Photonic SAR High-Precision Imaging Based on Optimal Subaperture Division, in IEEE Transactions on Geoscience and Remote Sensing, vol.60, pp.1-17, 2022" proposes a high-precision imaging algorithm for microwave photonic SAR based on optimal subaperture division. This algorithm calculates the sub-aperture image offset using image registration and reconstructs the trajectory based on the relationship between the image offset and trajectory error parameters, effectively ensuring imaging accuracy.

[0005] The methods described above can estimate the linear error of the trajectory, including the intercept and slope, but the error model is insufficient to describe the trajectory error of distributed SAR. In summary, due to model mismatch and computational burden, existing trajectory reconstruction methods cannot be directly applied to distributed SAR error reconstruction and compensation, and there is a lack of accurate modeling and reconstruction methods for multi-source distributed trajectory errors in distributed SAR. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention provides a high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction.

[0007] The technical solution adopted in this invention is: a high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction, the specific steps of which are as follows:

[0008] S1. Initialize the parameters of the distributed synthetic aperture radar system, then perform echo acquisition and demodulation to obtain the distributed SAR echo signal;

[0009] S2. Perform range pulse compression on the distributed SAR echo signal obtained in step S1, then perform autofocus on each platform to reconstruct the high-order trajectory error of the distributed SAR, and use the trajectory reconstruction results to obtain the focused sub-image of each platform.

[0010] S3. Based on step S2, perform platform-by-platform low-order trajectory error reconstruction to obtain the final distributed SAR focusing imaging result.

[0011] Furthermore, step S1 is specifically as follows:

[0012] The system parameters include: radar pulse parameters and configuration parameters.

[0013] Among them, the radar pulse parameters include: transmitted signal bandwidth B r Distance oversampling coefficient γ r Pulse width T r Observation time T s Transmitted signal carrier frequency f0, azimuth sampling frequency PRF, range sampling frequency f s Azimuth time η, range time τ, number of azimuth sampling points K, number of range sampling points N r Configuration parameters include: number of radar platforms N, platform speed v.

[0014] Set at azimuth time η, the nth radar platform UAV n The ideal trajectory center position is denoted as P. n (x n y n The actual trajectory center position P′ n (x n ′,y n ′); In azimuth time η, by UAV n The ideal distance to the target is historically denoted as R. n,η The actual distance from the historical record is denoted as R′. n,η .

[0015] Then, the echoes from each radar platform are recorded and demodulated to baseband to obtain the distributed SAR echo signal S. n(τ,η), expressed as follows:

[0016]

[0017] Where n represents the radar platform number, and n = 1, 2, 3, ..., N, T represents the target coordinates, and σ T Let T represent the scattering intensity of the target T, and c represent the speed of light. The expression representing the history of two-way distance is as follows:

[0018]

[0019] Among them, P n,η This represents the position of the nth radar platform in the azimuth direction with time η, and the range direction time variable τ = [-N]. r / 2:N r / 2] / f s The range of the azimuth time variable is η = [-T s / 2:T s / 2].

[0020] Furthermore, step S2 is specifically as follows:

[0021] S21, Distance-directed pulse compression;

[0022] Echo S n (τ,η) and reference signal S ref (τ) is processed to obtain the result after range-directed pulse compression. The expression is as follows:

[0023]

[0024] S22, platform-by-platform autofocus, and reconstruction of higher-order trajectory errors;

[0025] First, the phase center space is determined. Then, the phase error is estimated using the back projection principle. The phase centers of the N platforms at different times are denoted as the phase center space. The phase center space of the nth platform is then denoted as Q. n =[q n,1 ,q n,2 ,…,q n,k ,…,q n,K ].

[0026] Where, q n,k Let k represent the k-th azimuth phase center of the nth radar platform, where k∈[1,K].

[0027] Then, the back projection algorithm is used to calculate the coherent accumulation component for each selected phase center, as follows:

[0028] A1. Initialize the backward projection imaging space;

[0029] The back projection imaging space, denoted as Ω, is generated using the ground plane and unit vectors perpendicular to the ground plane. The back projection imaging space is then divided into M pixel units, denoted as Ω = [1, 2, ..., M].

[0030] Wherein, the grid interval d Ω It should be smaller than the system's imaging resolution.

[0031] A2. Calculate the delay;

[0032] Calculate the phase center q n,k To the m-th pixel T in the back projection imaging space m (x m ,y m ,z m The distance history R(n,k,m) is used to calculate the corresponding time delay for pixel T. m q n,k The time delay of each phase center is The expression is as follows:

[0033]

[0034] Where m∈[1,M], η k This represents the azimuth sampling time corresponding to the k-th azimuth sampling point.

[0035] A3. Calculate the BP basis matrix for the partial phase centers formed by the coherent accumulation components;

[0036] By combining the range-compressed echo signal with backscattering projection characteristics, the q-th... n,k The coherent accumulation component from the phase center to the m-th pixel unit is q n,k The basis vectors of the back projection imaging algorithm are The partial basis matrix obtained from the basis vectors of the K back projection imaging algorithms is denoted as... B n The dimension is N r ×K.

[0037] in, Indicates the nth r The distance sampling time corresponding to each distance sampling point, n r ∈[1,N r ]. (.) T This indicates the matrix transpose.

[0038] A4. Coherently superimpose the basis vectors from the K phase centers to the target scene to obtain the coherent accumulation component of each selected phase center;

[0039] The coherent accumulation component is represented as Among them, z n Represents an N r Dimensional vector.

[0040] Based on steps A1-A4, a self-focusing algorithm based on maximum sharpness is then used to estimate the higher-order trajectory error of each platform, i.e., to estimate the phase error of the selected phase center, as follows:

[0041] B1. Initialize parameters;

[0042] Set the maximum number of iterations MAX, initialize the number of iterations i = 1, and initialize the estimated phase error Φ0 = [0,...,0] for K phase centers.

[0043] B2. Calculate image sharpness;

[0044] The phase error Φ estimated using i-1 iterations i-1 The coherent accumulation after phase compensation is calculated as follows:

[0045] in, N represents r 3D vector, elements This represents the accumulated value of the m-th pixel in the imaging space.

[0046] Then calculate the intensity of the m-th pixel: (.) * The conjugate calculation is used to obtain image sharpness:

[0047] in, Represents the image sharpness function, and

[0048] B3. Estimate the phase error of the selected phase center using the coordinate descent method. Initialize j = 1;

[0049] B4. To maximize image sharpness, estimate the j-th parameter of the phase error.

[0050] The results of step B2 yield two M-dimensional vectors, u and v, expressed as follows:

[0051] u = zb j v=b j

[0052] Where v and b jLet represent the back projection result corresponding to the current j-th phase selection center, u represent the sum of the back projection results of all phase selection centers other than the j-th phase selection center, and z represent the sum of the back projection results of all phase selection centers.

[0053] Based on the BP principle, z = u + ve jφ φ represents the phase error term to be estimated. For the m-th pixel in the backprojection result, i.e., the m-th element u in the M-dimensional vector u, v, z m v m z m , there is z m =u m +v m e jφ Then the intensity f of the m-th pixel m The expression is as follows:

[0054]

[0055] (f con ) m =|u m | 2 +|v m | 2

[0056] (f Φ ) m =2Re(u m v m * e jΦ )

[0057] (f φ ) m =2Re(v m * u m )cosφ+2Im(v m * u m sinφ

[0058] Among them, (f con ) m and (f Φ ) m They represent f respectively m The constant part and the variable part, Im and Re respectively represent operations on imaginary and real numbers.

[0059] The image intensity can then be expressed as f = f con +f φ The maximum sharpness principle will be used to calculate maxψ(f) = ||f|| 2 The problem is transformed into finding the longest vector f in M-dimensional space. Φ=acosΦbsinΦ, which is the equation of an M-dimensional elliptic curve.

[0060] in, f represents the vector from the center of the ellipse to a point on the ellipse. con Let a = [a1, a2, ..., a] represent the vector from the origin to the center of the elliptical plane. M ],b=[b1,b2,...,b M ], and a m =2Re(v m * u m ), b m =2Im(v m * u m ).

[0061] Then, let u0 be the foot of the perpendicular from f to the two-dimensional plane spanned by a and b. The problem of finding the longest f is transformed into finding the farthest distance from a point on the ellipse to the point u0 outside the ellipse. For the entire image a and b, vectors a and b are orthogonalized to unit, and QR decomposition is used to obtain E = [e1 e2], generating the two-dimensional plane spanned by a and b, where e1 and e2 are expressed as follows:

[0062]

[0063] In the new coordinate system, a and b are represented as The expression is as follows:

[0064]

[0065] in, These represent the components of a and b transformed into the new coordinate system, respectively.

[0066] Then the expression for u0 in the new coordinate system e1, e2 is as follows:

[0067] u0 = -E T f con

[0068] The equation of the ellipse in the new coordinate system e1, e2 is expressed as follows: Rewritten as a quadratic form: f(h) = h T Rh = 1.

[0069] Using the properties of an ellipse, f con to a point on the ellipse If there is a longest f, then The expression satisfies the following:

[0070]

[0071] Where α represents the unknown parameter, R represents the symmetric positive definite matrix, and I represents the identity matrix, the expression is as follows:

[0072]

[0073] Then perform eigenvalue decomposition on R: R = VΔV T , where V represents the eigenvector matrix, Δ=diag{λ1,λ2} represents the eigenvalue matrix, and λ1 and λ2 represent the eigenvalues ​​of matrix R.

[0074] The maximum sharpness problem shares the same roots with the following fourth-order polynomial. α is obtained by solving the following equation, as detailed below:

[0075]

[0076] in,

[0077]

[0078] γ3=-2λ1λ2(λ1+λ2)

[0079] γ4=-(λ1λ2) 2

[0080] Where, [β1 β2] T =V T u0, u0 = -E T f con Let α be the smallest real root of the equation. When the maximum sharpness problem obtains the optimal solution, u is the point on the ellipse farthest from u0. The expression is as follows:

[0081]

[0082] in, Represent the real roots of the equation. This represents the estimated error phase.

[0083] Will Substitute R Finally, the j-th estimated phase error Φ is obtained. j The expression is as follows:

[0084]

[0085] B5. Termination determination of parameter estimation in the i-th iteration;

[0086] If j = K, terminate the i-th parameter estimation and obtain the i-th estimated phase error Φ based on maximum sharpness. i =[Φ1,...,Φ K Alternatively, return to step B4 and let j = j + 1.

[0087] B6. Termination determination of maximum sharpness phase error estimation iteration;

[0088] If i = MAX, the estimation result is obtained. The iteration terminates; otherwise, i = i+1, return to step B2, and continue the (i+1)th iteration.

[0089] Based on steps B1-B6, the higher-order trajectory error of each platform is estimated, and finally the higher-order trajectory error of the nth platform is reconstructed. That is, by utilizing the characteristic that the higher-order error of the platform does not cross the distance gate, the higher-order trajectory error of the nth platform is inverted.

[0090] Among them, the trajectory reconstruction result with higher-order trajectory errors is: λ = c / f0 represents the signal wavelength.

[0091] S23. Based on the trajectory reconstruction results from step S22, obtain the focused sub-images IM1′, IM2′, ... IM′ for each platform. N .

[0092] Furthermore, step S3 is specifically as follows:

[0093] S31. Initialize parameters, initialize the loop count n = 1, that is, initialize the radar platform number n = 1;

[0094] S32. Image registration: Solve for the affine transformation matrix of the focused sub-image of the platform.

[0095] First, IM′ is constructed using the SAR-Harris matrix. n The scale space of IM1′ and IM1′, the SAR-Harris matrix expression is as follows:

[0096]

[0097] Where x and y represent the two-dimensional pixel indices of the image, The standard deviation is expressed as Gaussian kernel, G x,γ =log(R) 1,γ ) and G y,γ =log(R) 3,γ ) represent the horizontal and vertical gradients, respectively, γ represents the exponential weight used to calculate the local mean, and R i,γ This represents the ratio of the exponentially weighted mean operator along the i-th direction.

[0098] Then, the multi-scale SAR-Harris function R is constructed. SH (x,y,γ)=det(C SH (x,y,γ)-tr(C SH(x,y,γ)) 2 The SAR-Harris function is used to threshold and filter edges and low-contrast points, keypoints are detected and feature descriptors are generated; the NNDR method is used to find matching point pairs, and the IM is obtained through the matching point pairs. n The offset between IM1 and IM2 is represented by the affine transformation matrix. F 1,1 F 1,2 F 2,1 F 2,2 Represents linear transformations between images, including: rotation, scaling, and shearing, F. 1,3 F 2,3 This indicates the image offset.

[0099] Where F represents the affine transformation relationship between the two images. They represent IM1 and IM respectively. n The imaging location of the same target point T in the image.

[0100] S33. Based on step S32, select l points in image IM1′, denoted as... And calculate its value in image IM′ based on the affine transformation matrix. n The corresponding position in is represented as

[0101] Where l≥3.

[0102] S34, According to position and Calculate the trajectory P with error n ′n;

[0103] First, initialize the loop variable con = 1, based on the two points in IM1. And its role in IM n The corresponding position in Combination Solve for the image position of point j by setting up equations. The distance equation is then listed and simplified as follows:

[0104]

[0105] in, and The quantity is known.

[0106] Then, the intersection point of the two circles is defined using a distance equation. The two solutions to the equation are the center of the trajectory with error, i.e., the actual trajectory center position P′. n With the center point of the false trajectory Based on the scenario constraints, by removing false trajectory center points, the trajectory center P′ with error can be obtained. n .

[0107] Finally, determine whether con≥l-1 is satisfied. If not, let con=con+1 and repeat step S34. If yes, go to step S35.

[0108] S35. Solving for trajectory parameters;

[0109] First, the translation amount is calculated for the UAV. n The average value of the trajectory offset relative to UAV1 is calculated as follows:

[0110] Then solve for the rotation amount according to the formula. Solve UAV n trajectory rotation angle θ n .

[0111] in, λ represents the signal wavelength, and v represents the platform speed;

[0112] S36. Determination of termination of trajectory parameter solution loop;

[0113] Determine if n≥N. If not, let n=n+1, return to step S32, and continue the (n+1)th iteration. If yes, terminate the iteration and reconstruct the linear trajectory error.

[0114] Based on the reconstructed trajectory parameters, distributed platform trajectory reconstruction is achieved by reconstructing the linear trajectory error, as expressed below:

[0115]

[0116] Where, Δx n Δy n This represents the nth radar trajectory offset; the reconstructed trajectory parameters include: two-dimensional offset parameter Δx = [x′1, x′2, ..., x′...]. n ···,x′ N ], Δy=[y′1,y′2,···y′ n ···,y′ N ], rotation parameter θ = [θ1, θ2, θ3,···θ n ···,θ N ].

[0117] S37. Based on the trajectory reconstruction results of step S36, perform platform-by-platform re-imaging to obtain the focused sub-images IM1, IM2, ... IM for each platform. N Finally, the images from multiple platforms are coherently superimposed to obtain the final distributed SAR focusing imaging result;

[0118] The final expression for the distributed SAR focusing imaging result is as follows:

[0119]

[0120] The beneficial effects of this invention are as follows: The method of this invention first solves for the higher-order errors of each platform in the distributed system to ensure the focusing degree of the imaging results of each platform. Then, it solves for the linear error of each platform. Using the quantitative relationship between the linear trajectory error parameters of each platform and the image offset, the linear error parameters are inverted, and finally, the high-precision motion trajectory of the distributed platforms is reconstructed to compensate for motion errors, thus achieving distributed SAR focusing imaging. This invention solves the motion error compensation problem in distributed synthetic aperture radar imaging, enabling high-precision reconstruction of motion errors across multiple distributed SAR platforms. It provides a prerequisite for multi-platform collaborative focusing imaging and has the advantages of high imaging accuracy and wide applicability to various scenarios. Attached Figure Description

[0121] Figure 1 This is a flowchart of a high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction according to the present invention.

[0122] Figure 2 This is a schematic diagram of the geometric configuration of distributed SAR in an embodiment of the present invention.

[0123] Figure 3 This is a schematic diagram illustrating the error model in an embodiment of the present invention.

[0124] Figure 4 This is a geometric interpretation diagram of the process of solving the true position of UAVn in an embodiment of the present invention.

[0125] Figure 5 This is a diagram showing the motion trajectory reconstruction result in an embodiment of the present invention.

[0126] Figure 6 This is a radar imaging result diagram from an embodiment of the present invention. Detailed Implementation

[0127] The method of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0128] like Figure 1 The flowchart of a high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction according to the present invention is shown below. The specific steps are as follows:

[0129] S1. Initialize the parameters of the distributed synthetic aperture radar system, then perform echo acquisition and demodulation to obtain the distributed SAR echo signal;

[0130] S2. Perform range pulse compression on the distributed SAR echo signal obtained in step S1, then perform autofocus on each platform to reconstruct the high-order trajectory error of the distributed SAR, and use the trajectory reconstruction results to obtain the focused sub-image of each platform.

[0131] S3. Based on step S2, perform platform-by-platform low-order trajectory error reconstruction to obtain the final distributed SAR focusing imaging result.

[0132] In this embodiment, step S1 is specifically as follows:

[0133] The system parameters include: radar pulse parameters and configuration parameters.

[0134] Among them, the radar pulse parameters include: transmitted signal bandwidth B r Distance oversampling coefficient γ r Pulse width T r Observation time T s Transmitted signal carrier frequency f0, azimuth sampling frequency PRF, range sampling frequency f s Azimuth time η, range time τ, number of azimuth sampling points K, number of range sampling points N r Configuration parameters include: number of radar platforms N, platform speed v.

[0135] In this embodiment, the geometric configuration of the distributed SAR is as follows: Figure 2 As shown in Table 1, the basic parameters are as follows.

[0136] Table 1

[0137] <![CDATA[Range oversampling coefficient (γ r )]]> 1.20 <![CDATA[pulse width (T r )]]> 0.5μs <![CDATA[Observation time (T s )]]> 0.6s <![CDATA[Carrier frequency of transmitted signal (f0)]]> 9.5GHz Azimuth sampling frequency (PRF) 500Hz <![CDATA[Number of range sampling points N r > 4096 Number of radar platforms (N) 16 Platform speed (v) 86km / h

[0138] In this embodiment, the range sampling frequency f s =γ r B r =1.2GHz, the number of azimuth sampling points per platform is K=T s PRF=300; Sets the error for the distributed system platform. The error model is explained as follows: Figure 3 As shown.

[0139] Set at azimuth time η, the nth radar platform UAV n The ideal trajectory center position is denoted as P. n (x n y n The actual trajectory center position P n ′(x n ′,y n ′); In azimuth time η, by UAV n The ideal distance to the target is historically denoted as R. n,η The actual distance from the historical record is denoted as R′.n,η .

[0140] Then, the echoes from each radar platform are recorded and demodulated to baseband to obtain the distributed SAR echo signal S. n (τ,η), expressed as follows:

[0141]

[0142] Where n represents the radar platform number, and n = 1, 2, 3, ..., N, T represents the target coordinates, and σ T Let T represent the scattering intensity of the target T, and c represent the speed of light. The expression representing the history of two-way distance is as follows:

[0143]

[0144] Among them, P n,η This represents the position of the nth radar platform in the azimuth direction with time η, and the range direction time variable τ = [-N]. r / 2:N r / 2] / f s = [-2.05:2.05]μs, the range of the azimuth time variable is η = [-T s / 2:T s / 2]=[-0.3:0.3]s.

[0145] In this embodiment, step S2 is specifically as follows:

[0146] S21, Distance-directed pulse compression;

[0147] Echo S n (τ,η) and reference signal S ref (τ) is processed to obtain the result after range-directed pulse compression. The expression is as follows:

[0148]

[0149] S22, platform-by-platform autofocus, and reconstruction of higher-order trajectory errors;

[0150] First, the phase center space is determined. Using the back projection principle, the phase error is estimated. The phase centers of N=16 platforms at different times are denoted as the phase center space. Then, the phase center space of the nth platform is denoted as Q. n =[q n,1 ,q n,2 ,…,q n,k ,…,q n,K ].

[0151] Where, q n,kLet k represent the k-th azimuth phase center of the nth radar platform, where k∈[1,K].

[0152] Then, the back projection algorithm is used to calculate the coherent accumulation component for each selected phase center, as follows:

[0153] A1. Initialize the backward projection imaging space;

[0154] The back projection imaging space, denoted as Ω, is generated using the ground plane and unit vectors perpendicular to the ground plane. The back projection imaging space is then divided into a grid of M = 2048 * 2048 pixel units, denoted as Ω = [1, 2, ..., M].

[0155] In order to distinguish between two adjacent point targets, the grid interval d Ω It should be slightly smaller than the system's imaging resolution.

[0156] A2. Calculate the delay;

[0157] Calculate the phase center q n,k To the m-th pixel T in the back projection imaging space m (x m ,y m ,z m The distance history R(n,k,m) is used to calculate the corresponding time delay for pixel T. m q n,k The time delay of each phase center is The expression is as follows:

[0158]

[0159] Where m∈[1,M], η k This represents the azimuth sampling time corresponding to the k-th azimuth sampling point.

[0160] A3. Calculate the BP basis matrix for the partial phase centers formed by the coherent accumulation components;

[0161] By combining the range-compressed echo signal with backscattering projection characteristics, the q-th... n,k The coherent accumulation component from the phase center to the m-th pixel unit is q n,k The basis vectors of the back projection imaging algorithm are The partial basis matrix obtained from the basis vectors of the K back projection imaging algorithms is denoted as... B n The dimension is N r ×K.

[0162] in, Indicates the nth rThe distance sampling time corresponding to each distance sampling point, n r ∈[1,N r ]. (.) T This indicates the matrix transpose.

[0163] A4. Coherently superimpose the basis vectors from K=300 phase centers to the target scene to obtain the coherent accumulation component of each selected phase center;

[0164] The coherent accumulation component is represented as Among them, z n Represents an N r Dimensional vector.

[0165] Based on steps A1-A4, a self-focusing algorithm based on maximum sharpness is then used to estimate the higher-order trajectory error of each platform, i.e., to estimate the phase error of the selected phase center, as follows:

[0166] B1. Initialize parameters;

[0167] Set the maximum number of iterations MAX = 5, initialize the number of iterations i = 1, and initialize the estimated phase error Φ0 = [0,...,0] for the K phase centers.

[0168] B2. Calculate image sharpness;

[0169] The phase error Φ estimated using i-1 iterations i-1 The coherent accumulation after phase compensation is calculated as follows:

[0170] in, N represents r 3D vector, elements This represents the accumulated value of the m-th pixel in the imaging space.

[0171] Then calculate the intensity of the m-th pixel: (.) * The conjugate calculation is used to obtain image sharpness:

[0172] in, Represents the image sharpness function, and

[0173] B3. Estimate the phase error of the selected phase center using the coordinate descent method. Initialize j = 1;

[0174] B4. To maximize image sharpness, estimate the j-th parameter of the phase error.

[0175] The results of step B2 yield two M-dimensional vectors, u and v, expressed as follows:

[0176] u = zb j v=b j

[0177] Where v and b j Let represent the back projection result corresponding to the current j-th phase selection center, u represent the sum of the back projection results of all phase selection centers other than the j-th phase selection center, and z represent the sum of the back projection results of all phase selection centers.

[0178] Based on the BP principle, z = u + ve jφ φ represents the phase error term to be estimated. For the m-th pixel in the backprojection result, i.e., the m-th element u in the M-dimensional vector u, v, z m v m z m , there is z m =u m +v m e jφ Then the intensity f of the m-th pixel m The expression is as follows:

[0179]

[0180] (f con ) m =|u m | 2 +|v m | 2

[0181] (f Φ ) m =2Re(u m v m * e jΦ )

[0182] (f φ ) m =2Re(v m * u m )cosφ+2Im(v m * u m sinφ

[0183] Among them, (f con ) m and (f Φ ) m They represent f respectively m The constant part and the variable part, Im and Re respectively represent operations on imaginary and real numbers.

[0184] The image intensity can then be expressed as f = f con +f φ The maximum sharpness principle will be used to calculate maxψ(f) = ||f|| 2 The problem is transformed into finding the longest vector f in M-dimensional space. Φ =acosΦbsinΦ, which is the equation of an M-dimensional elliptic curve.

[0185] in, f represents the vector from the center of the ellipse to a point on the ellipse. con Let a = [a1, a2, ..., a] represent the vector from the origin to the center of the elliptical plane. M ],b=[b1,b2,...,b M ], and a m =2Re(v m * u m ), b m =2Im(v m * u m ).

[0186] Then, let u0 be the foot of the perpendicular from f to the two-dimensional plane spanned by a and b. The problem of finding the longest f is transformed into finding the farthest distance from a point on the ellipse to the point u0 outside the ellipse. For the entire image a and b, vectors a and b are orthogonalized to unit, and QR decomposition is used to obtain E = [e1 e2], generating the two-dimensional plane spanned by a and b, where e1 and e2 are expressed as follows:

[0187]

[0188] In the new coordinate system, a and b are represented as The expression is as follows:

[0189]

[0190] in, These represent the components of a and b transformed into the new coordinate system, respectively.

[0191] Then the expression for u0 in the new coordinate system e1, e2 is as follows:

[0192] u0 = -E T f con

[0193] The equation of the ellipse in the new coordinate system e1, e2 is expressed as follows: Rewritten as a quadratic form: f(h) = h T Rh = 1.

[0194] Using the properties of an ellipse, f con to a point on the ellipse If there is a longest f, then The expression satisfies the following:

[0195]

[0196] Where α represents the unknown parameter, R represents the symmetric positive definite matrix, and I represents the identity matrix, the expression is as follows:

[0197]

[0198] Then perform eigenvalue decomposition on R: R = VΔV T , where V represents the eigenvector matrix, Δ=diag{λ1,λ2} represents the eigenvalue matrix, and λ1 and λ2 represent the eigenvalues ​​of matrix R.

[0199] The maximum sharpness problem shares the same roots with the following fourth-order polynomial. α is obtained by solving the following equation, as detailed below:

[0200]

[0201] in,

[0202]

[0203] γ3=-2λ1λ2(λ1+λ2)

[0204] γ4=-(λ1λ2) 2

[0205] Where, [β1 β2] T =V T u0, u0 = -E T f con To obtain α, we solve the equation where α is the smallest real root. When the maximum sharpness problem yields the optimal solution, u is the point on the ellipse farthest from u0. The expression is as follows:

[0206]

[0207] in, Represent the real roots of the equation. This represents the estimated error phase.

[0208] Will Substitute R Finally, the j-th estimated phase error Φ is obtained. j The expression is as follows:

[0209]

[0210] B5. Termination determination of parameter estimation in the i-th iteration;

[0211] If j = K, terminate the i-th parameter estimation and obtain the i-th estimated phase error Φ based on maximum sharpness. i =[Φ1,...,Φ K Alternatively, return to step B4 and let j = j + 1.

[0212] B6. Termination determination of maximum sharpness phase error estimation iteration;

[0213] If i = MAX, the estimation result is obtained. The iteration terminates; otherwise, i = i+1, return to step B2, and continue the (i+1)th iteration.

[0214] Based on steps B1-B6, the higher-order trajectory error of each platform is estimated, and finally the higher-order trajectory error of the nth platform is reconstructed. That is, by utilizing the characteristic that the higher-order error of the platform does not cross the distance gate, the higher-order trajectory error of the nth platform is inverted.

[0215] Among them, the trajectory reconstruction result with higher-order trajectory errors is: λ = c / f0 represents the signal wavelength.

[0216] S23. Based on the trajectory reconstruction results from step S22, image the focused sub-images IM′1, IM′2, ... IM′ of each platform. N .

[0217] In this embodiment, step S3 is specifically as follows:

[0218] S31. Initialize parameters, initialize the loop count n = 1, that is, initialize the radar platform number n = 1;

[0219] S32. Image registration: Solve for the affine transformation matrix of the focused sub-image of the platform.

[0220] First, IM′ is constructed using the SAR-Harris matrix. n The scale space of IM1′ and IM1′, the SAR-Harris matrix expression is as follows:

[0221]

[0222] Where x and y represent the two-dimensional pixel indices of the image, The standard deviation is expressed as Gaussian kernel, G x,γ =log(R) 1,γ ) and G y,γ =log(R) 3,γ ) represent the horizontal and vertical gradients, respectively, γ represents the exponential weight used to calculate the local mean, and R i,γ This represents the ratio of the exponentially weighted mean operator along the i-th direction.

[0223] Then, the multi-scale SAR-Harris function R is constructed. SH (x,y,γ)=det(C SH (x,y,γ)-tr(C SH (x,y,γ)) 2 The SAR-Harris function is used to threshold and filter edges and low-contrast points, keypoints are detected and feature descriptors are generated; the NNDR (Nearest Neighbor Ratio) method is used to find matching point pairs, and the IM (Instantaneous Motion) is obtained through these matching point pairs. n The offset between IM1 and IM2 is represented by the affine transformation matrix. F 1,1 F 1,2 F 2,1 F 2,2 Represents linear transformations between images, including: rotation, scaling, and shearing, F. 1,3 F 2,3 This indicates the image offset.

[0224] Where F represents the affine transformation relationship between the two images. They represent IM1 and IM respectively. n The imaging location of the same target point T in the image.

[0225] S33. Based on step S32, select l = 5 points in image IM1′, represented as follows: And calculate its value in image IM′ based on the affine transformation matrix. n The corresponding position in is represented as

[0226] Where l≥3.

[0227] S34, According to position and Calculate the trajectory P′ with error n ;

[0228] First, initialize the loop variable con = 1, based on the two points in IM1. And its role in IM n The corresponding position in Combination Solve for the image position of point j by setting up equations. The distance equation is then listed and simplified as follows:

[0229]

[0230] in, and The quantity is known.

[0231] The process of solving the true position of UAVn in this embodiment is as follows: Figure 4 As shown, the intersection point of two circles is defined by a distance equation, and the two solutions to the equation are the centers of the trajectories with error (the two circles). and The intersection point P′ n ), that is, the actual trajectory center position P′ n With the center point of the false trajectory Based on the scenario constraints, false trajectory center points are removed, resulting in the trajectory center P′ with error. n That is, by solving for vectors The angle β between the P' and the Y-axis can be used to determine P'. n (x n ′,y n ′), where β=β0+β1, β0 is The angle between β1 and the Y-axis is a triangle. The interior angle of . Based on geometric relationships, the expression can be derived as follows:

[0232]

[0233] Figure 4 The two circles actually have two intersection points, therefore a false trajectory center point exists. Removing this interference point is relatively simple. On the one hand, in practical applications... and The distance between them is constrained by the size of the scene, typically ranging from tens to hundreds of meters, while the distance between different radar platforms and the detection scene can range from hundreds to thousands of meters. and Much larger On the one hand, and Setting along the azimuth direction makes and P′ n Located on both sides of the scene. On the other hand, considering P n With P′ n The distance between them is within the positioning accuracy range, approximately on the order of centimeters, much smaller than P′. n and The distance between them can be easily avoided by setting a threshold. Interference;

[0234] Finally, determine whether con≥l-1 is satisfied. If not, let con=con+1 and repeat step S34. If yes, go to step S35.

[0235] S35. Solving for trajectory parameters;

[0236] First, the translation amount is calculated for the UAV. n The average value of the trajectory offset relative to UAV1 is calculated as follows:

[0237] Then solve for the rotation amount according to the formula. Solve UAV n trajectory rotation angle θ n .

[0238] in, λ represents the signal wavelength, and v represents the platform speed;

[0239] S36. Determination of termination of trajectory parameter solution loop;

[0240] Determine if n≥N. If not, let n=n+1, return to step S32, and continue the (n+1)th iteration. If yes, terminate the iteration and reconstruct the linear trajectory error.

[0241] Based on the reconstructed trajectory parameters, distributed platform trajectory reconstruction is achieved by reconstructing the linear trajectory error, as expressed below:

[0242]

[0243] Where, Δx n Δy n This represents the nth radar trajectory offset; the reconstructed trajectory parameters include: two-dimensional offset parameter Δx = [x′1, x′2, ..., x′...]. n ···,x′ N ], Δy=[y′1,y′2,···y′ n ···,y′ N ], rotation parameter θ = [θ1, θ2, θ3,···θ n ···,θ N ].

[0244] The reconstruction result of this embodiment is as follows: Figure 5 As shown, Figure 5 High-precision trajectory reconstruction results from distributed SAR are presented, including the plotting of the ideal trajectory, the actual trajectory, and the reconstructed trajectory using the method proposed in this paper. The reconstructed trajectory matches the actual trajectory very well, indicating that the actual trajectory was reconstructed with high quality.

[0245] S37. Based on the trajectory reconstruction results of step S36, perform platform-by-platform re-imaging to obtain the focused sub-images IM1, IM2, ... IM for each platform. N Finally, the images from multiple platforms are coherently superimposed to obtain the final distributed SAR focusing imaging result;

[0246] The final expression for the distributed SAR focusing imaging result is as follows:

[0247]

[0248] The imaging results of this embodiment are as follows: Figure 6 As shown, the preset strong points A, B, and C were magnified, and the imaging results without error compensation were obtained from... Figure 6 As shown in (a), in distributed SAR, due to trajectory errors, the images of each platform exhibit significant offset and defocusing, and the trajectory errors of multiple platforms affect the coherent fusion of images from all receivers, resulting in severely out-of-focus imaging. However, by using the method of this invention for multi-platform error reconstruction, re-imaging, and re-coherent superposition, an accurate estimate of the multi-platform trajectory errors is obtained. Figure 6 As shown in (b), the imaging results of the method of the present invention have good focusing effect.

[0249] In summary, the method of this invention utilizes the offset information between the images of each platform to solve for the linear error of each platform during the platform-by-platform backward projection imaging process, thereby accurately reconstructing the trajectory error of distributed synthetic aperture radar and achieving high-precision error compensation. This method solves the motion error compensation problem in distributed synthetic aperture radar imaging, enabling high-precision reconstruction of motion errors across multiple platforms in distributed SAR, providing a prerequisite for multi-platform collaborative focusing imaging, and offering advantages such as high imaging accuracy and wide applicability to various scenarios.

[0250] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction, the specific steps of which are as follows: S1. Initialize the parameters of the distributed synthetic aperture radar system, then perform echo acquisition and demodulation to obtain the distributed SAR echo signal; S2. Perform range pulse compression on the distributed SAR echo signal obtained in step S1, then perform autofocus on each platform to reconstruct the high-order trajectory error of the distributed SAR, and use the trajectory reconstruction results to obtain the focused sub-image of each platform. S3. Based on step S2, perform platform-by-platform low-order trajectory error reconstruction to obtain the final distributed SAR focusing imaging result; Step S3 is as follows: S31. Initialize parameters, initialize the loop count n=1, that is, initialize the radar platform number n=1; S32. Image registration: Solve for the affine transformation matrix of the focused sub-image of the platform. First, construct the SAR-Harris matrix separately. and In the scale space, the SAR-Harris matrix expression is as follows: ； in, , These represent the two-dimensional pixel indices of the image. The standard deviation is expressed as Gaussian kernel, and These represent the horizontal and vertical gradients, respectively. This represents the exponential weight used to calculate the local mean. This represents the ratio of the exponentially weighted mean operator along the i-th direction; Then, a multi-scale SAR-Harris function is constructed. The SAR-Harris function is used to threshold and filter edges and low-contrast points, key points are detected and feature descriptors are generated; the NNDR method is used to find matching point pairs, and through the matching point pairs, features are obtained. and The offset between them is represented by the affine transformation matrix. , Linear transformations between images include: rotation, scaling, and cropping. Indicates image offset; in, This represents the affine transformation relationship between two images. , They represent and The imaging location of the same target point T in the image; S33. Based on step S32, in the image Select 1 point, represented as And calculate its representation in the image based on the affine transformation matrix. The corresponding position in is represented as ; Where, l≥3; S34, According to position and Calculate the trajectory with error ; First, initialize the loop variable. =1, according to Two points in , and its The corresponding position in , , combined , Solve for the image position of point j by setting up equations. , Then, the distance equation is written and simplified as follows: ； in, and The quantity is known; Then, the intersection point of the two circles is defined using a distance equation. The two solutions to the equation are the center of the trajectory with error, i.e., the actual location of the trajectory center. With the center point of the false trajectory Based on the scenario constraints, false trajectory center points are removed to obtain the trajectory center with errors. ; Finally, determine whether the condition is satisfied. ≥1, otherwise let = +1, repeat step S34, if so, go to step S35; S35. Solving for trajectory parameters; First, the translation amount is calculated. Compared to Calculate the average value of the trajectory offset: ; Then solve for the rotation amount according to the formula. Solve trajectory rotation angle ; in, , Indicates the signal wavelength. Indicates platform speed; S36. Determination of termination of trajectory parameter solution loop; Determine if it satisfies , Indicate the number of radar platforms; otherwise, set... Return to step S32 and continue to step S32. If the next iteration is positive, the iteration terminates and the linear trajectory error is reconstructed. Based on the reconstructed trajectory parameters, distributed platform trajectory reconstruction is achieved by reconstructing the linear trajectory error, as expressed below: ； in, This represents the nth radar trajectory offset; the reconstructed trajectory parameters include: two-dimensional offset parameters. , Rotation parameters ; S37. Based on the trajectory reconstruction results from step S36, perform platform-by-platform re-imaging to obtain the focused sub-image for each platform. Finally, the images from multiple platforms are coherently superimposed to obtain the final distributed SAR focusing imaging result; The final expression for the distributed SAR focusing imaging result is as follows: 。 2. The high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction according to claim 1, characterized in that, The specific steps of S1 are as follows: The system parameters include: radar pulse parameters and configuration parameters; The radar pulse parameters include: transmitted signal bandwidth. Distance oversampling coefficient Pulse width Observation time Transmitted signal carrier frequency Azimuth sampling frequency Distance sampling frequency Direction to time Distance to time Number of azimuth sampling points K, number of range sampling points Configuration parameters include: number of radar platforms Platform speed ; Set in direction and time The nth radar platform The ideal trajectory center position is denoted as The actual center position of the trajectory In terms of location and time ,Depend on The ideal distance to the target is recorded in history. The actual distance from historical records ; Then, the echoes from each radar platform are recorded and demodulated to baseband to obtain distributed SAR echo signals. The expression is as follows: ； Where n represents the radar platform number, and n = 1, 2, 3, ..., N. Indicates the target coordinates. Indicate target The scattering intensity, Represents the speed of light. The expression representing the history of two-way distance is as follows: ； in, Indicates the azimuth time of the nth radar platform Location, distance to time variable The range of the azimuth-time variable is .

3. The high-precision imaging method for distributed synthetic aperture radar based on trajectory reconstruction according to claim 2, characterized in that, Step S2 is as follows: S21, Distance-directed pulse compression; echo With reference signal The process is performed to obtain the result after range-directed pulse compression. The expression is as follows: ； S22, platform-by-platform autofocus, and reconstruction of higher-order trajectory errors; First, determine the phase center space. Then, using the back projection principle, estimate the phase error. Denote the phase centers of N platforms at different times as the phase center space. The phase center space of the nth platform is then denoted as... ; in, This represents the k-th azimuth phase center of the nth radar platform. ; Then, the back projection algorithm is used to calculate the coherent accumulation component for each selected phase center, as follows: A1. Initialize the backward projection imaging space; The back-projected imaging space is generated using the ground plane and unit vectors perpendicular to the ground plane, denoted as . The back projection imaging space is divided into a grid, consisting of M pixel units, denoted as M0. ; Among them, grid interval It should be smaller than the system's imaging resolution; A2. Calculate the delay; Calculate the phase center To the rear projection imaging space 1 pixel Distance from history The corresponding time delay is calculated by tracing the distance from the history, for each pixel. No. The time delay of each phase center is The expression is as follows: ； ； in, , This represents the azimuth sampling time corresponding to the k-th azimuth sampling point; A3. Calculate the BP basis matrix for the partial phase centers formed by the coherent accumulation components; By combining the range-compressed echo signal with backscattering projection characteristics, the first... From the phase center to the first The coherent accumulation component of each pixel unit is , No. The basis vectors of the back projection imaging algorithm are The partial basis matrix obtained from the basis vectors of the K back projection imaging algorithms is denoted as... , The dimension is ; in, Indicates the first The distance sampling time corresponding to each distance sampling point Indicates matrix transpose; A4. Coherently superimpose the basis vectors from the K phase centers to the target scene to obtain the coherent accumulation component of each selected phase center; The coherent accumulation component is represented as ,in, Represent a dimensional vector; Based on steps A1-A4, a self-focusing algorithm based on maximum sharpness is then used to estimate the higher-order trajectory error of each platform, i.e., to estimate the phase error of the selected phase center, as follows: B1. Initialize parameters; Set the maximum number of iterations (MAX) and initialize the number of iterations. Initialize the estimated phase error of K phase centers ; B2. Calculate image sharpness; use Phase error estimated in the next iteration The coherent accumulation after phase compensation is calculated as follows: ; in, express 3D vector, elements Represents the first in the imaging space The accumulated value of each pixel; Then calculate the first... Intensity of each pixel: , The conjugate calculation is used to obtain image sharpness: ; in, Represents the image sharpness function, and ; B3. Estimate the phase error of the selected phase center using the coordinate descent method. ,initialization ; B4. To maximize image sharpness, estimate the j-th parameter of the phase error. ); Two M-dimensional vectors are obtained from the calculation results of step B2. The expression is as follows: ； in, and This represents the backprojection result corresponding to the j-th phase selection center. This represents the sum of the backprojection results of all phase selection centers except the j-th phase selection center. This represents the summation of all phase-selected center backprojection results; Based on the BP principle , This represents the phase error term to be estimated, for the m-th pixel in the backprojection result, i.e., the M-dimensional vector. The m-th element , ,have Then the intensity of the m-th pixel The expression is as follows: ； ； ； ； in, and They represent The constant part and the variable part, and These represent operations involving imaginary and real numbers, respectively. The image intensity can then be expressed as The maximum sharpness principle will be used to calculate The problem is transformed into Find the longest vector in 3D space , Right now Equations of 12-dimensional elliptic curves; in, This represents the vector from the center of the ellipse to a point on the ellipse. This represents the vector from the origin to the center of the elliptical plane. , ,and , ; Then set for exist The foot of the perpendicular from Zhang Cheng's two-dimensional plane will be used to find the longest... The problem is transformed into finding the distance from a point on the ellipse to a point outside the ellipse. The furthest distance; for the entire image , will vector and Unity orthogonalization, obtained using QR decomposition ,generate Zhang Cheng's two-dimensional plane, in which , The expression is as follows: ； In the new coordinate system, Represented as , The expression is as follows: ； in, , , They represent Components transformed to the new coordinate system; but In the new coordinate system , The following expression is: ； Equation of an ellipse in a new coordinate system , The following expression is Rewritten as a quadratic expression: ; Using the properties of an ellipse, to a point on the ellipse The longest ,but The expression satisfies the following: ； in, Indicates an unknown parameter. Describes a symmetric positive definite matrix. The identity matrix is ​​represented by the following expression: ； ； Then to Perform eigenvalue decomposition: ,in Represents the eigenvector matrix, Represents the eigenvalue matrix. , Representation matrix eigenvalues; The maximum sharpness problem shares roots with the following fourth-order polynomial, which can be obtained by solving the following equation. The details are as follows: ； in, ； ； ； ； ； in, , Let be the smallest real root of the equation, when the maximum sharpness problem obtains the optimal solution. That is, the distance on the ellipse The expression is as follows: ； in, Represent the real roots of the equation. This represents the estimated error phase; Will , Substitution Finally, the number was obtained. One estimated phase error The expression is as follows: ； B5, No. Termination determination of parameter estimation in the next iteration; If j=K, terminate the first... The second parameter estimation yields the third... Sub-phase error estimation based on maximum sharpness Conversely, return to step B4 and let 1; B6. Termination determination of maximum sharpness phase error estimation iteration; like The estimation results were obtained. The iteration terminates; otherwise, Return to step B2 and continue. The next iteration; Based on steps B1-B6, the higher-order trajectory error of each platform is estimated, and finally the higher-order trajectory error of the nth platform is reconstructed. That is, by utilizing the characteristic that the higher-order error of the platform does not cross the distance gate, the higher-order trajectory error of the nth platform is inverted. Among them, the trajectory reconstruction result with higher-order trajectory errors is: , Indicates the signal wavelength; S23. Based on the trajectory reconstruction results from step S22, obtain the focused sub-images for each platform. .

Images