Satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing

Abstract

The present invention discloses a satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing, comprising: performing distance compression on an echo signal to obtain heterodyne data; compensating the heterodyne data using a PGA algorithm, and forming a marker vector using the compensated data; cyclically shifting each distance unit in the heterodyne data, and adaptively windowing the cyclically shifted data; performing phase gradient estimation on the windowed data according to the marker vector to obtain a phase gradient value; constructing a characteristic matrix, a phase error estimation matrix, and a signal-to-noise ratio weighting matrix for each azimuth unit according to the phase gradient value; constructing a phase error compensation matrix and a discrete Fourier transform matrix according to the characteristic matrix, the phase error estimation matrix, and the signal-to-noise ratio weighting matrix; and compensating the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result. The present invention achieves highly intelligent, efficient, and high-precision compensation.

Description

Technical Field

[0001] The invention belongs to the technical field of radar signal processing, and in particular relates to a satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing. Background Art

[0002] Inverse Synthetic Aperture Ladar (ISAL) is an active super-resolution imaging technology. With the advantage of the short wavelength in the optical band, ISAL can achieve higher imaging resolution in a shorter imaging time, while the imaging quality is close to that of optical images, which is conducive to target analysis and identification. The use of space-based ISAL can not only avoid the influence of atmospheric turbulence on laser transmission, but also realize the imaging of high-orbit satellite targets at long distances. It is of great significance in many defense and scientific research fields such as reconnaissance and early warning, classification and identification, and military surveillance in space. However, the short wavelength of the laser is both an advantage and a disadvantage. Because the operating wavelength of ISAL is at the μm level, it is very sensitive to the vibration of the satellite target. Vibration at the μm level will cause the imaging result to be defocused. Therefore, the phase error caused by the vibration of the satellite target must be compensated to achieve high-resolution imaging.

[0003] There are two main types of existing solutions for compensating for vibration phase errors: the first type is based on radar system design. For example, the Institute of Optoelectronics Technology of the Chinese Academy of Sciences proposed in its patent application "A device and method for estimating vibration phase errors in inverse synthetic aperture lidar imaging" (application publication number CN111693995A) that a multi-channel signal receiving system is used to simultaneously receive echoes. The spacing and position of each receiver are obtained through optimization design, and the collected data is subjected to differential and averaging processing to extract and compensate for vibration phase errors. The second type is the self-focusing imaging algorithm based on echo data processing. The most typical and best focusing method is Phase Gradient Autofocus (PGA). Through steps such as cyclic shifting, windowing, phase gradient estimation and iterative compensation, it can compensate for the space-invariant phase error caused by vibration and has good robustness. In the paper "Phase Gradient Matrix Autofocus for ISAL Space-Time-Varied Phase Error Correction" (IEEE PHOTONICS TECHNOLOGY LETTERS, 2020, vol. 32, no. 6, pp. 353-356), Z. Song et al. also proposed a phase gradient matrix autofocus (PGMA) imaging method. This method models the space-variant phase error as a two-dimensional time-frequency matrix, combines the phase error information of all distance units, and uses the weighted least squares fitting method to fit the linear function of the phase error with respect to the azimuth frequency, thereby constructing a compensation matrix to compensate for the space-variant phase error.

[0004] However, the compensation scheme designed based on radar system has high system complexity and great difficulty in implementation, and cannot compensate for the space-varying phase error caused by satellite angular vibration, and is not suitable for space-based ISAL. The PGA algorithm based on echo data processing is limited by its principle of utilizing the redundancy of phase error in azimuth, and cannot compensate for the space-varying phase error caused by satellite angular vibration. The PGMA algorithm can compensate for the space-varying phase error caused by satellite vibration, but this method combines all distance units for estimation. The incorrectly estimated distance unit not only reduces the estimation accuracy, but also has a large number of distance units for data with a large sampling rate in the range direction, which will greatly reduce the iterative compensation efficiency of the algorithm. In addition, since this method requires the user to manually input the window length parameter, and the subsequent iteration is very dependent on the windowed data of the previous iteration, and it is difficult for the user to properly input reasonable parameters, the compensation effect of the algorithm is unstable, which is not conducive to the needs of fast imaging of space-based ISAL. Summary of the Invention

[0005] In order to solve the above problems existing in the prior art, the present invention provides a method for satellite target ISAL vibration phase compensation based on improved phase gradient autofocusing. The technical problem to be solved by the present invention is achieved through the following technical solutions:

[0006] An embodiment of the present invention provides a satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing, the method comprising:

[0007] The ISAL echo signal of the satellite target containing vibration is processed by range compression to obtain heterodyne data;

[0008] Compensating the heterodyne data using a PGA algorithm, and using the compensated data to form a label vector as prior information;

[0009] cyclically shifting each distance unit in the heterodyne data, and performing adaptive windowing processing on the cyclically shifted data;

[0010] According to the label vector, the phase gradient of the windowed data is estimated by using the maximum likelihood estimation method to obtain the phase gradient value corresponding to each distance unit;

[0011] Constructing a characteristic matrix, a phase error estimation matrix, and a signal-to-noise ratio weighting matrix for each azimuth unit according to the phase gradient value, and constructing a phase error compensation matrix and a discrete Fourier transform matrix according to the characteristic matrix, the phase error estimation matrix, and the signal-to-noise ratio weighting matrix;

[0012] The heterodyne data is compensated according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result, so as to complete the vibration phase compensation of the satellite target ISAL.

