Mono No.1 distributed wheel formation InSAR (Interferometric Synthetic Aperture Radar) satellite elevation inversion method and equipment

Abstract

The invention belongs to the field of synthetic aperture radar interferometry, and discloses a macro map No.1 distributed wheel formation InSAR satellite elevation inversion method, which can effectively improve the elevation inversion precision of a digital elevation model (DEM) of a complex terrain. The method comprises the following steps: firstly, establishing an equivalent phase center, and simplifying an original receiving and transmitting bistatic imaging geometric model of a cosatellite into a self-transmitting and self-receiving single-base multi-pass InSAR geometric model; then, registering the main and auxiliary images based on a single-base multi-pass InSAR geometric model, generating three main and auxiliary image pair interferograms with different baseline lengths, and performing flat-ground phase removal processing by combining parameters such as orbit information and the like; performing multi-baseline phase unwrapping on the interference pattern after the flat ground phase is removed, and recovering an absolute terrain phase; and finally, converting the unwrapped phase into a target scene elevation, and generating a target area DEM after geocoding processing.

Description

Technical Field

[0001] The present invention belongs to the field of synthetic aperture radar interferometry, and in particular relates to a macro Figure 1 Distributed wheel formation InSAR satellite elevation inversion method and equipment. Background Art

[0002] As an important part of basic geographic information data, Digital Elevation Model (DEM) is widely used in terrain analysis, visibility calculation and other multi-dimensional spatial analysis. At present, optical image stereo matching, Laser Radar (LiDAR) and spaceborne synthetic aperture radar interferometry (InSAR) are the three main technologies for generating large-scale DEM. Among them, the spaceborne InSAR method can quickly and accurately obtain high-precision DEM all day and all weather, and is particularly suitable for stable acquisition of DEM under complex meteorological and weather conditions such as cloudy and foggy weather conditions, and has important application value for geological disaster monitoring and early warning.

[0003] The spaceborne InSAR method uses interference analysis on SAR images with a certain parallax angle acquired in the same area to extract and reconstruct the elevation information of the surface based on the interference phase information. Traditional single-satellite InSAR systems are mainly divided into single-antenna repeating orbit interferometer systems and dual-antenna single-orbit interferometer systems. The accuracy of the DEMs obtained is restricted by many factors.

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

[0005] In order to solve the problem of limited accuracy of existing DEM, the present 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, there is provided a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion method includes:

[0007] The geometric model of single-base multi-flyover InSAR is corrected by establishing equivalent phase center;

[0008] Based on the corrected single-base multi-flyover InSAR geometric model, the primary and secondary images are registered to generate three primary and secondary image pairs with different baseline lengths, and then the ground-flattening phase processing is performed.

[0009] Perform multi-baseline phase unwrapping on the interferogram after removing the flat ground phase to restore the absolute terrain phase;

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

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

[0012] According to the orbit node data of the main satellite, the flight orbit of the main satellite during the imaging period is fitted using a polynomial with time as the independent variable;

[0013] Substitute the time corresponding to the auxiliary satellite orbit node into the polynomial to 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 acquired polynomial parameters are used to replace the auxiliary satellite orbit parameters to complete the auxiliary satellite geometric model correction.

[0015] As a further technical solution, the primary and secondary images are registered based on the corrected secondary satellite imaging geometry model, including:

[0016] Rough registration: By determining the positions of the same-name points of 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 primary satellite image. The search range of the corresponding same-name points in the auxiliary satellite image is determined according to the offset obtained by coarse registration. The precise position of the same-name points in the auxiliary satellite image is determined according to the principle of maximum coherence coefficient. The coordinate mapping relationship between the primary and auxiliary satellites is established. Finally, the auxiliary satellite image is coordinate transformed and resampled by interpolation.

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

[0019] The differential interference pattern is obtained by complex conjugate multiplication, and the interference pattern is subjected to de-flattening phase processing.

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

[0021] The following objective function is constructed and solved to obtain the optimal phase ambiguity number gradient:

[0022]

[0023] Among them, (s, s-1) represents a pair of adjacent pixels, which exists in two directions: azimuth and slant range; represents the winding phase gradient; Δk i (s, s-1) is the phase ambiguity number gradient;

