Macrograph No. 1 distributed wheel formation insar satellite height inversion method and equipment

Abstract

This invention belongs to the field of synthetic aperture radar interferometry and discloses a method for elevation inversion of the Hongtu-1 distributed wheel formation InSAR satellite, which can effectively improve the elevation inversion accuracy of digital elevation models (DEMs) for complex terrain. The method first establishes an equivalent phase center, simplifying the original bistatic imaging geometric model of the auxiliary satellite into a self-transmitting and self-receiving single-baseline multi-pass InSAR geometric model. Then, based on the single-baseline multi-pass InSAR geometric model, primary and secondary images are registered, generating three primary and secondary image pairs interferograms with different baseline lengths. These are then combined with parameters such as orbital information for terrain de-flattening phase processing. The de-flattened interferograms undergo multi-baseline phase unwrapping to recover the absolute terrain phase. Finally, the unwrapped phase is converted into the target scene elevation, and after geocoding, a DEM of the target area is generated.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar interferometry, and particularly relates to a macro... Figure 1 Method and equipment for InSAR satellite elevation inversion of distributed wheel formations. Background Technology

[0002] Digital Elevation Models (DEMs), as a crucial component of fundamental geographic information data, are widely used in terrain analysis, visibility calculation, and other multidimensional spatial analyses. Currently, stereo matching of optical imagery, Laser Radar (LiDAR), and spaceborne Interferometric Synthetic Aperture Radar (InSAR) are the three main technologies for generating large-scale DEMs. Among them, the spaceborne InSAR method can rapidly and accurately acquire high-precision DEMs in all weather conditions, making it particularly suitable for stable DEM acquisition under complex meteorological and weather conditions such as cloudy and foggy conditions. This has significant application value for geological disaster monitoring and early warning.

[0003] Spaceborne InSAR methods perform interferometric analysis on SAR images of the same region with a certain parallax angle, and extract and reconstruct the surface elevation information based on the interferometric phase information. Traditional single-satellite InSAR systems are mainly divided into single-antenna repeating orbit interferometric systems and dual-antenna single-track interferometric systems, and the accuracy of the acquired DEM is limited by a variety of factors.

[0004] Macro Figure 1 The HT-1SAR satellite adopts a "wheel" formation configuration of one main and three auxiliary satellites for Earth observation, which has the advantages of a large number of baselines, high mapping accuracy, and high mapping efficiency. Research on HT-1SAR data processing methods is of great significance. Summary of the Invention

[0005] To address the limitation of existing DEM accuracy, this invention provides a macro... Figure 1 The distributed wheel formation InSAR satellite elevation inversion method and equipment can effectively improve the DEM elevation inversion accuracy of complex terrain.

[0006] According to one aspect of the present invention, a macro is provided. Figure 1 The distributed wheel formation InSAR satellite elevation inversion method includes:

[0007] By establishing an equivalent phase center, the geometric model of single-base multi-pass InSAR is corrected;

[0008] Based on the corrected monobase multi-pass InSAR geometric model, the master and slave images are registered, and three master and slave image pairs interferograms with different baseline lengths are generated and then the flat-ground phase processing is performed.

[0009] Multi-baseline phase unwrapping is performed on the interferogram after removing the flat terrain phase to restore the absolute terrain phase;

[0010] The unwrapped phase is converted into the elevation of the target scene, and after geocoding, a DEM of the target area is generated.

[0011] As a further technical solution, the auxiliary satellite imaging geometric model is corrected by establishing an equivalent phase center, including:

[0012] Based on the orbital node data of the main satellite, the flight trajectory of the main satellite during the imaging period is fitted using a polynomial with time as the independent variable.

[0013] Substituting the time corresponding to the auxiliary satellite orbit node into the polynomial, we obtain the position and velocity of the main satellite at the same time, and then obtain the equivalent phase center position corresponding to the auxiliary satellite orbit node.

[0014] The polynomial is used for fitting, and the obtained polynomial parameters are used to replace the auxiliary satellite orbital parameters to complete the correction of the auxiliary satellite geometric model.

[0015] As a further technical solution, the registration of primary and secondary images is based on the corrected secondary satellite imaging geometric model, including:

[0016] Coarse registration: By determining the positions of corresponding points in the primary and secondary satellite images, the row and column coordinate offsets of the secondary satellite image relative to the primary satellite image are obtained;