[0013] In one embodiment of the present invention, the heterodyne data is compensated using a PGA algorithm, and a label vector is formed using the compensated data as prior information, including:

[0014] Extracting the phase error of each distance unit in the heterodyne data using a PGA algorithm;

[0015] For the phase error of each distance unit, the execution steps include: compensating the heterodyne data according to the phase error of the distance unit, and calculating the entropy value of the compensated data;

[0016] The average entropy value is calculated based on all entropy values, and it is determined whether each entropy value is greater than the average entropy value. If it is, the mark of the corresponding distance unit is set to 1, otherwise the mark of the corresponding distance unit is set to 0 to form a mark vector as prior information.

[0017] In one embodiment of the present invention, cyclically shifting each range unit in the heterodyne data includes:

[0018] Performing azimuth FFT transformation on the heterodyne data;

[0019] The strongest scattering point of each distance unit in the FFT transformed data is selected and moved to the center position corresponding to each distance unit. The data shifted out from one side is shifted in from the other side to form a cyclic shift.

[0020] In one embodiment of the present invention, performing adaptive windowing processing on the cyclically shifted data includes:

[0021] Normalizing the amplitude of all azimuth data of each range unit in the cyclically shifted data to obtain an azimuth energy distribution curve;

[0022] performing a maximum value search on the azimuthal energy distribution curve to extract a plurality of azimuthal units from the azimuthal energy distribution curve;

[0023] Perform interpolation fitting processing based on the extracted several azimuth units to obtain a fitting curve;

[0024] Searching from both sides of the edge of the fitting curve toward the center, detecting the first point with a sharp change on both sides of the fitting curve, and calculating the distance between the two points as the window length;

[0025] Adaptively windowing the cyclically shifted data according to the window length.

[0026] In one embodiment of the present invention, according to the label vector, a phase gradient estimation is performed on the windowed data using a maximum likelihood estimation method to obtain a phase gradient value corresponding to each range unit, including:

[0027] Perform azimuth IFFT transformation on the windowed data;

[0028] According to the label vector, selecting the distance unit marked as 1 from the data after IFFT transformation;

[0029] The maximum likelihood estimation method is used to estimate the phase gradient of all distance cells marked as 1 to obtain the phase gradient value corresponding to each distance cell.

[0030] In one embodiment of the present invention, the characteristic matrix of each orientation unit is constructed, and the formula is expressed as follows:

[0031]

[0032] Where H represents the characteristic matrix, n ranges from 1 to N, and N represents the number of orientation units. Represents the phase gradient value of the nth azimuth unit of the first distance unit marked as 1 in the marker vector, Express Find the average value of all elements in , represents the phase gradient value of the nth azimuth unit of the M1th range unit marked as 1 in the marking vector, Express The average value of all elements in is calculated, and M1 represents the number of distance units marked as 1 in the label vector;

[0033] The phase error estimation matrix of each azimuth unit is constructed as follows:

[0034]

[0035] Where φ(n) represents the phase error estimation matrix of the nth azimuth unit, Indicates the phase error value corresponding to the nth azimuth unit of the first distance unit marked as 1 in the mark vector, Indicates the phase error value corresponding to the nth azimuth unit of the distance unit marked as 1 in the mark vector,

[0036] The signal-to-noise ratio weighted matrix of each azimuth unit is constructed as follows:

[0037]

[0038] Where W represents the signal-to-noise ratio weighting matrix, Represents the normalized weight corresponding to the first distance unit marked as 1 in the mark vector, Indicates the phase gradient value The variance of Indicates the phase gradient value The variance of represents the normalized weight corresponding to the M1th distance unit marked as 1 in the mark vector, Indicates the phase gradient value The variance of , diag represents the function that generates a diagonal matrix.

[0039] In one embodiment of the present invention, the constructed phase error compensation matrix is ​​expressed as follows:

[0040]

[0041] Among them, PE N represents the phase error compensation matrix, j represents the imaginary unit, N represents the number of azimuth units, represents the translational vibration error parameter of the first azimuth unit, represents the translational vibration error parameter of the second azimuth unit, represents the translational vibration error parameter of the Nth azimuth unit, represents the angular vibration error parameter of the first azimuth unit, represents the angular vibration error parameter of the second azimuth unit, represents the angular vibration error parameter of the Nth azimuth unit, n ranges from 1 to N, where N represents the number of azimuth units, and f() represents the azimuth frequency function.

[0042] The constructed discrete Fourier transform matrix is ​​expressed as follows:

[0043]

[0044] in,

[0045] In one embodiment of the present invention, the heterodyne data is compensated according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result, which is expressed as follows:

[0046]

[0047] Where I(m,f) represents the compensation result, m ranges from 1 to M, M represents the number of distance units, n ranges from 1 to N, N represents the number of azimuth units, f represents the azimuth frequency, S'(m,n) represents the heterodyne data, FFT N represents the discrete Fourier transform matrix, PE N represents the phase error compensation matrix, Represents the dot product operation.

[0048] In one embodiment of the present invention, after compensating the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix, the method further includes:

[0049] Determine whether the window length of the windowed data is less than the empirical value. If so, output the compensation result to complete the vibration phase compensation of the satellite target ISAL. Otherwise, perform azimuth IFFT on the compensation result as new heterodyne data, return each range unit in the heterodyne data, and perform an adaptive windowing processing step on the cyclically shifted data until the window length of the windowed data is less than the empirical value.

