A DEM imaging method and device based on navigation satellite bistatic InSAR

Abstract

The present invention discloses a method and device for DEM imaging based on navigation satellite bistatic InSAR, which can improve the temporal correlation of PS points and enhance the registration accuracy of the system. The method includes: using the sliding least squares fitting method to fit the DEM in the local coordinate system; performing second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM; obtaining the smoothing kernels to be processed corresponding to the discontinuous units after fitting; for each smoothing kernel to be processed, generating a Gaussian smoothing kernel corresponding to the smoothing kernel to be processed and centered symmetrically transformed with the smoothing kernel to be processed as the reference; for each smoothing kernel to be processed, performing a weighted summation operation on it and the corresponding Gaussian smoothing kernel, thereby obtaining the fitting smoothing result of the discontinuous unit corresponding to each smoothing kernel to be processed; using the fitting smoothing results of all discontinuous units as the finally smoothed DEM, and using the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result.

Description

Technical Field

[0001] The invention relates to the technical field of bistatic synthetic aperture radar, and in particular to a DEM imaging method and device based on bistatic InSAR of a navigation satellite. Background Art

[0002] Interferometric Synthetic Aperture Radar (InSAR) can be used for elevation inversion. It requires the use of an external DEM (Digital Elevation Model) to estimate and compensate for the flat ground phase and terrain phase in the interference phase, thereby making the interference fringes sparser, which is conducive to subsequent phase unwrapping and high-precision elevation inversion for deformation monitoring.

[0003] Compared with traditional airborne, satellite-based or ground-based systems, the dual-base InSAR system based on navigation satellites has three main characteristics: 1) low signal-to-noise ratio of echoes; 2) low resolution in the range direction; 3) asymmetry of the system geometric configuration. These three characteristics lead to a great difference between its imaging results and those of traditional systems, which is mainly manifested in the combination of point spread functions (PSF) of strong scattering points, poor resolution, and lack of texture information. When the system is applied to deformation inversion, if the method adopted by the traditional system is used, that is, directly introducing an external DEM as the imaging plane for processing, if the scene elevation fluctuation is complex and the DEM accuracy is high, the resolution unit of the permanent scatterer (PS) in the imaging result will jump in amplitude and phase, distort the edge, and have a discontinuous shape, which is not conducive to the subsequent extraction and registration of PS points, and affects the final deformation inversion accuracy. Summary of the invention

[0004] In view of this, the present invention provides a DEM imaging method and device based on dual-base InSAR of navigation satellites, which can make the amplitude and phase of the PS point resolution unit continuous and retain its PSF edge theoretical characteristics, ensure its shape continuity, improve the time correlation of the PS point, and enhance the registration accuracy of the system when the target scene has complex elevation fluctuations and the external source DEM has high accuracy.

[0005] In order to solve the above technical problems, the present invention is implemented as follows.

[0006] A DEM imaging method based on a navigation satellite dual-base InSAR, the method comprising the following steps:

[0007] Step S1: Use the sliding least squares fitting method to fit the DEM in the local coordinate system;

[0008] Step S2: Perform second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM;

[0009] Step S3: For each fitted discontinuous unit, perform boundary expansion with the fitted discontinuous unit as the center, and the expanded result is used as the smoothing kernel to be processed corresponding to the fitted discontinuous unit;

[0010] Step S4: For each smoothing kernel to be processed, generate a Gaussian smoothing kernel corresponding to the smoothing kernel to be processed and obtained by performing a central symmetry transformation based on the smoothing kernel to be processed;

[0011] Step S5: For each smoothing kernel to be processed, perform a weighted summation operation on it and the corresponding Gaussian smoothing kernel, and then obtain the fitted smoothing result of the discontinuous unit corresponding to each smoothing kernel to be processed;

[0012] Step S6: Use the fitted smoothing results of all discontinuous units as the finally smoothed DEM, and use the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result.

[0013] Preferably, the Step S1: Fitting the DEM in the local coordinate system using the sliding least squares fitting method includes:

[0014] Step S11: N units of the DEM are denoted as F T =[f(x1,y1) f(x2,y2) … f(x N ,y N )],