[0017] Fine registration: Control points are evenly distributed on the main satellite image. The search range of corresponding points in the auxiliary satellite image is determined based on the offset obtained from coarse registration. The precise position of corresponding points in the auxiliary satellite image is determined according to the principle of maximizing coherence coefficient. The coordinate mapping relationship between the main and auxiliary satellites is established. Finally, coordinate transformation and interpolation resampling are performed on the auxiliary satellite image.

[0018] As a further technical solution, after completing the registration of the primary and secondary images, it also includes:

[0019] Complex conjugate multiplication is performed to obtain a differential interferogram, and the interferogram is then processed to remove the flat phase.

[0020] As a further technical solution, multi-baseline phase unwrapping is performed on the interferogram after removing the flat terrain phase to recover the absolute terrain phase, including:

[0021] Construct the following objective function and solve for the optimal phase ambiguity number gradient:

[0022]

[0023] Where (s, s-1) represents a pair of adjacent pixels, which have two directions: azimuth and slant range. Indicates the entanglement phase gradient; Δk i (s, s-1) is the gradient of the phase ambiguity number;

[0024] After determining the optimal phase ambiguity gradient, the LP minimum norm optimization model is used as the objective function to perform single-baseline phase unwrapping on each interferogram. After obtaining the phase ambiguity number of all pixels, the absolute terrain phase is restored.

[0025] As a further technical solution, when the ground point P=1, the problem of solving the LP minimum norm optimization model is transformed into a linear minimization problem; when the ground point P=2, the problem of solving the LP minimum norm optimization model is transformed into a weighted least squares problem.

[0026] As a further technical solution, the unwrapped phase is converted into the elevation of the target scene, and after geocoding, a DEM of the target area is generated, including:

[0027] Obtain the ground elevation corresponding to each pixel in the radar coordinate system;

[0028] Based on the vertical baseline length corresponding to the interferogram and the ground height, the final elevation value is obtained;

[0029] The final elevation value is geocoded to obtain the DEM product.

[0030] According to one aspect of the present invention, a macro is provided. Figure 1 The distributed wheel formation InSAR satellite elevation inversion equipment includes:

[0031] The first main module is used to correct the single-base multi-pass InSAR geometric model by establishing an equivalent phase center;

[0032] The second main module is used to register master and slave images based on the corrected single-base multi-pass InSAR geometric model, generate three master and slave image pairs interferograms with different baseline lengths, and perform flat-ground phase processing.

[0033] The third main module is used to perform multi-baseline phase unwrapping on the interferogram after removing the flat phase, and to restore the absolute terrain phase.

[0034] The fourth main module is used to convert the unwrapped phase into the elevation of the target scene, and after geocoding, generate the DEM of the target area.

[0035] According to one aspect of the present invention, a macro is provided. Figure 1The distributed wheel formation InSAR satellite elevation inversion device includes a memory and a processor. The memory stores program instructions that are executed by the processor, which calls the program instructions to execute the macro. Figure 1 The steps of the InSAR satellite elevation inversion method for distributed wheel formation.

[0036] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided, the non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the macro. Figure 1 The steps of the InSAR satellite elevation inversion method for distributed wheel formation.

[0037] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0038] This invention discloses a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion method can effectively improve the elevation inversion accuracy of digital elevation models (DEMs) for complex terrain. This method first establishes an equivalent phase center, simplifying the original dual-baseline imaging geometric model with separate transmitter and receiver on the auxiliary satellite into a self-transmitting and self-receiving single-baseline multi-pass InSAR geometric model. Then, based on the single-baseline multi-pass InSAR geometric model, primary and auxiliary images are registered, generating three primary and auxiliary image pairs interferograms with different baseline lengths. These are then combined with orbital information and other parameters for terrain de-flattening phase processing. The de-flattened interferograms undergo multi-baseline phase unwrapping to restore the absolute terrain phase. Finally, the unwrapped phase is converted into the target scene elevation, and after geocoding, a DEM of the target area is generated. Attached Figure Description

[0039] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0040] Figure 1 Macros provided for embodiments of the present invention Figure 1 A flowchart of the InSAR satellite elevation inversion method for distributed wheel formations.

[0041] Figure 2 This is a schematic diagram illustrating the principle of equivalent phase center processing provided in an embodiment of the present invention.