[0050] In one embodiment of the present invention, when the window length of the windowed data is greater than or equal to the empirical value, it is continued to determine whether the current window length is equal to the window length of the previous windowed data. If it is, the statistical number is increased by 1, and it is determined whether the statistical number reaches a preset statistical threshold. If it is reached, the current windowed data is updated to be equal to 0.8 times the window length of the previous windowed data. The compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data. If it is not reached, the compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data. If it is not equal, the statistical number is cleared to 0, the compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data.

[0051] Beneficial effects of the present invention:

[0052] The present invention proposes a satellite target ISAL vibration phase compensation method based on improved phase gradient self-focusing. Aiming at the problem that the existing technology is difficult to compensate for the space-varying phase error caused by satellite angular vibration, and is difficult to simultaneously ensure estimation accuracy and iteration efficiency, and cannot realize space-based ISAL fast imaging, a satellite target ISAL vibration phase compensation method based on phase gradient self-focusing based on prior information and adaptive windowing technology is proposed. The method comprehensively considers estimation accuracy and iteration efficiency, adopts adaptive windowing so that users do not need to manually input parameters, avoids interference of external input parameters on the algorithm, and makes the algorithm more robust. At the same time, the PGA algorithm is innovatively used to compensate the data to obtain prior information. Based on the prior information, the distance unit for subsequent algorithm estimation is selected. The method can compensate for the space-varying and non-space-varying phase errors caused by satellite vibration with high intelligence, high efficiency and high precision, and meets the needs of space-based ISAL fast imaging.

[0053] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 1 is a flow chart of a satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing provided by an embodiment of the present invention;

[0055] Figure 2 1 is a schematic diagram of an azimuth energy distribution curve after normalizing amplitude processing of cyclically shifted complex image data provided by an embodiment of the present invention;

[0056] Figure 3Schematic diagram of the fitting curve and window length after interpolation fitting processing provided by an embodiment of the present invention;

[0057] Figure 4(a) to Figure 4(c) 1 is a schematic diagram of a scattering point model of a satellite target ISAL provided by an embodiment of the present invention, and imaging results of an RD algorithm in the absence and presence of vibration;

[0058] Figure 5(a) to Figure 5(d) 1 is a schematic diagram of an imaging result obtained by compensating for a true value of a vibration phase error according to an embodiment of the present invention, and an imaging result obtained by compensating for a true value of a vibration phase error according to a PGA algorithm, a PGMA algorithm, and the method proposed by the present invention;

[0059] Figure 6(a) to Figure 6(b) Graph showing the relationship between window length and width, and image contrast as a function of the number of iterations for the PGMA algorithm provided in an embodiment of the present invention and the method proposed in the present invention;

[0060] Figure 7(a) to Figure 7(b) This is a schematic diagram of imaging results using only prior information and only adaptive windowing technology provided by an embodiment of the present invention;

[0061] Figures 8(a) to 8(c) The PGMA algorithm provided by the embodiment of the present invention and the method proposed by the present invention are respectively graphs showing the relationship between the image contrast using only prior information and the window length and width and the image contrast as a function of the number of iterations when using only the adaptive windowing technology. DETAILED DESCRIPTION

[0062] The present invention will be further described in detail below with reference to specific examples, but the embodiments of the present invention are not limited thereto.

[0063] See Figure 1 The embodiment of the present invention provides a satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing, which specifically includes the following steps:

[0064] S10. Perform range compression processing on the ISAL echo signal of the satellite target containing vibration to obtain heterodyne data.

[0065] The present embodiment assumes that envelope alignment and initial phase correction have been completed, that is, imaging is performed on a turntable model. When the satellite target ISAL contains vibration, the vibration phase will be introduced into the echo signal. The heterodyne data obtained after range compression is expressed as follows:

[0066]

[0067] Wherein, S'(m,n) represents the heterodyne data of echo data containing vibration after range compression, S(m,n) represents the heterodyne data of echo data without vibration after range compression, m ranges from 1 to M, M represents the number of range units, n ranges from 1 to N, N represents the number of azimuth units, and j represents the imaginary unit. represents the Doppler frequency, represents the angular vibration error parameter, Represents the translational vibration error parameter. Usually in formula (1), It is known that S(m,n), and unknown.

[0068] S20. Compensate the heterodyne data using the PGA algorithm, and use the compensated data to form a label vector as prior information.

[0069] The purpose of calculating the prior information is to select distance units that can be used to estimate the phase error. If too few distance units are selected, the statistical regularity of the clutter background will be lost, affecting the estimation accuracy. If too many distance units are selected, more estimation errors will be introduced, while increasing the amount of calculation and slowing down the iteration efficiency. The embodiment of the present invention provides an optional solution that innovatively uses the PGA algorithm to preprocess the heterodyne data to obtain prior information. Specifically, the PGA algorithm is used to compensate the heterodyne data, and the compensated data is used to form a label vector as the prior information, including:

[0070] The phase error of each range unit in the heterodyne data is extracted using the PGA algorithm;

[0071] For the phase error of each distance unit, the execution steps include: compensating heterodyne data according to the phase error of the distance unit, and calculating the entropy value of the compensated data;

[0072] The average entropy value is calculated based on all entropy values, and each entropy value is judged to be greater than the average entropy value. If it is, the mark of the corresponding distance unit is set to 1, indicating that the estimation is correct. Otherwise, the mark of the corresponding distance unit is set to 0, indicating that the estimation is wrong, so as to form a mark vector as prior information.