[0015] where f(x i ,y i ) is the i-th unit of the DEM, i = 1, 2, …, N; x i is the abscissa of the i-th unit, and y i is the ordinate of the i-th unit;

[0016] Step S12: For each f(x i ,y i ) among the N units of the DEM, adopt second-order weighted least squares fitting, including: for the i-th unit f(x i ,y i ), let the basis function be b T =[1 x y x 2 xy y 2 , assuming the number of reference units for fitting this unit is m, and denote f T =[f(x1,y1) f(x2,y2) … f(x m ,y m)], the weights of each reference unit follow a Gaussian distribution, denoted as:

[0017]

[0018] Then, according to the weighted least squares calculation, the fitting result of the i-th unit f(x i , y i ) is:

[0019]

[0020] Wherein,

[0021]

[0022] Wherein, x i and y i are two variables of the basis function b i , w i (x j , y j ) is the weighting coefficient of the j-th reference unit, 1 ≤ j ≤ m; f T is the transpose of the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit, and f is the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit;

[0023] Step S13: Based on the fitting results of the N units of the DEM, obtain the global fitting result of the DEM

[0024] Preferably, the step S2 includes: performing a second-order difference processing on the global fitting result of the DEM, and using the mean value of the global fitting result of the DEM as a threshold; traversing the N units of the DEM, finding the units with the second-order difference result greater than the threshold, and determining them as the fitting discontinuous units;

[0025] The second-order difference processing is: finding the second-order partial derivative of the global fitting result of the DEM, taking the partial derivative with respect to the x direction or the y direction.

[0026] Preferably, in the step S4, the method for generating a Gaussian smoothing kernel corresponding to the to-be-smoothed kernel and centrosymmetrically transformed based on the to-be-smoothed kernel is: generating a matrix with the same size as the to-be-smoothed kernel, the values in the matrix follow a two-dimensional Gaussian distribution, and then, with the central element of the matrix as the center, performing a centrosymmetric transformation on all elements of the matrix to generate a matrix with the smallest central value, gradually increasing surrounding values, and following a Gaussian distribution, which is the required Gaussian smoothing kernel.

[0027] An DEM imaging device based on navigation satellite bistatic InSAR provided by the present invention includes:

[0028] Fitting module: Configured to fit the DEM in the local coordinate system using the sliding least squares fitting method;

[0029] Incoherent unit determination module: Configured to perform second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM;

[0030] Expansion module: Configured to perform boundary expansion for each fitted discontinuous unit with the fitted discontinuous unit as the center, and the expanded result is used as the smoothing kernel to be processed corresponding to the fitted discontinuous unit;

[0031] Gaussian smoothing kernel determination module: Configured to generate, for each smoothing kernel to be processed, a Gaussian smoothing kernel corresponding to the smoothing kernel to be processed and centered symmetrically transformed with the smoothing kernel to be processed as the reference;

[0032] Smoothing result generation module: Configured to perform weighted summation operations on each smoothing kernel to be processed and its corresponding Gaussian smoothing kernel, and then obtain the fitted smoothing result of the discontinuous unit corresponding to each smoothing kernel to be processed;

[0033] Imaging result determination module: Configured to use the fitted smoothing results of all discontinuous units as the finally smoothed DEM, and use the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result.

[0034] Preferably, the fitting module includes:

[0035] First fitting sub-module: Configured for N units of the DEM, denoted as F T =[f(x1,y1) f(x2,y2) … f(x N ,y N )],

[0036] where f(x i ,y i ) is the i-th unit of the DEM, i = 1, 2, …, N; x i is the abscissa of the i-th unit, and y i is the ordinate of the i-th unit;

[0037] Second fitting sub-module: Configured to perform second-order weighted least squares fitting for each of the N units of the DEM, including: for the i-th unit f(x i ,y i ), setting the basis function as b i ,y i ) and setting the basis function as b T =[1 x y x 2 xy y 2, assuming the number of reference units for fitting this unit is m, denoted as f T = [f(x1, y1) f(x2, y2) … f(x m , y m )], and the weights of each reference unit follow a Gaussian distribution, denoted as:

[0038]

[0039] Then, according to the weighted least squares calculation, the fitting result of the i-th unit f(x i , y i ) is:

[0040]

[0041] Among them,

[0042]

[0043] Among them, x i and y i are two variables of the basis function b i , w i (x j , y j ) is the weighting coefficient of the j-th reference unit, 1 ≤ j ≤ m; f T is the transpose of the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit, and f is the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit;

[0044] The third fitting sub-module: configured to obtain the global DEM fitting result based on the fitting results of N units of the DEM

[0045] Preferably, the discontinuous unit determination module includes: performing second-order difference processing on the global DEM fitting result and using the mean value of the global DEM fitting result as the threshold; traversing the N units of the DEM to find the units with second-order difference results greater than the threshold and determining them as fitting discontinuous units;

[0046] The second-order difference processing is: obtaining the second-order partial derivative of the global DEM fitting result, taking the partial derivative with respect to the x direction or the y direction.

[0047] Preferably, the Gaussian smoothing kernel determination module generates a Gaussian smoothing kernel corresponding to the kernel to be smoothed and centered symmetrically transformed based on the kernel to be smoothed in the following manner: generating a matrix with the same size as the kernel to be smoothed, where the values within the matrix follow a two-dimensional Gaussian distribution, and then centering on the central element of the matrix, performing a central symmetric transformation on all elements of the matrix to generate a matrix with the smallest central value, gradually increasing surrounding values, and following a Gaussian distribution, which is the required Gaussian smoothing kernel.

[0048] Beneficial effects: The present invention solves the problems of amplitude and phase jumps within the PS point resolution unit, edge distortion, and shape discontinuity caused by directly introducing exogenous high-precision DEM information when the navigation satellite bistatic InSAR system performs complex scene imaging, improves the temporal correlation of PS points, and is conducive to improving the accuracy of the system applied to deformation inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 Schematic flow chart of the DEM imaging method based on navigation satellite bistatic InSAR according to the embodiment of the present invention;

[0050] Figure 2 Schematic configuration diagram of the navigation satellite bistatic SAR system according to the embodiment of the present invention;

[0051] Figure 3 Two-dimensional schematic diagram of Gaussian symmetric transformation according to the embodiment of the present invention.

[0052] Figure 4 Schematic flow chart of the fitting edge discontinuity smoothing process according to the embodiment of the present invention.

[0053] Figure 5 Complex scene high-precision DEM imaging method of the navigation satellite bistatic InSAR system according to the embodiment of the present invention.

[0054] Figure 6 Direct imaging result of the simulation original DEM according to the embodiment of the present invention.

[0055] Figure 7 Original terrain in the measured scene according to the embodiment of the present invention.

[0056] Figure 8 Phase error amplitude fitting result according to the embodiment of the present invention.

[0057] Figure 9 DEM fitting result in the measured scene according to the embodiment of the present invention.

[0058] Figure 10 DEM fitting and smoothing result in the measured scene according to the embodiment of the present invention.

[0059] Figure 11 The original DEM imaging result under the measured scenario of the embodiment cited in the present invention.

[0060] Figure 12 The fitted DEM imaging result under the measured scenario of the embodiment cited in the present invention.

[0061] Figure 13 The DEM imaging result after fitting and smoothing under the measured scenario of the embodiment cited in the present invention. Specific implementation manners

[0062] The present invention will be described in detail below in conjunction with the accompanying drawings and embodiments.

[0063] The present invention provides a DEM imaging method and device based on navigation satellite bistatic InSAR. After performing sliding least squares second-order fitting on the DEM in the local coordinate system, second-order differences are made, and the mean value of the fitted DEM is used as a threshold. If the second-order difference result of each unit is greater than this threshold, it is determined that the unit is a discontinuous unit in the fitted DEM. Taking each fitted discontinuous unit as the center, combining the DEM coverage range and the proportion of the discontinuous unit, the number of extended units is determined, so as to obtain the rectangular smoothing kernel to be smoothed corresponding to all fitted discontinuous units. Taking the size of the smoothing kernel to be smoothed as the standard, a Gaussian smoothing kernel after central symmetry transformation corresponding to it is generated. Weighted summation operations are performed on all smoothing kernels to be smoothed and their corresponding Gaussian smoothing kernels, and the obtained result is used as the fitted smoothing result of the discontinuous unit corresponding to the smoothing kernel to be smoothed. Taking this result as the back-projection imaging plane for imaging, finally an imaging result with no jump in amplitude and phase, no distortion at the edge, and continuous shape within the PS point resolution unit is obtained.

[0064] The DEM imaging method based on navigation satellite bistatic InSAR provided by the present invention, as Figure 1-2 shown, includes the following steps:

[0065] Step S1: Fitting the DEM in the local coordinate system using the sliding least squares fitting method;

[0066] Step S2: Making second-order differences for each item of the fitted DEM to find all discontinuous units existing in the fitted DEM;

[0067] Step S3: For each fitted discontinuous unit, performing boundary expansion with the fitted discontinuous unit as the center, and the expanded result is used as the smoothing kernel to be smoothed corresponding to the fitted discontinuous unit;

[0068] Step S4: For each smoothing kernel to be smoothed, generating a Gaussian smoothing kernel corresponding to the smoothing kernel to be smoothed and after central symmetry transformation with the smoothing kernel to be smoothed as the benchmark;

[0069] Step S5: For each kernel to be smoothed, perform a weighted summation operation with the corresponding Gaussian smoothing kernel, thereby obtaining the fitting and smoothing result of the discontinuous unit corresponding to each kernel to be smoothed;

[0070] Step S6: Use the fitting and smoothing results of all discontinuous units as the finally smoothed DEM, and use the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result.

[0071] The said Step S1: Fitting the DEM in the local coordinate system using the sliding least squares fitting method, includes:

[0072] Step S11: The N units of the DEM are denoted as F T = [f(x1,y1) f(x2,y2) … f(x N ,y N )],

[0073] wherein, f(x i ,y i ) is the i-th unit of the DEM, i = 1, 2, …, N; x i is the abscissa of the i-th unit, and y i is the ordinate of the i-th unit;

[0074] Step S12: For each f(x i ,y i ) among the N units of the DEM, adopt the second-order weighted least squares fitting, including: For the i-th unit f(x i ,y i ), set the basis function as b T = [1 x y x 2 xy y 2 , assume that the number of reference units for fitting this unit is m, and denote f T = [f(x1,y1) f(x2,y2) … f(x m ,y m ), and the weights of each reference unit follow a Gaussian distribution, denoted as:

[0075]

[0076] Then, according to the weighted least squares calculation, the fitting result of the i-th unit f(x i ,y i ) is:

[0077]

[0078] wherein,

[0079]

[0080] Among them, x i and y i are two variables of the basis function b i ; w i (x j , y j ) is the weighting coefficient of the j-th reference unit, where 1 ≤ j ≤ m; f T is the transpose of the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit. f is the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit.

[0081] Step S13: Obtain the global DEM fitting result based on the fitting results of N units of the DEM

[0082] The global DEM fitting result obtained by the present invention can not only retain the distribution characteristics of the local DEM, but also reduce the strong jump of adjacent DEM values. It can achieve the fitting effect of high-precision DEM in complex scenarios while ensuring the accuracy.

[0083] However, since this method cannot consider the numerical distribution outside the reference points during local least squares fitting, when the variation range of the full-scene DEM is large, there will be a problem of discontinuous local fitting edges, resulting in discontinuous imaging results. Therefore, smoothing processing needs to be introduced.

[0084] The step S2 includes: performing second-order difference processing on the global DEM fitting result and using the mean value of the global DEM fitting result as the threshold; traversing the N units of the DEM to find the units whose second-order difference results are greater than the threshold, and determining them as the discontinuous fitting units. The second-order difference processing is: finding the second-order partial derivative of the global DEM fitting result, which can be the partial derivative in the x direction or the y direction, aiming to find the positions where the adjacent DEM changes greatly, that is, the discontinuous positions in the fitting result.

[0085] The step S3: Taking the original DEM resolution as the reference basis, with the discontinuous fitting unit after fitting as the center, expand to the surrounding to ensure the distance that can ensure the existence of continuous units within the obtained smoothing kernel to be processed, and use it as the smoothing kernel to be processed corresponding to the discontinuous fitting unit after fitting.

[0086] The step S4, as Figure 3 shown, when directly introducing a Gaussian smoothing kernel for weighted summation smoothing processing at the discontinuous points, since the central weight value of the conventional Gaussian smoothing kernel is the largest and the surrounding weight values show a Gaussian distribution, as Figure 3As shown on the left side, during actual processing, the smoothing effect is poor. Therefore, the Gaussian smoothing kernel is subjected to a central symmetry transformation, that is, the central weight value is the smallest, the edge weight value is the largest, and the middle part presents a symmetric Gaussian distribution. As shown on the right side of Figure 3, it is the required Gaussian smoothing kernel. The processed Gaussian smoothing kernel can fully consider the data around the discontinuous pixel points of the fitting edge, and the smoothing effect is better.

[0087] The method for generating a Gaussian smoothing kernel corresponding to the kernel to be smoothed and subjected to a central symmetry transformation based on the kernel to be smoothed is as follows: Generate a matrix with the same size as the kernel to be smoothed. The values in this matrix follow a two-dimensional Gaussian distribution. Then, with the central element of this matrix as the center, perform a central symmetry transformation on all elements of this matrix to generate a matrix with the smallest central value, gradually increasing surrounding values, and following a Gaussian distribution, which is the required Gaussian smoothing kernel.

[0088] In this embodiment, the overall process of the discontinuous smoothing process of the fitting edge is as Figure 4 shown.

[0089] In step S6, based on the obtained one or more fitting smoothing results, use the one or more fitting smoothing results as the back-projection imaging plane for imaging, and finally obtain an imaging result with no jumps in amplitude and phase inside the PS point resolution unit, no edge distortion, and continuous shape.

[0090] In this embodiment, the main process of the high-precision DEM imaging algorithm in a complex scene is as Figure 5 shown.

[0091] In this embodiment, according to the parameters in Table 1, a slope terrain of 450 meters × 450 meters is used as the imaging scene, and 9 points are selected as target points. As Figure 6 shown, perform simulation on it, and the obtained results are as Figure 6 shown.

[0092] Table 1 Simulation system parameters

[0093]

[0094] Figure 6 In it, the edges of the echo PS point resolution units of the 9 target points are severely distorted, and there are also discontinuous problems. Using the method described in the present invention for imaging, the results are as Figure 7 shown. Figure 7 In it, there are no longer edge distortion problems in the resolution units of each target point, and the inside is also continuous.

[0095] Imaging was performed using the orbit data of Beidou MEO satellites on November 22, November 29, December 6, and December 13 to obtain 3 groups of registration pairs for the last 3 days and the first day. The imaging results in the above two cases were used to extract PS points, and the correlation coefficients of the PS points were calculated respectively and averaged. The obtained results are shown in Table 2.

[0096] Table 2 Comparison of Correlation Coefficients of PS Point Pairs in Simulation Results

[0097]

[0098] The results in Table 2 show that after fitting the high-precision DEM, the correlation of the image pairs for PS points will be improved, thereby improving the subsequent deformation inversion accuracy.

[0099] The processing results of the measured data are described below. A certain side terrain of the Malanzhuang mining area was selected as the imaging scene, as Figure 8 shown. The result of using sliding least squares fitting is as Figure 9 shown. There is a problem of discontinuous fitting edges in the area pointed by the arrow. Therefore, a Gaussian smoothing kernel after central symmetry transformation was used to smooth the DEM fitting result, and the processing result is as Figure 10 shown. Figure 9 And Figure 10 After comparison, it can be seen that the problem of discontinuous fitting has been solved.

[0100] Figure 11 、 Figure 12 、 Figure 13 are the imaging results of the original DEM, the fitted DEM, and the DEM after fitting and smoothing processing respectively. Among them, Figure 11 There are a total of 10 points to be evaluated in the small box, which are used for subsequent correlation coefficient calculation. Comparing the left and right frames of Figure 12 and Figure 13 respectively, it can be seen that unsmoothed processing will bring problems of discontinuous imaging; comparing the middle frames, it can be seen that the edge distortion problem of the PS point resolution unit in the imaging result of the processed DEM has been alleviated. Imaging was performed using the orbit data of Beidou MEO satellites on November 22, November 29, December 6, and December 13 to obtain 3 groups of registration pairs for the last 3 days and the first day. The 10 PS points to be evaluated selected from the imaging results of the original DEM and the fitted and smoothed DEM were extracted, and the correlation coefficients of these 3 registration pairs were calculated respectively and averaged. The obtained results are shown in Table 3.

[0101]

[0102]

[0103] As can be seen from the results in Table 3, the correlation coefficient in the DEM imaging results after processing has increased compared with that before processing, which proves the effectiveness of the present invention.

[0104] The present invention also provides a DEM imaging device based on navigation satellite bistatic InSAR, and the device includes:

[0105] A fitting module: configured to fit the DEM in the local coordinate system using the sliding least squares fitting method;

[0106] A discontinuous unit determination module: configured to perform second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM;

[0107] An expansion module: configured to perform boundary expansion with each fitted discontinuous unit as the center, and the expanded result is used as the smoothing kernel to be corresponding to the fitted discontinuous unit;

[0108] A Gaussian smoothing kernel determination module: configured to generate, for each smoothing kernel to be, a Gaussian smoothing kernel corresponding to the smoothing kernel to be and centered symmetrically transformed based on the smoothing kernel to be;

[0109] A smoothing result generation module: configured to perform weighted summation operation on each smoothing kernel to be and its corresponding Gaussian smoothing kernel, so as to obtain the fitting smoothing result of the discontinuous unit corresponding to each smoothing kernel to be;

[0110] An imaging result determination module: configured to use the fitting smoothing results of all discontinuous units as the finally smoothed DEM, and use the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result.

[0111] The above specific embodiments only describe the design principle of the present invention. The shapes and names of the components in this description can be different and are not limited. Therefore, those skilled in the art of the present invention can modify or make equivalent replacements to the technical solutions recorded in the foregoing embodiments; and these modifications and replacements do not depart from the spirit and technical solutions of the present invention, and shall fall within the protection scope of the present invention.

Claims

1. A DEM imaging method based on navigation satellite bistatic InSAR, characterized in that, The method includes: Step S1: Fitting the DEM in the local coordinate system using the sliding least squares fitting method; Step S2: Performing second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM; Step S3: For each fitted discontinuous unit, performing boundary expansion with the fitted discontinuous unit as the center, and using the expanded result as the smoothing kernel to be corresponding to the fitted discontinuous unit; Step S4: For each smoothing kernel to be, generating a Gaussian smoothing kernel corresponding to the smoothing kernel to be and centrosymmetrically transformed based on the smoothing kernel to be; Step S5: For each smoothing kernel to be, performing a weighted summation operation on it and the corresponding Gaussian smoothing kernel, so as to obtain the fitting smoothing result of the discontinuous unit corresponding to each smoothing kernel to be; Step S6: Using the fitting smoothing results of all discontinuous units as the finally smoothed DEM, and using the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result; Step S2 includes: performing second-order difference processing on the global fitting result of the DEM, and using the mean value of the global fitting result of the DEM as the threshold; traversing N units of the DEM to find the units with the second-order difference result greater than the threshold, and determining them as the fitting discontinuous units; The second-order difference processing is: finding the second-order partial derivative of the global fitting result of the DEM, and taking the partial derivative in the x direction or the y direction.

2. The method according to claim 1, wherein Step S1: Fitting the DEM in the local coordinate system using the sliding least squares fitting method includes: Step S11: N cells of the DEM, denoted as F T = [f(x1,y1) f(x2,y2) … f(x N ,y N )] where f(x i , y i ) is the i-th unit of the DEM, i = 1, 2, …, N; x i is the abscissa of the i-th unit, and y i is the ordinate of the i-th unit; Step S12: For each of the N cells of the DEM, f(x i , y i ), second-order weighted least squares fitting is adopted, including: for the i-th cell f(x i , y i ), let the basis function be b T = [1 x y x 2 xy y 2 , assuming that the number of reference cells for fitting this cell is m, denote f T = [f(x1, y1) f(x2, y2) … f(x m , y m ), and the weights of each reference cell follow a Gaussian distribution, denoted as: Then, according to the weighted least squares calculation, the fitting result of the i-th unit f(x i , y i ) is as follows: Wherein, where, w i (x j , y j ) is the weighting coefficient of the j-th reference unit, 1 ≤ j ≤ m; f T is the transpose of the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit, and f is the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit; Step S13: Obtain the global fitting result of the DEM based on the fitting results of N units of the DEM 3. The method according to claim 2, characterized in that, In Step S4, the method for generating a Gaussian smoothing kernel corresponding to the smoothing kernel to be and centrosymmetrically transformed based on the smoothing kernel to be is: generating a matrix with the same size as the smoothing kernel to be, where the values in the matrix follow a two-dimensional Gaussian distribution, and then taking the central element of the matrix as the center, performing centrosymmetric transformation on all elements of the matrix to generate a matrix with the minimum central value, gradually increasing surrounding values, and following a Gaussian distribution, which is the required Gaussian smoothing kernel.

4. A DEM imaging device based on navigation satellite bistatic InSAR, characterized in that The device includes: A fitting module: configured to fit the DEM in the local coordinate system using the sliding least squares fitting method; A discontinuous unit determination module: configured to perform second-order differences on each item of the fitted DEM to find all discontinuous units existing in the fitted DEM; An expansion module: configured to, for each fitted discontinuous unit, perform boundary expansion with the fitted discontinuous unit as the center, and using the expanded result as the smoothing kernel to be corresponding to the fitted discontinuous unit; A Gaussian smoothing kernel determination module: configured to, for each smoothing kernel to be, generate a Gaussian smoothing kernel corresponding to the smoothing kernel to be and centrosymmetrically transformed based on the smoothing kernel to be; A smoothing result generation module: configured to, for each smoothing kernel to be, perform a weighted summation operation on it and the corresponding Gaussian smoothing kernel, so as to obtain the fitting smoothing result of the discontinuous unit corresponding to each smoothing kernel to be; Imaging result determination module: configured to use the fitting and smoothing results of all discontinuous units as the finally smoothed DEM, and use the finally smoothed DEM as the back-projection imaging plane for imaging to obtain the final imaging result; The incoherence unit determination module includes: performing second-order difference processing on the global fitting result of the DEM, and using the mean value of the global fitting result of the DEM as the threshold; traversing N units of the DEM to find the units with the second-order difference result greater than the threshold, and determining them as fitting discontinuous units; The second-order difference processing is: finding the second-order partial derivative of the global fitting result of the DEM, taking the partial derivative with respect to the x direction or the y direction.

5. The device according to claim 4, wherein The fitting module includes: The first fitting sub-module: configured with N units of the DEM, denoted as F T = [f(x1, y1) f(x2, y2) … f(x N , y N )] Among them, f(x i , y i ) is the i-th unit of the DEM, where i = 1, 2, …, N; x i is the abscissa of the i-th unit, and y i is the ordinate of the i-th unit; Second fitting sub-module: Configured to perform second-order weighted least squares fitting for each of the N cells of the DEM f(x i , y i ), including: for the i-th cell f(x i , y i ), let the basis function be b T = [1 x y x 2 xy y 2 , assuming the number of reference cells for fitting this cell is m, denote f T = [f(x1, y1) f(x2, y2) … f(x m , y m ), and the weights of each reference cell follow a Gaussian distribution, denoted as: Then, according to the weighted least squares calculation, the fitting result of the i-th unit f(x i ,y i ) is as follows: Wherein, where, w i (x j , y j ) is the weighting coefficient of the j-th reference unit, 1 ≤ j ≤ m; f T is the transpose of the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit, and f is the matrix composed of the true DEMs of all reference units participating in the fitting of the i-th unit; The third fitting sub-module: configured to obtain the global fitting result of the DEM based on the fitting results of N units of the DEM 6. The device according to claim 5, characterized in that The Gaussian smoothing kernel determination module, the method for generating a Gaussian smoothing kernel corresponding to the kernel to be smoothed and centrosymmetrically transformed based on the kernel to be smoothed is: generating a matrix with the same size as the kernel to be smoothed, the values in this matrix follow a two-dimensional Gaussian distribution, and then centering on the central element of this matrix, performing centrosymmetric transformation on all elements of this matrix to generate a matrix with the smallest central value, gradually increasing surrounding values, and following a Gaussian distribution, which is the required Gaussian smoothing kernel.

Images