[0024] 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 numbers of all pixels, the absolute terrain phase is restored.

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

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

[0027] Get the ground height corresponding to each pixel in the radar coordinate system;

[0028] According to the vertical baseline length corresponding to the interference pattern and the ground height, a final elevation value is obtained;

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

[0030] According to one aspect of the present invention, there is provided a macro Figure 1 Distributed wheel formation InSAR satellite elevation inversion equipment, including:

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

[0032] The second main module is used to register the primary and secondary images based on the corrected single-base multi-pass InSAR geometric model, generate three primary and secondary image pair interferograms with different baseline lengths, and perform ground-free phase processing;

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

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

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

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

[0037] Compared with the prior art, the present invention has the following beneficial effects:

[0038] The present invention discloses a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion method can effectively improve the elevation inversion accuracy of the digital elevation model (DEM) in complex terrain. The method first establishes the equivalent phase center and simplifies the dual-base imaging geometric model with the auxiliary satellite's original transmission and reception separated into a single-base multi-flyover InSAR geometric model with self-transmission and self-reception; then, the primary and auxiliary images are registered based on the single-base multi-flyover InSAR geometric model to generate three primary and auxiliary image pairs with different baseline lengths, and the ground-removing phase processing is performed in combination with orbit information and other parameters; the interferogram after ground-removing phase is subjected to multi-baseline phase unwrapping to restore the absolute terrain phase; finally, the unwrapped phase is converted into the target scene elevation, and after geocoding processing, the target area DEM is generated. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, a brief introduction is given below to the drawings used in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0040] Figure 1 The macro provided by the embodiment of the present invention Figure 1 Flowchart of the distributed wheel formation InSAR satellite elevation inversion method.

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

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

[0043] Figure 4A schematic diagram of the minimum cost flow phase unwrapping provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0044] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. In addition, the technical features in the various embodiments or single embodiments provided by the present invention are arbitrarily combined with each other to form a new technical solution. This combination is not restricted by the sequence of steps and / or the structural composition mode, but must be based on the ability of ordinary technicians in this field to achieve. When the combination of technical solutions is contradictory or cannot be achieved, it should be considered that this combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0045] The embodiment of the present invention provides a macro Figure 1 The distributed wheel formation InSAR satellite elevation inversion method is proposed. First, the equivalent phase center is established, and the dual-base model with the auxiliary satellite's original transmission and reception separated is simplified into a single-base multi-flyover InSAR geometric model with self-transmission and self-reception; then, the primary and auxiliary images are registered according to the single-base multi-flyover InSAR geometric model to generate three main and auxiliary image pairs with different baseline lengths. The interferograms are then processed to remove the flattening phase in combination with orbital information and other parameters; multi-baseline phase unwrapping is performed on the interferogram after the flattening phase is removed to restore the absolute terrain phase; finally, the unwrapped phase is converted into the target scene elevation, and the DEM of the target area is generated through geocoding processing.

[0046] The known data required by the embodiment of the present invention include: HT-1SAR image (single-view complex data) and header file, external DEM (such as open source SRTM DEM). Figure 1 As shown, the process of the embodiment of the present invention includes the following steps:

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

[0048] like Figure 2 As shown in the figure, when the SAR satellite is in the one-transmit-two-receive working mode, the position 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 formula (1):

[0049]

[0050] Among them, p t Indicates the ground point position; p m,v m represents the position and velocity vector of the main satellite ground point at imaging time t; p s ,v s represents the auxiliary satellite position and velocity vector; λ is the radar wavelength; f dc Indicates the auxiliary satellite Doppler center frequency; R 1 ,R 2 Respectively represent the slant distances from the primary and auxiliary satellites to the ground point.

[0051] The midpoint of the line connecting the phase centers of the main and auxiliary radar antennas is constructed as the equivalent phase center of the auxiliary radar antenna, specifically as formula (2):

[0052]