[0073] It can be seen that the embodiment of the present invention uses the traditional PGA algorithm to extract the phase error of each range unit and compensates the heterodyne data respectively, and forms a label vector as prior information based on the entropy value of the compensated data. Based on this prior information, the range unit that can be used for subsequent algorithm estimation is selected, which improves the iterative efficiency and error estimation accuracy of the subsequent algorithm. It can compensate for the space-varying and non-space-varying phase errors caused by satellite vibration with high intelligence, high efficiency and high precision, and meet the needs of space-based ISAL fast imaging.

[0074] S30: cyclically shift each range unit in the heterodyne data, and perform adaptive windowing processing on the cyclically shifted data.

[0075] Cyclic shifting shifts the data in all azimuths of each range bin in the heterodyne data, moving the strongest scattering point to the center. This shifts the strongest scattering point to zero Doppler frequency, which can reduce the subsequent estimated phase gradient error. An embodiment of the present invention provides an optional solution for cyclically shifting each range bin in the heterodyne data, including:

[0076] Perform azimuth FFT transformation on the heterodyne data; select the strongest scattering point of each range unit in the FFT-transformed data, and move the strongest scattering point to the center position corresponding to each range unit. The data shifted out from one side is shifted in from the other side to form a cyclic shift.

[0077] As can be seen, the embodiment of the present invention performs azimuth FFT on the heterodyne data to transform it into the image domain, selects the strongest scattering point in each range cell in the image domain, and cyclically shifts it to the center of the range cell. Data shifted out from one side needs to be shifted in from the other side and cannot be directly discarded.

[0078] Furthermore, in order to retain all information of the strongest scattering point, that is, the area with a high signal-to-noise ratio, while removing the area with a low signal-to-noise ratio, the embodiment of the present invention performs adaptive windowing processing on the cyclically shifted data, including:

[0079] Normalized amplitude processing is performed on all azimuth data of each distance unit in the cyclically shifted data to obtain an azimuth energy distribution curve; a maximum search is performed on the azimuth energy distribution curve to extract several azimuth units from the azimuth energy distribution curve; interpolation fitting processing is performed based on the extracted several azimuth units to obtain a fitting curve; a search is performed from both sides of the edge of the fitting curve to the center to detect the first point with a sharp change on both sides of the fitting curve, and the distance between the two points is calculated as the window length for windowing; adaptive windowing processing is performed on the cyclically shifted data based on the window length. Specifically:

[0080] First, the energy of all azimuth directions in each range of the cyclically shifted data is accumulated, that is, normalized amplitude processing is performed. The calculation formula is:

[0081]

[0082] Among them, I′(m,n) represents the data after cyclic shift, l(n) represents the data after normalized amplitude processing, that is, the azimuth energy distribution curve is obtained as follows Figure 2 shown.

[0083] Next, a method for estimating the adaptive window length is used to extract the waveform profile based on curve fitting. Specifically, a maximum search is performed on the azimuth energy distribution curve. For example, a maximum search can be performed once, preferably twice, to extract several azimuth units from the azimuth energy distribution curve. Based on the extracted azimuth units, a piecewise cubic Hermite interpolation fitting process is performed using the makima function in Matlab to obtain the fitting curve as shown below. Figure 3 As shown, Figure 3 The red dotted line in the figure shows the curve, and the red circle in the figure represents the extracted orientation unit used in the fitting. Search from the edge of the fitting curve to the center, detect the first point with drastic changes on both sides of the fitting curve, and calculate the Euclidean distance between the two points as the window length, that is, Figure 3 The straight-line distance is shown by the solid red line.

[0084] Finally, according to the window length shown in Figure 3, the cyclically shifted data is adaptively windowed to extract data within the window length range for subsequent estimation and compensation.

[0085] It can be seen that the embodiment of the present invention adopts a method of extracting adaptive windowing based on curve fitting waveform contour, specifically using the makima function to perform segmented cubic Hermite interpolation fitting on the normalized amplitude data after cyclic shift, and determines the window length of adaptive windowing according to the fitting curve, so that the algorithm can adaptively select a suitable window length, avoiding the interference of external input parameters on the algorithm, further improving the estimation accuracy and iteration efficiency, and making the algorithm more intelligent and efficient.

[0086] S40 , performing phase gradient estimation on the windowed data using a maximum likelihood estimation method according to the label vector to obtain a phase gradient value corresponding to each range unit.

[0087] In the embodiment of the present invention, based on the marker vector, the maximum likelihood estimation method is used to perform phase gradient estimation on the windowed data to obtain the phase gradient value corresponding to each distance unit, including:

[0088] Perform azimuth IFFT transformation on the windowed data; select the range bins marked as 1 from the IFFT transformed data according to the marker vector; perform phase gradient estimation on all range bins marked as 1 using the maximum likelihood estimation method to obtain the phase gradient value corresponding to each range bin. Specifically:

[0089] First, perform azimuth IFFT transformation on the windowed data to obtain the range Doppler domain data, which is recorded as K(m,n). K(m,n)=[K1(n),K2(n),…K m(n)]; then the maximum likelihood estimation method is used to estimate the phase gradient value of the distance unit marked as 1 in the marking vector, and the formula is expressed as:

[0090]

[0091] in, It represents the phase gradient value of the nth azimuth unit of the m1th distance unit marked as 1 in the marker vector. m1 ranges from 1 to M1. M1 represents the total number of distance units marked as 1 in the marker vector. Represents the nth range Doppler number domain data of the range unit with the m1th mark as 1 in the mark vector, It represents the n+1th range-Doppler data in the range unit with the m1th mark as 1 in the mark vector. * represents the conjugate operation, and arg[] represents the phase function.