[0042] Figure 3 The schematic diagram of the InSAR interferometric phase geometry model provided in the embodiment of the present invention.

[0043] Figure 4The diagram illustrates the principle of minimum cost flow phase unwrapping provided in this embodiment of the invention. Detailed Implementation

[0044] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0045] This invention provides a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion method first establishes an equivalent phase center, simplifying the original dual-baseline model of separate transmission and reception of the auxiliary satellite into a self-transmitting and self-receiving single-baseline multi-pass InSAR geometric model. Then, based on the single-baseline multi-pass InSAR geometric model, the master and auxiliary images are registered to generate three master and auxiliary image pairs interferograms with different baseline lengths. The terrain phase is then removed by combining parameters such as orbit information. The multi-baseline phase is unwrapped on the interferogram after terrain phase removal to restore the absolute terrain phase. Finally, the unwrapped phase is converted into the target scene elevation and geocoded to generate the target area DEM.

[0046] The known data required for embodiments of this invention include: HT-1SAR imagery (single-look complex data) and header files, and external DEMs (such as open-source SRTM DEMs). Figure 1 As shown, the process of this embodiment of the invention includes the following steps:

[0047] Step 1: Correction of auxiliary satellite imaging geometric model.

[0048] like Figure 2 As shown, in the dual-receiver mode of SAR satellite, the positional relationship between the auxiliary satellite (HT-1B / C / D) and the ground point during imaging can be described by the range-Doppler model, specifically as shown in equation (1):

[0049]

[0050] Where, p t Indicates the location of a point on the ground; p m,v m p represents the position and velocity vector of the ground point imaged by the main satellite at time t; s ,v s The vector represents the position and velocity of the auxiliary satellite; λ is the radar wavelength; f dc R1 and R2 represent the center frequency of the auxiliary satellite's Doppler signal; R1 and R2 represent the slant distances from the primary and auxiliary satellites to the ground point, respectively.

[0051] The midpoint of the line connecting the phase centers of the main and auxiliary radar antennas is used as the equivalent phase center of the auxiliary radar antenna, as shown in equation (2):

[0052]

[0053] Where, p se (t) represents the equivalent phase center of the auxiliary satellite; R 2se This represents the slant distance between the secondary satellite and the ground point target after the equivalent phase center. The HT-1 image header file provides precise satellite orbit data at 1-second intervals. The specific method for correcting the secondary satellite imaging geometry model is to fit the primary satellite's flight orbit during the imaging time period using a polynomial with time (using UTC time) as the independent variable, based on the primary satellite's orbital node data.

[0054]

[0055] Among them, (X) m ,Y m Z m ),(V mx V my V mz The values ​​) represent the orbital position and velocity of the primary satellite at time t; a, b, c, d, e, and f are polynomial parameters obtained by fitting orbital node data. The time t corresponding to the orbital nodes of the secondary satellite is then considered. s Substituting into equation (3), the position and velocity of the primary satellite at the same time are obtained. Equation (2) can then be used to obtain the equivalent phase center position corresponding to the orbital node of the secondary satellite. Furthermore, the polynomial shown in equation (3) is used for fitting, and the obtained polynomial parameters are used to replace the orbital parameters of the secondary satellite to complete the correction of the secondary satellite's geometric model.

[0056] After replacing the phase center of the auxiliary satellite with the equivalent phase center, the bistatic imaging model described by equation (1) can be approximated as the classical single-base imaging model, specifically as equation (4):

[0057]

[0058] Step 2: Registration of primary and secondary SAR images and generation of deflating interferograms.

[0059] The registration of primary and secondary SAR images consists of two steps: coarse registration and fine registration. Coarse registration determines the positions of corresponding points in both primary and secondary satellite images and obtains the row and column coordinate offsets of the secondary satellite image relative to the primary satellite image. Using the center point of the primary satellite image as a reference point, the latitude and longitude coordinates of the corresponding ground point P are obtained from the primary satellite header file. Assuming the geodetic height of P is 0, the coordinates are converted to spatial rectangular coordinates (X, Y, X) in WGS-84. P ,Y P Z P In the coarse registration step, it is assumed that the single-look complex image data used has been focused onto the zero Doppler plane, that is, the radar side-looking direction and the satellite flight direction are perpendicular. By minimizing equation (5), the imaging time of ground point P in the auxiliary satellite can be obtained.

[0060]