[0053] Among them, p se (t) represents the equivalent phase center of the auxiliary satellite; R 2se It represents the slant distance between the auxiliary satellite and the ground point target after the equivalent phase center. The HT-1 image header file provides the satellite precise orbit data at a time interval of 1 second. The specific method of correcting the auxiliary satellite imaging geometric model is to fit the flight orbit of the main satellite during the imaging period using a polynomial with the independent variable being time (UTC time) according to the orbit node data of the main satellite:

[0054]

[0055] Among them, (X m ,Y m ,Z m ),(V mx ,V my ,V mz ) represents the orbital position and flight speed of the main satellite at time t; a, b, c, d, e, and f are polynomial parameters obtained by fitting the orbital node data. s Substitute into formula (3) to obtain the position and velocity of the main satellite at the same time, and use formula (2) to obtain the equivalent phase center position corresponding to the auxiliary satellite orbit node. Further, use the polynomial shown in formula (3) for fitting, replace the auxiliary satellite orbit parameters with the obtained polynomial parameters, and complete the auxiliary satellite geometric model correction.

[0056] After the auxiliary satellite phase center is replaced by the equivalent phase center, the dual-base imaging model described by equation (1) can be approximated as the classic single-base imaging model, specifically equation (4):

[0057]

[0058] Step 2: Primary and secondary SAR image registration and flattened interferogram generation.

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

[0060]

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

[0062] After obtaining the imaging time t of the auxiliary satellite for point P and the corresponding satellite position, the pixel coordinates (r s ,c s ):

[0063]

[0064] Among them, R near is the slant range corresponding to the first column of pixels of the auxiliary satellite; dR is the sampling interval in the slant range direction; t 0 is the starting imaging time of the auxiliary satellite image; dt is the azimuth sampling time interval. near ,dR,t 0 ,dt parameters can be obtained from the image header file. After determining the position of the same-name point in the primary and secondary satellite images, obtain the row and column coordinate offset of the secondary satellite image relative to the primary satellite image.

[0065] In the fine registration step, control points are first evenly distributed on the primary satellite image. The search range of the corresponding same-name points in the auxiliary satellite image is determined according to the offset obtained by the rough registration. The precise position of the same-name points in the auxiliary satellite image is determined according to the principle of maximum coherence coefficient, and the coordinate mapping relationship between the primary and auxiliary satellites is established. Finally, the coordinate transformation and interpolation resampling of the auxiliary satellite image are performed. After completing the registration of the primary and auxiliary images, complex conjugate multiplication is performed to obtain the differential interferogram.

[0066] The principle diagram of the single-base repetitive orbit InSAR interferometry phase is as follows: Figure 3 As shown, O represents the center of the earth, the main satellite radar incident angle θ is used to image the ground point P from the side, and the slant distance is R 1 , M represents the phase center position of the satellite antenna during imaging, H represents the height of the main satellite antenna flight platform, P is the surface target, h is the surface elevation value of point P that needs to be solved, B represents the spatial baseline between the main and auxiliary satellites, and the angle between the baseline B and the horizontal direction is recorded as α. Parallel baseline B parallel to the main satellite radar slant range MP || and the vertical baseline B perpendicular to the radar slant range direction SP ⊥ , calculated by formula (7):

[0067]

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

[0069]

[0070] Among them, P 0 is the equidistant projection point of point P on the reference ellipsoid, that is, r 1 =R 1 , R 2 Indicates the auxiliary satellite slant distance, the main satellite radar to P 0 The side angle is θ 0 By removing the reference ellipsoid phase φ from the interference pattern ref , we can get the terrain phase φ topo , that is, to flatten the interference pattern.

[0071] Step 3: multi-baseline phase unwrapping, recover the wrapped terrain phase to obtain the absolute terrain phase.

[0072] For a set of HT-1 satellite data, HT-1A data is used as the main image, and B / C / D data are used as auxiliary images for registration interference, generating a total of three flattened interferograms. The corresponding vertical baseline is B AB ,B AC ,B AD . Original winding phase of pixel s The relationship between the true absolute phase ψ(s) and the true absolute phase ψ(s) is as follows:

[0073]

[0074] Where k(s) represents the phase ambiguity number of pixel s. For pixel s and the corresponding adjacent pixel s-1 in the interference pattern, the objective function shown in formula (10) is constructed:

[0075]

[0076] Among them, (s, s-1) represents a pair of adjacent pixels, which exists in two directions: azimuth and slant range; represents the winding phase gradient; Δk i (s, s-1) is the phase ambiguity gradient, which can be any integer. i (s,s-1),Δ j (s, s-1) represents the pixel pair variable consisting of pixel s and its pixel s-1 in the row and column directions. In order to determine the optimal phase fuzzy number gradient, a brute force solution is used: the fuzzy number gradient Δk is set according to prior knowledge such as external DEM i The search range of (s,s-1) [min i (s),max i (s)], and exhaustively traverse all fuzzy number gradient combinations (Δk AB (s,s-1),Δk AC (s,s-1),Δk AD (s,s-1)), the combination that minimizes equation (10) is taken as the optimal solution of the phase ambiguity number gradient. When setting the ambiguity search range, it is necessary to consider the range [min i (s),max i (s)] has the only optimal solution of formula (10). According to 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] In order to obtain the optimal phase ambiguity gradient Δk i After (s,s-1)(i∈{AB,AC,AD}), the LP minimum norm optimization model is constructed to perform single baseline phase unwrapping on each interferogram. The objective function is shown in formula (11):

[0078]

[0079] Among them, w i (s,s-1)∈[0,1] is the weight coefficient, which can be determined by the interference quality map or other prior knowledge, k i (s) represents the phase ambiguity gradient of pixel s in interference pattern i, and Z represents the integer domain. For the case of P = 1, define Δ x k(i,j),Δ yk(i,j) is the phase ambiguity gradient 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); is the discrete partial derivative of the winding phase in the row and column directions. At this time, equation (11) is transformed into the nonlinear minimization problem as shown below:

[0080]

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

[0082] Since the real 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 ambiguity gradient of the pixel in the i-th row and j-th column in the row direction, k y (i, j) represents the phase ambiguity gradient of the pixel in the i-th row and j-th column in the column direction. Introducing the variable x + (i,j),x - (i,j),y + (i,j),y - (i, j), the nonlinear minimization problem expressed in equation (12) is converted 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 The hollow circles represent nodes in the network flow; the arrows represent the flow direction; 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 input of the minimum cost flow problem is the degree of each node and the unit cost of the flow, and the output is the flow rate of each flow that satisfies the cost and minimum. At present, the main algorithms for solving the minimum cost flow problem include the primal dual algorithm, the relaxation algorithm, and the network simplex algorithm. + (i,j),x- (i,j),y + (i,j),y - (i, j), the phase ambiguity number gradient can be further derived.

[0087] Similarly, for the case of P = 2, equation (11) is expressed as an L2 norm optimization model, which is a weighted least squares problem and can be solved using the Gauss-Seidel relaxation iterative algorithm. After obtaining the phase ambiguity number k(s) of all pixels, the absolute terrain phase is restored by equation (15).

[0088]

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

[0090] Step 4 calculates the surface elevation based on the terrain phase.

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

[0092]

[0093] It should be noted that the whole process generates three interference graphs, corresponding to three sets of elevation results h AB ,h AC ,h AD According to the vertical baseline length corresponding to the interference pattern, the final elevation value is obtained by weighted average:

[0094]

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

[0096] The implementation basis of each embodiment of the present invention is to implement programmed processing through a device with a processor function. Therefore, in engineering practice, the technical solutions and functions of each embodiment of the present invention are encapsulated into various modules. Based on this reality, on the basis of the above embodiments, the embodiment of the present invention provides a macro Figure 1A distributed wheel formation InSAR satellite elevation inversion device is used to execute the macro in the above method embodiment. Figure 1 A distributed wheel formation InSAR satellite elevation inversion method.

[0097] The device includes: a first main module, which is used to correct the single-base multi-pass InSAR geometric model by establishing an equivalent phase center; a second main module, which is used to align the primary and auxiliary images based on the corrected single-base multi-pass InSAR geometric model, generate three primary and auxiliary image pair interference maps with different baseline lengths, and perform ground-removing phase processing; a third main module, which is used to perform multi-baseline phase unwrapping on the interference map after ground-removing phase, and restore the absolute terrain phase; a fourth main module, which is used to convert the unwrapped phase into the target scene elevation, and generate the target area DEM after geocoding processing.