[0092] S50. Construct a characteristic matrix, a phase error estimation matrix, and a signal-to-noise ratio weighting matrix for each azimuth unit according to the phase gradient value, and construct a phase error compensation matrix and a discrete Fourier transform matrix according to the characteristic matrix, the phase error estimation matrix, and the signal-to-noise ratio weighting matrix.

[0093] In the embodiment of the present invention, the phase error value of each range unit can be obtained by accumulating the phase gradient values, which is expressed as follows:

[0094]

[0095] in, Indicates the phase gradient value of the nth azimuth unit of the first range unit marked as 1 in the marker vector.

[0096] In the embodiment of the present invention, the characteristic matrix of each orientation unit is constructed, and the formula is expressed as follows:

[0097]

[0098] Where H represents the characteristic matrix, n ranges from 1 to N, and N represents the number of orientation units. Represents the phase gradient value of the nth azimuth unit of the first distance unit marked as 1 in the marker vector, Express Find the average value of all elements in , represents the phase gradient value of the nth azimuth unit of the M1th range unit marked as 1 in the marking vector, Express The average value of all elements in is calculated, and M1 represents the number of distance units marked as 1 in the label vector;

[0099] The phase error estimation matrix of each azimuth unit is constructed as follows:

[0100]

[0101] Where φ(n) represents the phase error estimation matrix of the nth azimuth unit, Indicates the phase error value corresponding to the nth azimuth unit of the first distance unit marked as 1 in the mark vector, Indicates the phase error value corresponding to the nth azimuth unit of the distance unit marked as 1 in the mark vector,

[0102] The signal-to-noise ratio weighted matrix of each azimuth unit is constructed as follows:

[0103]

[0104] Where W represents the signal-to-noise ratio weighting matrix, Represents the normalized weight corresponding to the first distance unit marked as 1 in the mark vector, Indicates the phase gradient value The variance of Indicates the phase gradient value The variance of represents the normalized weight corresponding to the M1th distance unit marked as 1 in the mark vector, Indicates the phase gradient value The variance of , diag represents the function that generates a diagonal matrix.

[0105] In the embodiment of the present invention, the constructed phase error compensation matrix is ​​expressed as follows:

[0106]

[0107] Among them, PE N represents the phase error compensation matrix, j represents the imaginary unit, N represents the number of azimuth units, represents the translational vibration error parameter of the first azimuth unit, represents the translational vibration error parameter of the second azimuth unit, represents the translational vibration error parameter of the Nth azimuth unit, represents the angular vibration error parameter of the first azimuth unit, represents the angular vibration error parameter of the second azimuth unit, represents the angular vibration error parameter of the Nth azimuth unit, n ranges from 1 to N, where N represents the number of azimuth units, and f() represents the azimuth frequency function.

[0108] The constructed discrete Fourier transform matrix is ​​expressed as follows:

[0109]

[0110] in,

[0111] S60. Compensate the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result, so as to complete the vibration phase compensation of the satellite target ISAL.

[0112] In the embodiment of the present invention, the phase error compensation matrix and the discrete Fourier transform matrix are point-multiplied, and then the matrix multiplication is performed with the heterodyne data S'(m,n) to complete the compensation of the vibration phase error. The compensation result is expressed as follows:

[0113]

[0114] Where I(m,f) represents the compensation result, m ranges from 1 to M, M represents the number of distance units, n ranges from 1 to N, N represents the number of azimuth units, f represents the azimuth frequency, S'(m,n) represents the heterodyne data, FFT N represents the discrete Fourier transform matrix, PE N represents the phase error compensation matrix, Represents the dot product operation.

[0115] Furthermore, after compensating the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix, the embodiment of the present invention further includes:

[0116] A determination is made as to whether the windowed window length is less than an empirical value. If so, a compensation result is output to complete vibration phase compensation for the satellite target ISAL. Otherwise, the compensation result is subjected to an azimuth IFFT as new heterodyne data. Each range unit in the heterodyne data is then cyclically shifted, and the cyclically shifted data is subjected to an adaptive windowing process until the windowed window length is less than the empirical value. Thus, the present invention employs iterative compensation. As iterations proceed, the ISAL image becomes increasingly focused, i.e., the vibration phase error decreases, and the windowed window length also decreases. When the window length is determined to be less than the empirical threshold, a well-focused ISAL image can be output.

[0117] Furthermore, to ensure better iterative convergence, in an embodiment of the present invention, when the window length of the windowed data is greater than or equal to the empirical value, it is further determined whether the current window length is equal to the window length of the previous windowed data. If it is, the statistical number is increased by 1, and it is determined whether the statistical number reaches a preset statistical threshold. If it reaches, the current windowed data is updated to be equal to 0.8 times the window length of the previous windowed data. The compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data. If it is not reached, the compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data. If it is not equal, the statistical number is cleared to 0, the compensation result is subjected to azimuth IFFT as new heterodyne data, and each distance unit in the heterodyne data is returned to be cyclically shifted, and the adaptive windowing processing step is performed on the cyclically shifted data. It can be seen that when the window length is small enough, the width of the fitting curve envelope is greater than the window length, that is, the calculated window length of the current window is equal to the window length of the previous window. If the empirical value is continued to be used, the convergence will stop. At this time, the window length should be set to be equal to the previous window length * 0.8 to continue convergence. Preferably, the number of times the current window length is equal to the previous window length is counted. After exceeding the preset statistical threshold, the window length should be set to be equal to the previous window length * 0.8 to continue convergence. To ensure the convergence of the window length, the starting position of the search window length in each iteration starts from the previous window length position.