[0061] Among them, the imaging time and satellite position at time t are corresponding to each row of pixels of the auxiliary satellite. s (t),Y s (t),Z s (t) and flight velocity vector (v) X (t),v Y (t),v z (t) can be obtained directly from the header file or through interpolation.

[0062] After obtaining the imaging time t of the auxiliary satellite at point P and the corresponding satellite position, the pixel coordinates (r) of P in the auxiliary satellite image are determined by equation (6). s ,c s ):

[0063]

[0064] Among them, R near dR is the slant range corresponding to the first column of pixels of the auxiliary satellite; t0 is the initial imaging time of the auxiliary satellite image; dt is the azimuth sampling time interval. near The parameters dR, t0, and dt can all be obtained from the image header file. After determining the position of the corresponding point in the main and auxiliary satellite images, obtain the row and column coordinate offsets of the auxiliary satellite image relative to the main satellite image.

[0065] In the fine registration step, control points are first evenly distributed on the master satellite image. Based on the offset obtained from the coarse registration, the search range for corresponding points in the slave satellite image is determined. Following the principle of maximizing coherence, the precise location of these corresponding points in the slave satellite image is determined, establishing the coordinate mapping relationship between the master and slave satellites. Finally, coordinate transformation and interpolation resampling are performed on the slave satellite image. After master-slave image registration is completed, complex conjugate multiplication is performed to obtain a differential interferogram.

[0066] The schematic diagram of single-base repeating orbit InSAR interferometric phase principle is shown below. Figure 3 As shown, O represents the Earth's center. The main satellite radar incident angle θ performs side-looking imaging of ground point P, with a slant range of R1. M represents the position of the satellite antenna phase center during imaging, H represents the altitude of the main satellite antenna flight platform, P is the ground target, h is the ground elevation value of point P to be solved, and B represents the spatial baseline between the main and auxiliary satellites. The angle between baseline B and the horizontal direction is denoted as α. The parallel baseline B is parallel to the slant range MP of the main satellite radar. || and the vertical baseline B perpendicular to the radar slant range direction SP ⊥ The result is obtained by calculation using equation (7):

[0067]

[0068] Assuming no change in surface and atmospheric conditions during the imaging process of the primary and secondary satellite images, the interferometric phase is simplified to the phase φ of the reference ellipsoid. ref Terrain phase φ topo and noise phase φ noi Three parts. The noise phase is removed using the Goldstein phase filtering method, and the ellipsoidal phase is calculated using equation (8):

[0069]

[0070] Where P0 is the equidistant projection of point P onto the reference ellipsoid, i.e., r1 = R1, R2 represents the slant range of the auxiliary satellite, and the side angle from the main satellite radar to P0 is θ0. The phase φ of the reference ellipsoid is removed from the interferogram. ref Terrain phase φ can be obtained topo That is, the flattened interference diagram.

[0071] Step 3: Unwrap the multi-baseline phase to recover the absolute terrain phase.

[0072] For a set of HT-1 satellite data, using HT-1A data as the primary image and B / C / D data as secondary images respectively, registration and interferometry are performed to generate three flattened interferograms. The corresponding vertical baseline is B AB B AC B AD Original entangled phase of pixel s The following relationship exists between the true absolute phase ψ(s):

[0073]

[0074] Where k(s) represents the phase ambiguity number of pixel s. For pixel s and its corresponding neighboring pixel s-1 in the interferogram, the objective function shown in equation (10) is constructed:

[0075]