[0098] The macro provided by the embodiment of the present invention Figure 1 The No. 1 distributed wheel formation InSAR satellite elevation inversion equipment aims to solve the problem of limited accuracy of existing DEM. It adopts several of the aforementioned modules and effectively improves the DEM elevation inversion accuracy of complex terrain through the coordination of single-base multi-pass InSAR geometric model correction, SAR image registration and de-flattened interferogram generation, multi-baseline phase unwrapping, phase conversion to elevation and geocoding.

[0099] It should be noted that the device embodiments provided by the present invention are not only used to implement the methods in the above-mentioned method embodiments, but also used to implement the methods in other method embodiments provided by the present invention. The only difference lies in the setting of corresponding functional modules, and the principles thereof are basically the same as the principles of the above-mentioned device embodiments provided by the present invention. As long as technical personnel in this field refer to the specific technical solutions in other method embodiments on the basis of the above-mentioned device embodiments, obtain corresponding technical means and technical solutions composed of these technical means by combining technical features, and on the premise of ensuring the practicality of the technical solutions, they will improve the devices in the above-mentioned device embodiments to obtain corresponding device class embodiments for implementing the methods in other method class embodiments.

[0100] Based on the same inventive concept as the above-mentioned embodiment, the embodiment of the present invention further provides a macro Figure 1 A distributed wheel formation InSAR satellite elevation inversion device, comprising a memory and a processor, wherein the memory stores program instructions executed by the processor, and the processor calls the program instructions to execute the macro Figure 1 The steps of the distributed wheel formation InSAR satellite elevation inversion method.

[0101] In an embodiment of the present invention, the memory may be a non-volatile memory, such as a hard disk drive (HDD) or a solid-state drive (SSD), etc., or a volatile memory (volatile memory), such as a random-access memory (RAM). The memory is any other medium that can be used to carry or store the desired program code in the form of instructions or data structures and can be accessed by a computer, but is not limited thereto. The memory in an embodiment of the present invention may also be a circuit or any other device that can implement a storage function, for storing program instructions and / or data.

[0102] In the embodiments of the present invention, the processor may be a general-purpose processor, a digital signal processor, an application-specific integrated circuit, a field programmable gate array or other programmable logic device, a discrete gate or transistor logic device, or a discrete hardware component, and may implement or execute the methods, steps, and logic block diagrams disclosed in the embodiments of the present invention. A general-purpose processor may be a microprocessor or any conventional processor, etc. The steps of the method disclosed in the embodiments of the present invention may be directly embodied as being executed by a hardware processor, or may be executed by a combination of hardware and software modules in the processor.

[0103] Based on the same inventive concept as the above embodiment, the embodiment of the present invention further provides a non-transitory computer-readable storage medium, wherein the non-transitory computer-readable storage medium stores computer instructions, wherein the computer instructions enable the computer to execute the macro Figure 1 The steps of the distributed wheel formation InSAR satellite elevation inversion method.

[0104] In summary, the present invention discloses a macro Figure 1 The invention discloses a method for inverting the elevation of a distributed wheel formation InSAR satellite. First, an equivalent phase center is established to simplify the dual-base imaging geometric model of the auxiliary satellite with separate transmission and reception into a single-base multi-pass InSAR geometric model of self-transmission and self-reception; then, the primary and secondary images are registered based on the single-base multi-pass InSAR geometric model to generate interference diagrams of three primary and secondary image pairs with different baseline lengths, and the ground-removing phase is processed in combination with parameters such as orbit information; the interference diagram after ground-removing phase is subjected to multi-baseline phase unwrapping to restore the absolute terrain phase; finally, the unwrapped phase is converted into the elevation of the target scene, and after geocoding processing, the DEM of the target area is generated. The method of the invention can effectively improve the elevation inversion accuracy of the digital elevation model (DEM) of complex terrain.

[0105] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the technical solutions of the embodiments of the present invention.

Claims

1. Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method, characterized by: include: The geometric model of single-base multi-flyover InSAR is corrected by establishing equivalent phase center; Based on the corrected single-base multi-flyover InSAR geometric model, the primary and secondary images are registered to generate three primary and secondary image pairs with different baseline lengths, and then the ground-flattening phase processing is performed. Perform multi-baseline phase unwrapping on the interferogram after removing the flat ground phase to restore the absolute terrain phase; The unwrapped phase is converted into the target scene elevation, and after geocoding processing, the target area DEM is generated.

2. According to claim 1, the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method is characterized in that: By establishing the equivalent phase center, the auxiliary satellite imaging geometric model is corrected, including: According to the orbit node data of the main satellite, the flight orbit of the main satellite during the imaging period is fitted using a polynomial with time as the independent variable; Substitute the time corresponding to the auxiliary satellite orbit node into the polynomial to 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; The polynomial is used for fitting, and the acquired polynomial parameters are used to replace the auxiliary satellite orbit parameters to complete the auxiliary satellite geometric model correction.

3. According to claim 1, the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method is characterized in that: The primary and secondary images are registered based on the corrected secondary satellite imaging geometry model, including: Rough registration: By determining the positions of the same-name points of 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 primary satellite image. The search range of the corresponding same-name points in the auxiliary satellite image is determined according to the offset obtained by coarse registration. The precise position of the same-name points in the auxiliary satellite image is determined according to the principle of maximum coherence coefficient. The coordinate mapping relationship between the primary and auxiliary satellites is established. Finally, the auxiliary satellite image is coordinate transformed and resampled by interpolation.

4. According to claim 3, the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method is characterized in that: After completing the primary and secondary image registration, it also includes: The differential interference pattern is obtained by complex conjugate multiplication, and the interference pattern is subjected to de-flattening phase processing.

5. According to claim 1, the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method is characterized in that: Perform multi-baseline phase unwrapping on the interferogram after removing the flat ground phase to restore the absolute terrain phase, including: The following objective function is constructed and solved to obtain the optimal phase ambiguity number gradient: Among them, (s, s-1) represents a pair of adjacent pixels, which exists in two directions: azimuth and slant range; represents the winding phase gradient; Δk i (s, s-1) 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 numbers of all pixels, the absolute terrain phase is restored.

6. According to claim 5, the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method is characterized in that: When the ground point P=1, the solution problem of the LP minimum norm optimization model is converted into a linear minimization problem; when the ground point P=2, the solution problem of the LP minimum norm optimization model is converted into a weighted least squares problem.

7. The method for inverting the altitude of the Hongtu-1 distributed wheel formation InSAR satellite according to claim 1 is characterized in that: The unwrapped phase is converted into the target scene elevation, and after geocoding processing, the target area DEM is generated, including: Get the ground height corresponding to each pixel in the radar coordinate system; According to the vertical baseline length corresponding to the interference pattern and the ground height, a final elevation value is obtained; The final elevation value is geocoded to obtain a DEM product.

8. Hongtu-1 distributed wheel formation InSAR satellite elevation inversion equipment, characterized by: include: The first main module is used to correct the single-base multi-flyover InSAR geometric model by establishing an equivalent phase center; The second main module is used to register the primary and secondary images based on the corrected single-base multi-pass InSAR geometric model, generate three primary and secondary image pair interferograms with different baseline lengths, and perform ground-free phase processing; The third main module is used to perform multi-baseline phase unwrapping on the interferogram after removing the flat ground phase and recover the absolute terrain phase; The fourth main module is used to convert the unwrapped phase into the target scene elevation, and generate the target area DEM after geocoding processing.

9. A Hongtu-1 distributed wheel formation InSAR satellite elevation inversion device, characterized in that: It includes a memory and a processor, the memory stores program instructions executed by the processor, and the processor calls the program instructions to execute the steps of the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method as described in any one of claims 1 to 7.

10. A non-transitory computer-readable storage medium, characterized in that: The non-transitory computer-readable storage medium stores computer instructions, which enable the computer to execute the steps of the Hongtu-1 distributed wheel formation InSAR satellite elevation inversion method as described in any one of claims 1 to 7.

Images