[0118] In order to verify the effectiveness of the satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing provided by an embodiment of the present invention, the following experiments were conducted for verification.

[0119] 1. Simulation conditions

[0120] The present invention was simulated using MATLAB 2020b developed by Mathworks, Inc., USA, on an Intel(R) Core(TM) i5-8250U 1.6 GHz CPU, an AMD Radeon 535 GPU, and a Microsoft Windows 10 operating system. The compared methods included the traditional Range-Doppler Algorithm (RD), the PGA algorithm, and the PGMA algorithm.

[0121] 2. Simulation content

[0122] The parameters of the inverse synthetic aperture lidar in the simulation experiment of the present invention are shown in Table 1.

[0123] Table 1 ISAL imaging radar and target parameters

[0124] parameter Numerical Light wavelength 1550nm Laser modulation bandwidth 4GHz Pulse repetition frequency 50kHz Single pulse duration 0.2μs Synthetic aperture pulse number 512 Target distance 50km Target angular velocity 0.1deg / s

[0125] The simulation experiment of the present invention uses the method proposed in the present invention to perform vibration phase error compensation imaging on a satellite target turntable model with an equivalent rotation speed of 0.1° / s, and finally obtains a target imaging result with a good focusing effect. The simulation image generated during the experiment is as follows: Figure 4(a) to Figure 4(c) 、 Figure 5(a) to Figure 5(d) 、 Figure 6(a) to Figure 6(b) ;in,

[0126] Figure 4(a) is a scattering point model diagram of a satellite target in a simulation experiment. The horizontal axis X represents the azimuth direction, in meters, and the vertical axis Y represents the range direction, in meters. Figure 4(b) is an imaging result diagram of ISAL imaging of a satellite target using the RD algorithm when there is no vibration. Figure 4(c) is an imaging result diagram of ISAL imaging of a satellite target using the RD algorithm when there is vibration. As can be seen from Figure 4(c), the RD algorithm has no ability to compensate for vibration, the image is completely defocused, and the details of the scattering points are completely submerged. Vibration phase error compensation imaging must be performed. In Figures 4(b) and 4(c), the horizontal axis represents the azimuth sampling points, and the vertical axis represents the range sampling points.

[0127] Figure 5(a) shows the result of compensation using the true value of the vibration phase error, which is the best theoretical result after compensation. Figure 5(b) shows the result of imaging using the PGA algorithm. It can be seen from Figure 5(b) that the PGA algorithm has a certain compensation ability for the non-space-varying phase error caused by vibration, but has no compensation ability for the space-varying phase error caused by diagonal vibration. The center of the image is focused, but the edge defocus is still serious. Figure 5(c) shows the result of imaging using the PGMA algorithm. It can be seen from Figure 5(c) that the PGMA algorithm can compensate for the phase error caused by micro-vibration of the satellite target, but the edge of the image still has a certain degree of defocus. Figure 5(d) shows the result of imaging using the method proposed in this invention. It can be seen from Figure 5(d) that the method proposed in this invention can compensate for the phase error caused by micro-vibration of the satellite target with good effect. The center and edge of the image are well focused. This method is also the result closest to compensation using the true value.

[0128] In order to quantitatively evaluate the focusing effect of the imaging results, image entropy and image contrast are used as image focusing indicators as shown in Table 2. The smaller the image entropy, the higher the image focusing quality, and the larger the image contrast, the higher the image quality.

[0129] Table 2 Comparison of image entropy and image contrast of different imaging methods

[0130] method truth value RD PGA PGMA Method of the present invention Image entropy 10.69 11.17 11.08 10.97 10.79 Image contrast 5.08 3.17 3.52 3.82 4.80

[0131] From the results in Table 2, it can be seen that: for the image entropy index, the image entropy gradually decreases from the RD algorithm to the method proposed in the present invention; for the image contrast index, the image contrast gradually increases from the RD algorithm to the method proposed in the present invention; the image entropy value and contrast after correction using the method proposed in the present invention are very close to the entropy value and contrast compensated with the true value, indicating that the compensation effect of the vibration phase error by the method proposed in the present invention is higher than that of the existing algorithm.

[0132] Figure 6(a) shows the relationship between the window length and width of the existing PGMA algorithm and the method proposed in the present invention and the number of iterations. It can be seen from Figure 6(a) that the window length of the PGMA algorithm depends on the windowed data of the previous iteration, and the width of the subsequent window length is 80% of the previous window length. In the first eight iterations, the method proposed in the present invention uses adaptive windowing based on curve fitting. While the window length converges, it is not completely affected by the previous window length. The selection of each window length is more intelligent. Figure 6(b) shows the relationship between the image contrast of the existing PGMA algorithm and the method proposed in the present invention and the number of iterations. It can be seen from Figure 6(b) that the image contrast improvement curve of the existing PGMA algorithm is flat, while the method proposed in the present invention can achieve a higher focusing effect with fewer iterations, thereby improving image quality and the correction efficiency of vibration phase error, verifying the advanced nature of the method.

[0133] In order to further verify the advancement of the prior information and curve fitting adaptive windowing technology of the method proposed in this invention, only the prior information and only the adaptive windowing method are used to compare with the existing technology. The simulation results are as follows: Figure 7(a) to Figure 7(b) and Figures 8(a) to 8(c) ;in,