[0076] Where (s, s-1) represents a pair of adjacent pixels, which have two directions: azimuth and slant range. Indicates the entanglement phase gradient; Δk i (s, s-1) is the phase ambiguity number gradient, which takes any integer value, Δ i (s,s-1),Δ j (s, s-1) represents the pixel pair variable consisting of pixel s and its corresponding pixel s-1 in the row and column directions. To determine the optimal phase ambiguity gradient, a brute-force method is used: the ambiguity gradient Δk is set based on prior knowledge such as the external DEM. i The search range of (s, s-1) is [min] i (s),max i (s)], within a given range, all fuzzy number gradient combinations (Δk) are traversed using an exhaustive method. AB (s,s-1),Δk AC (s,s-1),Δk AD (s,s-1)), the combination that minimizes equation (10) will be taken as the optimal solution for the phase ambiguity number gradient. When setting the ambiguity search range, it is necessary to consider the range [min i (s),max i The unique optimal solution of equation (10) is found within (s). Based on the nonlinear mixed integer programming criterion and the measured vertical baseline length, the search range of the fuzzy number gradient is determined to be [-4, 4].

[0077] To obtain the optimal phase ambiguity number gradient Δk i After (s,s-1)(i∈{AB,AC,AD}), the constructed LP minimum norm optimization model performs single-baseline phase unwrapping on each interferogram, and the objective function is shown in Equation (11):

[0078]

[0079] Among them, w i (s,s-1)∈[0,1] are weighting coefficients, which can be determined by the interferometric quality map or other prior knowledge, k i (s) represents the phase blur number gradient of pixel s in interferogram i, and Z represents the integer domain. For the case of P = 1, Δ is defined. x k(i,j),Δ y k(i,j) is the gradient of the phase blur number of the pixel in the i-th row and j-th column in the row and column directions, which has been obtained by solving equation (10); Let be the discrete partial derivatives of the winding phase in the row and column directions. Then, equation (11) is transformed into the following nonlinear minimization problem:

[0080]

[0081] Where M and N represent the number of rows and columns of the interferogram, respectively, w x (i,j) represents the row-direction weight coefficient of the pixel in the i-th row and j-th column, w y (i,j) represents the column direction weight coefficient of the pixel in the i-th row and j-th column.

[0082] Since the true terrain phase should belong to an irrotational vector field, equation (12) is subject to the constraints shown in equation (13):

[0083]

[0084] k x (i,j) represents the phase blur gradient of the pixel in the i-th row and j-th column along the row direction, k y (i,j) represents the phase blur gradient of the pixel in the i-th row and j-th column along the column direction. Introducing the variable x... + (i,j),x - (i,j),y + (i,j),y - (i,j), the nonlinear minimization problem represented by equation (12) is transformed into a linear minimization problem for solution:

[0085]

[0086] like Figure 4 As shown, the minimum cost flow algorithm is a commonly used algorithm for solving linear minimization problems. Figure 4 Hollow circles represent nodes in a network flow; arrows indicate the direction of flow; x + (i,j),x - (i,j),y + (i,j),y - (i,j) represents the flow size flowing into or out of a node; the unit cost of the flow is the weighting coefficient w. x (i,j),w y (i,j). The minimum cost flow problem takes the degree of each node and the unit cost of the flow as input, and outputs the flow rates of all flows that satisfy the minimum cost. Currently, the main algorithms for solving the minimum cost flow problem include the primal-dual algorithm, the relaxation algorithm, and the network simplex algorithm. x is calculated using existing minimum cost flow algorithms. + (i,j),x - (i,j),y + (i,j),y - After (i,j), the gradient of the phase ambiguity number can be further derived.

[0087] Similarly, for the case of P=2, Equation (11) represents the L2 paradigm optimization model, which is a weighted least squares problem that can be solved using the Gauss-Seidel relaxation iterative algorithm. After obtaining the phase blur number k(s) of all pixels, the absolute terrain phase is recovered using Equation (15).

[0088]

[0089] Where A is a constant. Take a ground control point P (usually the ground point corresponding to the center point of the satellite image), and use equation (8) to calculate the true terrain phase φ of the corresponding pixel at the ground point. real , will φ real and multi-baseline unwrapping phase ψ unwrapped The difference is the offset A = φ real -ψ unwrapped .

[0090] Step 4: Calculate the surface elevation based on the topographic phase.

[0091] First, the ground height corresponding to each pixel in the radar coordinate system is obtained using equation (16):

[0092]

[0093] It should be noted that the entire process generates three interferograms, resulting in three sets of elevation results h. AB ,h AC ,h AD Based on the vertical baseline lengths corresponding to the interferogram, a weighted average is used to obtain the final elevation value:

[0094]

[0095] Finally, after geocoding, the DEM product is obtained. Using satellite orbital parameters and imaging geometry (range-Doppler equation, slant range equation, and ellipsoid equation), the three-dimensional coordinates of the ground points corresponding to each pixel in the image are calculated, and after resampling using interpolation methods (such as Kriging interpolation), a DEM product in the form of a regular grid is obtained.

[0096] The implementation of the various embodiments of the present invention is based on programmed processing through a device with processor functionality. Therefore, in practical engineering, the technical solutions and functions of the various embodiments of the present invention are encapsulated into various modules. Based on this reality, and building upon the above embodiments, the embodiments of the present invention provide a macro... Figure 1 The distributed wheel formation InSAR satellite elevation inversion device is used to perform the macroscopic method described in the above embodiments. Figure 1 Distributed wheel formation InSAR satellite elevation inversion method.

[0097] The device includes: a first main module for correcting the single-base multi-pass InSAR geometric model by establishing an equivalent phase center; a second main module for registering master and slave images based on the corrected single-base multi-pass InSAR geometric model, generating three master and slave image pairs interferograms with different baseline lengths, and performing terrain removal phase processing; a third main module for unwrapping the multi-baseline phase of the terrain removal interferogram to restore the absolute terrain phase; and a fourth main module for converting the unwrapped phase into the target scene elevation, and generating a DEM of the target area after geocoding.

[0098] Macros provided in the embodiments of the present invention Figure 1 The distributed wheel formation InSAR satellite elevation inversion equipment addresses the limitation of existing DEM accuracy by employing several modules mentioned above. Through the coordinated use of single-baseline multi-pass InSAR geometric model correction, SAR image registration and deflating interferogram generation, multi-baseline phase unwrapping, phase-to-elevation conversion, and geocoding, it effectively improves the elevation inversion accuracy of DEMs in complex terrains.

[0099] It should be noted that the device embodiments provided by the present invention are used not only to implement the methods in the above method embodiments, but also to implement the methods in other method embodiments provided by the present invention. The only difference is that corresponding functional modules are set. The principle is basically the same as that of the above device embodiments provided by the present invention. As long as those skilled in the art can improve the device in the above device embodiments by referring to the specific technical solutions in other method embodiments, combining technical features to obtain corresponding technical means and technical solutions composed of these technical means, and ensuring the practicality of the technical solutions, they can obtain corresponding device-type embodiments for implementing the methods in other method-type embodiments.

[0100] Based on the same inventive concept as the foregoing embodiments, this embodiment of the invention also provides a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion device includes a memory and a processor. The memory stores program instructions that are executed by the processor, which calls the program instructions to execute the macro. Figure 1 The steps of the InSAR satellite elevation inversion method for distributed wheel formation.

[0101] In embodiments of the present invention, the memory can be non-volatile memory, such as a hard disk drive (HDD) or a solid-state drive (SSD), or it can be volatile memory, such as random-access memory (RAM). Memory is any other medium capable of carrying or storing desired program code having an instruction or data structure form and accessible by a computer, but is not limited thereto. The memory in embodiments of the present invention can also be a circuit or any other device capable of implementing a storage function for storing program instructions and / or data.

[0102] In this embodiment of the invention, the processor may be a general-purpose processor, a digital signal processor, an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components, capable of implementing or executing the methods, steps, and logic block diagrams disclosed in this embodiment of the invention. The general-purpose processor may be a microprocessor or any conventional processor. The steps of the methods disclosed in this embodiment of the invention can be directly manifested as being executed by a hardware processor, or executed by a combination of hardware and software modules within the processor.

[0103] Based on the same inventive concept as the foregoing embodiments, this embodiment of the invention also provides a non-transitory computer-readable storage medium that stores computer instructions that cause the computer to execute the macro. Figure 1 The steps of the InSAR satellite elevation inversion method for distributed wheel formation.

[0104] In summary, the present invention discloses a macro Figure 1 The method for elevation inversion of distributed wheel formation InSAR satellites first establishes an equivalent phase center, simplifying the original dual-baseline imaging geometric model of the auxiliary satellite (with separate transmitter and receiver) into a self-transmitting and self-receiving single-baseline multi-pass InSAR geometric model. Then, based on the single-baseline multi-pass InSAR geometric model, primary and secondary images are registered, generating three interferograms of primary and secondary image pairs with different baseline lengths. These are then processed to remove the terrain phase by incorporating orbital information and other parameters. The interferograms after terrain phase removal are then unwrapped using multi-baseline phase to recover the absolute terrain phase. Finally, the unwrapped phase is converted into the target scene elevation, and after geocoding, a DEM (Digital Elevation Model) of the target area is generated. This method effectively improves the elevation inversion accuracy of digital elevation models (DEMs) for complex terrain.

[0105] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for elevation inversion of InSAR satellites using a distributed wheel formation of the Hongtu-1 satellite, wherein the wheel formation comprises one primary satellite and three secondary satellites, characterized in that... include: For the three auxiliary satellites, equivalent phase centers are established for each satellite. The original bistatic imaging geometric model with separate transmission and reception for each auxiliary satellite is corrected to a single-base multi-pass InSAR geometric model with self-transmission and self-reception. This includes: fitting the flight trajectory of the main satellite within the imaging time period using a polynomial with time as the independent variable based on the orbital node data of the main satellite; substituting the time corresponding to the orbital nodes of the auxiliary satellites into the polynomial to obtain the position and velocity of the main satellite at the same time, and then obtaining the position of the equivalent phase center corresponding to the orbital nodes of the auxiliary satellites; using the polynomial for fitting, and replacing the orbital parameters of the auxiliary satellites with the obtained polynomial parameters to complete the correction of the auxiliary satellite geometric model. Based on the corrected single-base multi-pass InSAR geometric model, the master and slave images are registered, and three master and slave image pairs interferograms with different baseline lengths corresponding to the three slave satellites are generated. The flat phase of each interferogram is then processed. Multi-baseline phase unwrapping is performed on the interferogram after removing the flat terrain phase to recover the absolute terrain phase. This includes: constructing the following objective function and solving for the optimal phase ambiguity number gradient: ， in, Represents a pair of adjacent pixels, which exist in two directions: azimuth and slant range. Indicates the entanglement phase gradient; It is the phase ambiguity number gradient; after determining the optimal phase ambiguity number gradient, the LP minimum norm optimization model is used as the objective function to perform single-baseline phase unwrapping on each interferogram. After obtaining the phase ambiguity number of all pixels, the absolute terrain phase is restored. The unwrapped phase is converted into the elevation of the target scene, and after geocoding, a DEM of the target area is generated.

2. The Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method according to claim 1, characterized in that, Registration of primary and secondary images based on the corrected secondary satellite imaging geometric model includes: Coarse registration: By determining the positions of corresponding points in the primary and secondary satellite images, the row and column coordinate offsets of the secondary satellite image relative to the primary satellite image are obtained; Fine registration: Control points are evenly distributed on the main satellite image. The search range of corresponding points in the auxiliary satellite image is determined based on the offset obtained from coarse registration. The precise position of corresponding points in the auxiliary satellite image is determined according to the principle of maximizing coherence coefficient. The coordinate mapping relationship between the main and auxiliary satellites is established. Finally, coordinate transformation and interpolation resampling are performed on the auxiliary satellite image.

3. The Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method according to claim 2, characterized in that, After completing the registration of the master and slave images, the following is also included: Complex conjugate multiplication is performed to obtain a differential interferogram, and the interferogram is then processed to remove the flat phase.

4. The Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method according to claim 1, characterized in that, When the ground point P=1, the problem of solving the LP minimum norm optimization model is transformed into a linear minimization problem; when the ground point P=2, the problem of solving the LP minimum norm optimization model is transformed into a weighted least squares problem.

5. The Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method according to claim 1, characterized in that, The unwrapped phase is converted into the target scene elevation, and after geocoding, a DEM of the target area is generated, including: Obtain the ground elevation corresponding to each pixel in the radar coordinate system; Based on the vertical baseline length corresponding to the interferogram and the ground height, the final elevation value is obtained; The final elevation value is geocoded to obtain the DEM product.

6. The Hongtu-1 distributed wheel formation InSAR satellite elevation inversion device, used to implement the method described in any one of claims 1 to 5, characterized in that, include: The first main module is used to correct the single-base multi-pass InSAR geometric model by establishing an equivalent phase center; The second main module is used to register master and slave images based on the corrected single-base multi-pass InSAR geometric model, generate three master and slave image pairs interferograms with different baseline lengths, and perform flat-ground phase processing. The third main module is used to perform multi-baseline phase unwrapping on the interferogram after removing the flat phase, and to restore the absolute terrain phase. The fourth main module is used to convert the unwrapped phase into the elevation of the target scene, and after geocoding, generate the DEM of the target area.

7. A Hongtu-1 distributed wheel formation InSAR satellite elevation inversion device, characterized in that, The system includes a memory and a processor, the memory storing program instructions that are executed by the processor, the processor invoking the program instructions to perform the steps of the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method as described in any one of claims 1 to 5.

8. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions that cause the computer to perform the steps of the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method as described in any one of claims 1 to 5.

Images