[0134] Figure 7(a) shows the imaging result using only prior information, with the vibration conditions unchanged. As can be seen from Figure 7(a), both the center and edges of the image achieve good focusing. Figure 7(b) shows the imaging result using only the adaptive windowing technique. While the solar panel on the right side of the image exhibits some defocus, the overall image quality still exceeds that of existing algorithms. The results were also evaluated using image entropy and image contrast metrics, as shown in Table 3.

[0135] Table 3 Comparison of image entropy and image contrast of different imaging methods

[0136] method PGMA Adaptive windowing Prior information Image entropy 10.97 10.92 10.85 Image contrast 3.76 4.01 4.52

[0137] The results in Table 3 show that, compared with the existing PGMA algorithm, the image entropy obtained when using adaptive windowing alone and prior information alone is lower, and the image contrast is higher, indicating that the image quality is superior to that of the existing PGMA algorithm, verifying the effectiveness of both methods. Furthermore, analysis of the marker values ​​used in the proposed method shows that the total number of distance cells in the simulation is 480, while the number of distance cells corresponding to the marker vector is only 156. This means that based on this prior information, only 32.5% of the original number of distance cells needs to be estimated in each iteration to complete the compensation. This not only improves the estimation accuracy but also significantly increases the iteration rate.

[0138] Figure 8(a) is a graph showing the relationship between the window length and width of the existing PGMA algorithm and the method proposed in the present invention when only adaptive windowing is used, relative to the number of iterations. It can be seen from Figure 8(a) that as the number of iterations increases, the method proposed in the present invention and the existing PGMA algorithm have consistent window length and width; Figure 8(b) is a graph showing the relationship between the contrast of the image using only adaptive windowing and the number of iterations in the existing PGMA algorithm and the method proposed in the present invention. It can be seen from Figure 8(b) that the image contrast of the method proposed in the present invention is always better than that of the existing PGMA algorithm; Figure 8(c) is a graph showing the relationship between the contrast of the image using only prior information and the number of iterations in the existing PGMA algorithm and the method proposed in the present invention. It can be seen from Figure 8(c) that the image contrast of the method proposed in the present invention is always better than that of the existing PGMA algorithm, and as the iterations proceed, it is eventually much higher than that of the existing PGMA algorithm, proving the advanced nature of the method proposed in the present invention.

[0139] In summary, the embodiment of the present invention proposes a satellite target ISAL vibration phase compensation method based on improved phase gradient self-focusing, which addresses the problem that the existing technology is difficult to compensate for the space-varying phase error caused by satellite angular vibration, and is difficult to simultaneously ensure estimation accuracy and iteration efficiency, and cannot achieve space-based ISAL rapid imaging. A satellite target ISAL vibration phase compensation method based on phase gradient self-focusing based on prior information and adaptive windowing technology is proposed. This method comprehensively considers estimation accuracy and iteration efficiency, and adopts an adaptive windowing method so that the user does not need to manually input parameters, avoids interference of external input parameters on the algorithm, and makes the algorithm more robust. At the same time, the PGA algorithm is innovatively used to compensate the data to obtain prior information. Based on the prior information, the distance unit for subsequent algorithm estimation is selected, which can compensate for the space-varying and non-space-varying phase errors caused by satellite vibration with high intelligence, high efficiency and high precision, and meet the needs of space-based ISAL rapid imaging.

[0140] In the description of the present invention, it should be understood that the terms "first" and "second" are used for descriptive purposes only and should not be understood to indicate or imply relative importance or implicitly specify the number of the technical features indicated. Therefore, a feature specified as "first" or "second" may explicitly or implicitly include one or more of the features. In the description of the present invention, "plurality" means two or more, unless otherwise specifically defined.

[0141] Although the present invention is described herein in conjunction with various embodiments, those skilled in the art may understand and implement other variations of the disclosed embodiments by reviewing the specification and accompanying drawings in the process of implementing the claimed invention. In the specification, the word "comprising" does not exclude other components or steps, and "a" or "an" does not exclude multiple components or steps. The fact that certain measures are described in different embodiments does not mean that these measures cannot be combined to produce good results.

[0142] The above is a further detailed description of the present invention in conjunction with specific preferred embodiments, and the specific implementation of the present invention should not be considered to be limited to these descriptions. For those skilled in the art to which the present invention belongs, several simple deductions or substitutions can be made without departing from the concept of the present invention, and all of these should be considered to fall within the scope of protection of the present invention.

Claims

1. A satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing, characterized in that: The method comprises: The ISAL echo signal of the satellite target containing vibration is processed by range compression to obtain heterodyne data; Compensating the heterodyne data using a PGA algorithm, and using the compensated data to form a label vector as prior information; cyclically shifting each distance unit in the heterodyne data, and performing adaptive windowing processing on the cyclically shifted data; According to the label vector, the phase gradient of the windowed data is estimated by using the maximum likelihood estimation method to obtain the phase gradient value corresponding to each distance unit; Constructing a characteristic matrix, a phase error estimation matrix, and a signal-to-noise ratio weighting matrix for each azimuth unit according to the phase gradient value, and constructing a phase error compensation matrix and a discrete Fourier transform matrix according to the characteristic matrix, the phase error estimation matrix, and the signal-to-noise ratio weighting matrix; Compensating the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result, so as to complete the vibration phase compensation of the satellite target ISAL; The method of compensating the heterodyne data using a PGA algorithm and forming a label vector as prior information using the compensated data includes: Extracting the phase error of each distance unit in the heterodyne data using a PGA algorithm; For the phase error of each distance unit, the execution steps include: compensating the heterodyne data according to the phase error of the distance unit, and calculating the entropy value of the compensated data; The average entropy value is calculated based on all entropy values, and it is determined whether each entropy value is greater than the average entropy value. If it is, the mark of the corresponding distance unit is set to 1, otherwise the mark of the corresponding distance unit is set to 0 to form a mark vector as prior information.

2. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: Circularly shifting each range unit in the heterodyne data comprises: Performing azimuth FFT transformation on the heterodyne data; The strongest scattering point of each distance unit in the FFT transformed data is selected and moved to the center position corresponding to each distance unit. The data shifted out from one side is shifted in from the other side to form a cyclic shift.

3. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: Adaptive windowing is performed on the cyclically shifted data, including: Normalizing the amplitude of all azimuth data of each range unit in the cyclically shifted data to obtain an azimuth energy distribution curve; performing a maximum value search on the azimuthal energy distribution curve to extract a plurality of azimuthal units from the azimuthal energy distribution curve; Perform interpolation fitting processing based on the extracted several azimuth units to obtain a fitting curve; Searching from both sides of the edge of the fitting curve toward the center, detecting the first point with a sharp change on both sides of the fitting curve, and calculating the distance between the two points as the window length; Adaptively windowing the cyclically shifted data according to the window length.

4. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: According to the label vector, the phase gradient of the windowed data is estimated using the maximum likelihood estimation method to obtain the phase gradient value corresponding to each distance unit, including: Perform azimuth IFFT transformation on the windowed data; According to the label vector, selecting the distance unit marked as 1 from the data after IFFT transformation; The maximum likelihood estimation method is used to estimate the phase gradient of all distance cells marked as 1 to obtain the phase gradient value corresponding to each distance cell.

5. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: The characteristic matrix of each orientation unit constructed is expressed as follows: ； in, represents the feature matrix, n The value is 1~ N , N Indicates the number of orientation units, The first distance unit marked as 1 in the marker vector is n The phase gradient value of the azimuth unit, Express Find the average value of all elements in , Indicates the first The first distance cell marked as 1 n The phase gradient value of the azimuth unit, Express Find the average value of all elements in , Indicates the number of distance cells marked as 1 in the label vector; The phase error estimation matrix of each azimuth unit is constructed as follows: ； in, Indicates the n The phase error estimation matrix of azimuth units, The first distance unit marked as 1 in the marker vector is n The phase error value corresponding to the azimuth unit is: , Indicates the first The first distance cell marked as 1 n The phase error value corresponding to the azimuth unit is: ; The signal-to-noise ratio weighted matrix of each azimuth unit is constructed as follows: ； in, represents the signal-to-noise ratio weighting matrix, Represents the normalized weight corresponding to the first distance unit marked as 1 in the mark vector, , Indicates the phase gradient value The variance of Indicates the phase gradient value The variance of Indicates the first The normalized weight corresponding to the distance unit marked as 1, , Indicates the phase gradient value The variance of Represents a function that generates a diagonal matrix.

6. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: The constructed phase error compensation matrix is ​​expressed as follows: ； in, represents the phase error compensation matrix, represents the imaginary unit, N Indicates the number of orientation units, represents the translational vibration error parameter of the first azimuth unit, represents the translational vibration error parameter of the second azimuth unit, Indicates the N The translational vibration error parameters of the azimuth unit, represents the angular vibration error parameter of the first azimuth unit, represents the angular vibration error parameter of the second azimuth unit, Indicates the N Angular vibration error parameters of azimuth units, , n The value is 1~ N , N Indicates the number of orientation units, represents the azimuth frequency function, , ; The constructed discrete Fourier transform matrix is ​​expressed as follows: ； in, .

7. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1, characterized in that: The heterodyne data is compensated according to the phase error compensation matrix and the discrete Fourier transform matrix to obtain a compensation result, which is expressed as follows: ； in, Indicates the compensation result, m The value is 1~ M , M represents the number of distance units, n The value is 1~ N , N Indicates the number of orientation units, represents the azimuth frequency, represents heterodyne data, represents the discrete Fourier transform matrix, represents the phase error compensation matrix, Represents the dot product operation.

8. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 1 is characterized in that: After compensating the heterodyne data according to the phase error compensation matrix and the discrete Fourier transform matrix, the method further includes: Determine whether the window length of the windowed data is less than the empirical value. If so, output the compensation result to complete the vibration phase compensation of the satellite target ISAL. Otherwise, perform azimuth IFFT on the compensation result as new heterodyne data, return each range unit in the heterodyne data, and perform an adaptive windowing processing step on the cyclically shifted data until the window length of the windowed data is less than the empirical value.

9. The satellite target ISAL vibration phase compensation method based on improved phase gradient autofocusing according to claim 8, characterized in that: When the window length of the windowed data is greater than or equal to the empirical value, continue to determine whether the current window length is equal to the window length of the previous windowed data. If so, add 1 to the statistical number of times, and determine whether the statistical number reaches the preset statistical threshold. If so, update the current windowed data to be equal to 0.8 times the window length of the previous windowed data. Perform azimuth IFFT on the compensation result as new heterodyne data, return to perform cyclic shift on each distance unit in the heterodyne data, and perform adaptive windowing processing on the cyclically shifted data. If not, perform azimuth IFFT on the compensation result as new heterodyne data, return to perform cyclic shift on each distance unit in the heterodyne data, and perform adaptive windowing processing on the cyclically shifted data. If not, clear the statistical number to 0, perform azimuth IFFT on the compensation result as new heterodyne data, return to perform cyclic shift on each distance unit in the heterodyne data, and perform adaptive windowing processing on the cyclically shifted data.

Images