Method, device, equipment and medium for correcting plane position of insar dem based on lidar data

Abstract

The application discloses a kind of based on laser radar data correction InSAR DEM plane position method, device, equipment and medium, method is: obtaining InSAR data and interference processing, generates DEM;Obtain spaceborne laser radar data and combine DEM corresponding ground point, with spaceborne laser radar same track strip as unit, based on the relative error of the length of two terrain profiles projection Registration, obtain the offset of DEM relative to spaceborne laser radar data;Based on the offset of registration plane, utilize least square adjustment method to construct the offset polynomial model of DEM in longitude, latitude direction;Offset polynomial model is carried out to DEM by geometric transformation, finally realizes plane position correction.The application only uses spaceborne laser radar data to carry out plane position correction to DEM, overcome the problem that existing method relies on existing terrain product.

Description

Technical Field

[0001] This invention belongs to the field of geodesy in synthetic aperture radar interferometry (InSAR), specifically relating to a method, apparatus, equipment, and medium for correcting the planar position of an InSAR DEM based on lidar data. Background Technology

[0002] Digital Elevation Models (DEMs) are one of the fundamental data sources for surveying and mapping and geographic information construction, with important applications in basic engineering construction, disaster monitoring, agriculture, and forestry. Interferometric Synthetic Aperture Radar (InSAR) technology, due to its all-weather, wide-area, and high-resolution imaging capabilities, has become the primary technique for DEM production.

[0003] Currently, the commonly used method for generating DEMs using InSAR technology is DEM estimation based on differential interferometry (D-InSAR). Its core involves simulating the phase using an existing external DEM, differentially differentiating it from the InSAR interferometric phase, then unwrapping and performing phase height transformation on the differential phase before superimposing it onto the external DEM to generate a new DEM in the SAR coordinate system. Finally, geocoding from the external DEM to the geographic coordinate system generates the final DEM. However, the geocoding process, which transforms the DEM from the SAR coordinate system to the geographic coordinate system, relies on the positioning and registration of the SAR image and the external DEM. This process is typically affected by SAR image trajectory errors and DEM accuracy, and it easily introduces the trend-based planar errors of the external DEM into the InSAR-estimated DEM.

[0004] On the other hand, spaceborne lidar technology measures the height of the Earth's surface by emitting laser pulses into the Earth's surface and receiving echo signals. It can collect altimetry points densely around the globe and has high horizontal and vertical accuracy, and can be widely used in InSAR DEM control data.

[0005] In summary, existing InSAR DEM planar positioning methods heavily rely on existing terrain products. Therefore, it is necessary to design a high-precision DEM planar position correction method that does not depend on existing external DEM products and is based solely on spaceborne lidar data. Summary of the Invention

[0006] This invention provides a method, apparatus, device, and medium for correcting the planar position of an InSAR DEM based on lidar data, which can achieve high-precision DEM planar position correction without relying on existing external DEM products.

[0007] To achieve the above technical objectives, the present invention adopts the following technical solution:

[0008] A method for estimating DEM based on InSAR and lidar data includes the following steps:

[0009] Step 1: Acquire InSAR image data of the area to be studied and perform interferometric processing to generate an InSAR DEM;

[0010] Step 2: Acquire multiple co-track ground altimetry data from the spaceborne lidar in the area under study. Each co-track ground altimetry data is denoted as one lidar ground altimetry strip. Calculate the elevation projection length of each spaceborne lidar ground altimetry strip along the track direction on the O-xz plane of the spatial rectangular coordinate system, and denote it as the spaceborne lidar elevation projection length. ;

[0011] Step 3: Construct a registration model between the ground altimeter strips of the spaceborne lidar and the elevation strips of the DEM; wherein, the registration parameters of the registration model include the row and column offsets of the ground altimeter strips of the spaceborne lidar relative to the DEM in the geographic coordinate system; set the initial registration parameters;

[0012] Step 4: Based on the registration model with the current registration parameters and the geographical location of the ground altimetry strips of the spaceborne lidar, obtain the elevation strips corresponding to the DEM; then calculate the elevation projection length of the DEM elevation strips along the trajectory of the ground altimetry strips of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system, denoted as the DEM elevation projection length calculated based on the current registration parameters. ;

[0013] Step 5: Calculate the elevation projection length of the spaceborne lidar obtained in Step 2. and the DEM elevation projection length obtained in step 4 The relative error between them is used to construct an error function, and the registration parameters in the registration model are iteratively solved by the nonlinear least squares method.

[0014] Step 6: Based on the registration parameters obtained in Step 5, take the row and column numbers of all points in the ground altimetry strip of the spaceborne lidar in the geographic coordinate system and the offset of the DEM relative to them in the row and column directions as observations, and use the least squares adjustment method to solve the offset polynomials of the DEM relative to the ground altimetry strip of the spaceborne lidar in the row and column directions respectively, and obtain the offset lookup table of the DEM relative to the spaceborne lidar data.

[0015] Step 7: Using the offset lookup table obtained in Step 6, transform the coordinates of the initial DEM obtained in Step 1 to the new coordinate position to obtain the DEM after planar position correction.

[0016] Furthermore, the elevation projection length of each spaceborne lidar ground altimetry strip along its orbital direction on the O-xz plane of the spatial rectangular coordinate system is calculated. The calculation formula is:

[0017]

[0018] In the formula, This represents the number of altimeter points included in the ground altimeter strip of the spaceborne lidar, where x is the row number in the geographic coordinate system. It is the elevation of the measuring point.

[0019] Furthermore, in step 3, based on the planar positional deviation between the ground altimetry strips of the spaceborne lidar and the DEM elevation strips, the following registration model is constructed:

[0020]

[0021] In the formula, , It is the initial row and column number of the ground altimeter strip of the spaceborne lidar in the geographic coordinate system. It is the scaling factor. It is the rotation angle. , It represents the row and column offsets of the ground altimetry strips from the spaceborne lidar relative to the DEM elevation strips in the geographic coordinate system. , , , These are all registration parameters to be solved; , It is the row and column number after the registration parameter transformation of the ground altimeter strip of the spaceborne lidar.

[0022] Further, step 4 calculates the elevation projection length of the DEM elevation strip along the trajectory of the ground altimeter strip of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system. Specifically, it includes:

[0023] First, the terrain contour lines between adjacent ground points of the DEM are fitted with elevation change curves in the geographic coordinate system, as follows:

[0024]

[0025] In the formula, These are the row numbers of the DEM ground points in the geographic coordinate system. , , These are the fitting coefficients to be determined.

[0026] Then, the derivative of the elevation change curve at each surface point is expressed as the slope change rate of each surface point along the strip direction in the DEM. Furthermore, using the coordinates, elevation, and slope change rate of the terrain contour line at both ends of the ground as observed values, an expression for solving the fitting coefficients using the least squares adjustment method is constructed:

[0027]

[0028] In the formula, It is residual error. It is a matrix of row numbers for the ground points at both ends. It is the fitting coefficient matrix. These are the observation error vectors, represented as follows:

[0029]

[0030]

[0031]

[0032] In the formula, , These are the row numbers of the terrain outline at the two ends of the ground. , It refers to the elevation of the terrain outline at the two ends of the ground. , It is the rate of change of slope of the two ground points along the DEM elevation strip direction;

[0033] Then, the least squares estimation was used to fit the coefficient matrix. Perform the calculation:

[0034]

[0035] In the formula, This represents the weight matrix, which is an identity matrix with diagonal elements all being 1, indicating that the terrain information from the two ends of the ground has equal weight in fitting the function.

[0036] Then based on the fitting coefficients , , To solve for the given terrain contour, calculate the projection lengths of the terrain contour onto the O-xz plane between adjacent ground points by integration:

[0037]

[0038] In the formula, The path representing the terrain outline, For the first i The projection length of the terrain contour line onto the O-xz plane;

[0039] Finally, the total projected length of the terrain contour lines between all adjacent ground points in the DEM elevation strip is obtained by summing the projected lengths of the terrain contour lines in the DEM elevation strip.

[0040]

[0041] Furthermore, step 5 specifically includes:

[0042] First, construct the error function:

[0043]

[0044] In the formula, Indicates taking The maximum value in;

[0045] The radar elevation projection length calculated in step 2 As an observation, the DEM elevation projection length calculated by the registered model will be... Represented as a nonlinear function Then the error function is transformed into:

[0046]

[0047] Among them, the independent variable of the nonlinear function It is a ground point in the ground altimetry strip of the spaceborne lidar. i Initial row and column numbers in the geographic coordinate system; The registration parameters to be solved include scaling factors. Rotation angle and row and column offsets , ; Initial row and column numbers The corresponding elevation projection length error;

[0048] The nonlinear least squares method is used to solve for the registration parameters in the error function by minimizing the sum of squared errors, thus obtaining the optimal solution for the registration parameters. :

[0049]

[0050] In the formula, It represents the number of altimeter points included in the ground altimeter strip of the spaceborne lidar.

[0051] Step 6 includes the following specific steps:

[0052] Step 6.1: For each ground altimeter strip of all spaceborne lidars in the study area, run steps 2 to 5 to obtain the registration parameter vectors for different strips.

[0053] Step 6.2: Construct a polynomial fitting offset to correct the planar position of the entire DEM. The expression is as follows:

[0054]

[0055] In the formula, and These are the polynomial fitting coefficients. ; and The row and column numbers of the DEM in the geographic coordinate system. and This represents the offset of the DEM relative to the altimeter points of the spaceborne lidar in each ground altimeter strip;

[0056] Step 6.3: Determine the row and column numbers of the ground altimeter strips from the spaceborne lidar in the geographic coordinate system. and the row and column offsets in the registration parameter vector obtained in step 6.1. , Solve the polynomial fitting coefficients in equation (11) Given the fitting coefficient The expression (11) is the offset lookup table of DEM relative to the spaceborne lidar data.

[0057] An inversion device for correcting the planar position of an InSAR DEM based on lidar data, comprising:

[0058] The preprocessing module is used to: perform interferometric processing on the acquired InSAR data of the area under study to obtain the initial DEM of the image; and also to: calculate the elevation projection length along the orbital direction on the O-xz plane of the spatial rectangular coordinate system for each ground altimeter strip of the acquired spaceborne lidar in the area under study, denoted as the spaceborne lidar elevation projection length. ;

[0059] The registration module is used to: construct a registration model to register the DEM and spaceborne lidar data of the area under study according to the principle of terrain contour consistency, and obtain the planar offset of the DEM relative to the spaceborne lidar data in the geographic coordinate system.

[0060] The registration parameter calculation method for the registration model is as follows: Based on the geographical location of the ground altimeter strip of the spaceborne lidar, obtain the elevation strip corresponding to the DEM, and then calculate the elevation projection length of the DEM elevation strip along the trajectory of the ground altimeter strip of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system. And then according to and The relative error between them is calculated and the registration parameters in the registration model are iteratively solved using a nonlinear least squares method;

[0061] The lookup table construction module is used to: take the row and column numbers of all points in the ground altimetry strip of the spaceborne lidar in the geographic coordinate system and the offset of the DEM relative to them in the row and column directions as observations, and use the least squares adjustment method to solve the offset polynomials of the DEM relative to the ground altimetry strip of the spaceborne lidar in the row and column directions respectively, and obtain the offset lookup table of the DEM relative to the spaceborne lidar data.

[0062] The DEM correction module is used to: use the offset lookup table obtained by the lookup table construction module to transform the initial DEM coordinates obtained by the preprocessing module to a new coordinate position, and obtain the DEM after planar position correction.

[0063] An electronic device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, the processor can implement the method for correcting the InSAR DEM planar position based on lidar data as described above.

[0064] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for correcting the planar position of an InSAR DEM based on lidar data as described in any of the preceding claims.

[0065] Beneficial effects

[0066] This invention is based on InSAR data correction from spaceborne lidar. This method for DEM planar positioning correction achieves DEM planar positioning correction solely based on spaceborne lidar data. First, InSAR data is processed using common data processing methods to obtain an initial DEM. Then, the projected length of the elevation topographic contour line is calculated using spaceborne lidar topographic altimetry data in strip units. Next, initial registration parameters are set, and the corresponding DEM strips are obtained by moving the spaceborne lidar topographic altimetry data according to the registration formula. The topographic contour of adjacent point sets of the DEM strips is fitted, and the coefficients of the fitting function are solved according to the elevation and slope change rates. The projected length of the DEM strip topographic contour is then calculated by integration. Next, the relative error between the spaceborne lidar topographic altimetry data and the corresponding DEM strip topographic contour projection length is calculated and used as an error function to assist in iteratively solving the registration parameter vector. Then, based on the registration parameters of several spaceborne lidar topographic altimetry strips in the study area, an offset polynomial of the DEM relative to the spaceborne lidar data is fitted, and an offset lookup table is constructed. Finally, the DEM is geometrically transformed using the offset lookup table, transforming its coordinates to a new geographic location in the geographic coordinate system, estimating the final DEM result, and completing the DEM planar positioning correction.

[0067] The beneficial effects of this method are as follows: It constructs a method system for correcting the planar position of a DEM based solely on spaceborne lidar data, overcoming the problem that existing InSAR DEM geocoding heavily relies on existing terrain products; it proposes a method for registering spaceborne lidar data with the DEM based on the relative error of the terrain contour, achieving the correction of DEM planar errors. In terms of accuracy, this invention achieves a level of accuracy and operational efficiency comparable to existing methods without relying on existing terrain products, making it an effective method for DEM planar position correction that is autonomous, controllable, efficient, and highly accurate. Attached Figure Description

[0068] Figure 1 This is a flowchart of the method described in an embodiment of the present invention.

[0069] Figure 2 The location of the verification area and the coverage area of ​​the DEM selected for this invention.

[0070] Figure 3 This is a lookup table for the offset of the DEM relative to the spaceborne lidar in the row and column directions of the geographic coordinate system. Figure 3 (a) is the offset lookup table for DEM1 in the column direction; Figure 3 (b) is the offset lookup table for DEM1 in the row direction; Figure 3 (c) is the offset lookup table for DEM2 in the column direction; Figure 3 (d) is the offset lookup table for DEM2 in the row direction.

[0071] Figure 4 The image shows the estimated DEM after planar position correction and the difference between the DEM before and after planar position correction and the airborne DTM. Figure 4 (a) is the estimated DEM1 after planar position correction; Figure 4 (b) is a diagram showing the difference between DEM1 before planar position correction and airborne DTM; Figure 4 (c) is a diagram showing the difference between DEM1 after planar position correction and airborne DTM; Figure 4 (d) is the estimated DEM2 after planar position correction; Figure 4 (e) shows the difference between DEM2 before planar position correction and airborne DTM. Figure 4 (f) is a diagram showing the difference between the DEM2 after planar position correction and the airborne DTM. Detailed Implementation

[0072] The embodiments of the present invention will be described in detail below. These embodiments are based on the technical solutions of the present invention and provide detailed implementation methods and specific operation processes to further explain the technical solutions of the present invention.

[0073] To better illustrate the methods and steps of the present invention, the present invention is further described in detail in a test area located in central Spain using LT-1 InSAR data and ICESat-2 ATL03 data. Note that the specific embodiments described herein are only for explaining the present invention and are not intended to limit the present invention.

[0074] This embodiment provides a method for correcting the planar position of an InSAR DEM based on spaceborne lidar data, referring to... Figure 1 As shown, it includes the following steps:

[0075] Step 1: Based on the coordinate range of the selected experimental area, acquire LT-1 InSAR data, perform interferometric processing on the acquired data, and generate an InSAR DEM.

[0076] This implementation example selected two LT-1 data points located in the experimental region of central Spain. The location of the experimental region is as follows: Figure 2 As shown, the data is subjected to interferometry to generate an InSAR DEM of the area under study. Interferometry is an existing technique, and can be found in the following reference: PA Rosen et al., “Synthetic aperture radar interferometry,” in Proceedings of the IEEE, vol. 88, no. 3, pp. 333-382, March 2000, doi:10.1109 / 5.838084. This embodiment will not elaborate further.

[0077] Step 2: Based on the coordinate range of the selected experimental area, acquire the spaceborne lidar data. Using the same acquisition trajectory (i.e., the same orbital strip) passed by the satellite in one flight as the unit, project it from the spatial rectangular coordinate system to... Plane, and calculate the projected length of the terrain contour line. :

[0078]

[0079] In the formula, It represents the number of altimeter points included in the ground altimeter strip of the spaceborne lidar. , These are the row numbers of the elevation measurement points in the geographic coordinate system. It is the elevation of the measuring point.

[0080] Step 3: Construct a registration model. This model characterizes the planar positional offset between the ground altimeter strips of the spaceborne lidar and the DEM elevation strips. It includes parameters such as scaling factor, rotation angle, and row and column offsets of the two strips in the geographic coordinate system. The formula is expressed as:

[0081]

[0082] In the formula, , It is the initial row and column number of the ground altimeter strip of the spaceborne lidar in the geographic coordinate system. It is the scaling factor. It is the rotation angle. , It represents the row and column offsets of the ground altimeter strips from the spaceborne lidar relative to the DEM in the geographic coordinate system. , It is the row and column number after the registration parameter transformation of the ground altimeter strip of the spaceborne lidar. , , , Together they constitute the registration parameters to be solved in the registration model.

[0083] Step 4: Based on the registration model with the current registration parameters and the geographical location of the ground altimetry strips of the spaceborne lidar, obtain the elevation strips corresponding to the DEM; then calculate the elevation projection length of the DEM elevation strips along the trajectory of the ground altimetry strips of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system, denoted as the DEM elevation projection length calculated based on the current registration parameters. .

[0084] Set the initial registration parameters before step 4, and then in step 4, based on the current registration parameters... The registration model is used to calculate the geographical location of the spaceborne lidar strip after offset, so as to obtain the corresponding DEM elevation strip. Then, based on the topographic feature information such as the coordinates, elevation, slope, and slope change rate of the DEM in the geographic coordinate system, the elevation change curves of adjacent points of the DEM elevation strip are fitted, and the elevation projection length of the DEM elevation strip along the orbital direction of the spaceborne lidar strip is calculated. Since the spaceborne lidar data acquisition is dense, the topographic contour between adjacent ground points corresponding to the DEM ground point set can be fitted as a quadratic function, expressed as:

[0085]

[0086] These are the row numbers of the DEM ground points in the geographic coordinate system. , , The fitting coefficients are to be determined, representing the coefficients of the quadratic function used to fit the terrain profile between adjacent ground points.

[0087] The length of the projection line of the terrain contour of adjacent ground points onto the O-xz plane is the path integral of this function between the two endpoints, calculated as follows:

[0088]

[0089] In the formula, This represents the path between adjacent points on the DEM ground plane in the O-xz plane. It is the first point in the set. Segment path.

[0090] The rate of change of slope of a DEM ground point along the DEM strip direction can represent the derivative of the fitting function at that point. Therefore, the coordinates, elevations, and rate of change of slope of the terrain contour at both ends of the ground point can be used as observations to construct an expression for solving the fitting coefficients using the least squares adjustment method:

[0091]

[0092] In the formula, It is residual error. It is the overall coefficient matrix. It is the fitting coefficient matrix. This is the observation error vector. For a certain segment of terrain contour, the above matrices are represented as follows:

[0093]

[0094]

[0095]

[0096] In the formula, , These are the row numbers of the terrain outline at the two ends of the ground. , It refers to the elevation of the terrain outline at the two ends of the ground. , It is the rate of change of slope at both ends of the ground along the direction of the DEM strip;

[0097] Based on this, the least squares estimation is used to fit the coefficient matrix. Perform the calculation:

[0098]

[0099] In the formula, This represents the weight matrix, which is an identity matrix with diagonal elements of 1, indicating that the terrain information at both endpoints has equal weight in fitting the function.

[0100] Finally, by calculating the projected lengths of the terrain contour lines of adjacent point sets in the DEM strip, the total projected length of the terrain contour lines of the DEM strip is obtained. :

[0101]

[0102] Step 5: Calculate the elevation projection length of the spaceborne lidar obtained in Step 2. and the DEM elevation projection length obtained in step 4 The relative error between them is used to construct an error function, and the registration parameters in the registration model are iteratively solved by the nonlinear least squares method.

[0103] When transforming the coordinates of the spaceborne lidar data and its corresponding DEM ground points according to the registration parameter vector... The total projected length of the DEM strip terrain contour lines is a fixed value. The registration will change, and when the relative error between the two reaches its minimum, the registration is considered successful. Therefore, the relative error between the total projected length of the topographic contour line of the spaceborne lidar topographic altimeter strip and the corresponding DEM strip can be constructed under the current registration parameters, and this can be used as the error function:

[0104]

[0105] In the formula, This indicates that the maximum value is taken between the projected length of the spaceborne lidar and the DEM strip terrain contour line;

[0106] By incorporating the error function into the registration formula, it can be characterized as a nonlinear model function:

[0107]

[0108] In the formula, It is the first Each observation represents a response variable indicating the initial row and column number of the ground altimeter strip from the spaceborne lidar before registration in the geographic coordinate system. This is the length of the terrain contour line calculated from the spaceborne lidar altimeter data after registration parameter transformation, corresponding to... , ; The length of the terrain contour line corresponding to the DEM elevation strip is obtained through the registration model function, i.e., formula (2), where the independent variable is... It represents the row and column number of the i-th altimeter in the ground altimeter strip of the spaceborne lidar, before registration transformation, in the geographic coordinate system. This is a parameter vector containing scaling factors, rotation angles, and row and column offsets, represented as... .

[0109] The parameter vector is solved using a nonlinear least squares method. Specifically:

[0110] Based on the constructed registration model, calculate the sum of squared residuals:

[0111]

[0112] The optimal solution for the parameter vector is obtained by minimizing the sum of squared residuals.

[0113]

[0114] In the formula, This indicates that the parameters are obtained when the sum of squared residuals is minimized. The optimal solution is: ;

[0115] The optimal solution of registration parameters calculated in step 5 For reference, the relative error of the topographic contour projection line lengths of the mobile spaceborne lidar topographic altimetry strips and the corresponding DEM strips is calculated based on their planar positions. When the relative error is less than a given threshold (e.g., 0.0001), the iteration stops and the optimal solution with the currently calculated registration parameters is obtained. As the output, otherwise you need to return to step 3-5 to recalculate the registration parameter vector.

[0116] Step 6, Constructing the Offset Lookup Table: Based on the registration parameters obtained in Step 5, using the row and column numbers of all points in the ground altimeter strip of the spaceborne lidar in the geographic coordinate system and the offset of the DEM relative to them in the row and column directions as observations, the least squares adjustment method is used to solve the offset polynomials of the DEM relative to the ground altimeter strip of the spaceborne lidar in the row and column directions, respectively, to obtain the offset lookup table of the DEM relative to the spaceborne lidar data, which is used to characterize the degree of planar position offset of each pixel of the DEM in the row and column directions. This includes the following sub-steps:

[0117] Step 6.1: For each ground altimeter strip of the spaceborne lidar in the study area, calculate its registration parameters through steps 2 to 5.

[0118] Step 6.2: The expression for correcting the planar position of the entire DEM using polynomial fitting offset is:

[0119]

[0120] In the formula, and These are the polynomial fitting coefficients. ; and The row and column numbers of the DEM in the geographic coordinate system. and This represents the offset of the DEM relative to the altimeter point of the spaceborne lidar in each ground altimeter strip; for this polynomial, at least 5 spaceborne lidar altimeter strips are required to solve its coefficients.

[0121] Step 6.3: Using the row and column numbers of each altimeter point in the ground altimeter strip of the spaceborne lidar in the geographic coordinate system... and the row and column offsets in the registration parameter vector calculated in step 6.1. , As observed values, the polynomial coefficients in equation (12) are solved, and the fitting coefficients are known. The expression (12) is used to construct the offset lookup table of DEM relative to the spaceborne lidar data. Figure 3 The table shows the offset lookup values ​​for two DEMs of the area under study in the row and column directions, respectively. The size of the lookup table is the same as that of the DEM, and it stores the planar offset of each pixel of the DEM relative to the spaceborne lidar data.

[0122] Step 7: Using the offset lookup table constructed in Step 6, perform a geometric transformation on the DEM, transforming its coordinates to a new coordinate position in the geographic coordinate system, and output the final DEM, completing the DEM planar position correction. Figure 4 (a) and Figure 4As shown in (d). Using the airborne DTM as a reference, the InSAR DEM before and after the planar position shift is subtracted from the DTM, respectively. Figure 4 (b) and Figure 4 (e) shows the differences between the two DEM scenes and the airborne DTM before the planar position offset. Figure 4 (c) and Figure 4 (f) shows the differences between the two DEMs after the planar position shift and the airborne DTM.

[0123] Based on the above analysis, this embodiment uses spaceborne ICESat-2 ATL03 data as the ground altimetry data for the spaceborne lidar. It has dense sampling along the orbital direction, making it suitable for registration in strip units. The experimental area DEM estimated using spaceborne ICESat-2 ATL03 data and LT-1 InSAR data is registered, and planar position correction is performed to generate the final DEM of the experimental area. Figure 4 (a) and Figure 4 As shown in (d). Simultaneously, using the airborne DTM as a reference, comparisons were made with the DEM before and after planar position correction. Among them, Figure 4 (b) and Figure 4 (e) shows the differences between the two DEM scenes and the airborne DTM before the planar position offset. Figure 4 (c) and Figure 4 (f) shows the differences between the two DEM images after planar position offset and the airborne DTM. It can be seen that before planar position correction, there was a significant planar position deviation between the InSAR DEM and the airborne DTM. After registration and offset correction, the planar position misalignment between the two was effectively eliminated. In terms of accuracy, this invention achieves a level of accuracy comparable to existing methods without relying on existing terrain products. It is an effective method for autonomous, controllable, efficient, and high-precision DEM planar position correction.

[0124] The above embodiments are preferred embodiments of this application. Those skilled in the art can make various changes or improvements based on them. Without departing from the overall concept of this application, these changes or improvements should fall within the scope of protection claimed in this application.

Claims

1. A method for correcting the planar position of an InSAR DEM based on lidar data, characterized in that, Includes the following steps: Step 1: Acquire InSAR image data of the area to be studied and perform interferometric processing to generate an InSAR DEM; Step 2: Acquire ground altimetry data from multiple co-tracks of the spaceborne lidar in the area to be studied. Each co-track ground altimetry data is denoted as one lidar ground altimetry strip. Calculate the elevation projection length of each spaceborne lidar ground altimetry strip along the track direction on the O-xz plane of the spatial rectangular coordinate system, and denote it as the spaceborne lidar elevation projection length. ; Step 3: Construct a registration model between the ground altimeter strips of the spaceborne lidar and the elevation strips of the DEM; wherein, the registration parameters of the registration model include the row and column offsets of the ground altimeter strips of the spaceborne lidar relative to the DEM in the geographic coordinate system, the rotation angle and the scaling factor; set the initial registration parameters; Step 4: Based on the registration model with the current registration parameters and the geographical location of the ground altimetry strips of the spaceborne lidar, obtain the elevation strips corresponding to the DEM; then calculate the elevation projection length of the DEM elevation strips along the trajectory of the ground altimetry strips of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system, denoted as the DEM elevation projection length calculated based on the current registration parameters. ; Step 5: Calculate the elevation projection length of the spaceborne lidar obtained in Step 2. and the DEM elevation projection length obtained in step 4 The relative error between them is used to construct an error function, and the registration parameters in the registration model are iteratively solved by the nonlinear least squares method. Step 6: Based on the registration parameters obtained in Step 5, take the row and column numbers of all points in the ground altimetry strip of the spaceborne lidar in the geographic coordinate system and the offset of the DEM relative to them in the row and column directions as observations, and use the least squares adjustment method to solve the offset polynomials of the DEM relative to the ground altimetry strip of the spaceborne lidar in the row and column directions respectively, and obtain the offset lookup table of the DEM relative to the spaceborne lidar data. Step 7: Using the offset lookup table obtained in Step 6, transform the initial DEM coordinates obtained in Step 1 to the new coordinate position to obtain the DEM after planar position correction.

2. The method according to claim 1, characterized in that, Calculate the elevation projection length of each spaceborne lidar ground altimeter strip along its orbit on the O-xz plane of the Cartesian coordinate system. The calculation formula is: ； In the formula, It represents the number of altimeter points included in the ground altimeter strip of the spaceborne lidar. These are the row numbers of the elevation measurement points in the geographic coordinate system. It is the elevation of the measuring point.

3. The method according to claim 1, characterized in that, Step 3: Based on the planar positional deviation between the ground altimetry strips of the spaceborne lidar and the DEM elevation strips, the following registration model is constructed: ； In the formula, , It is the initial row and column number of the ground altimeter strip of the spaceborne lidar in the geographic coordinate system. It is the scaling factor. It is the rotation angle. , It represents the row and column offsets of the ground altimetry strips from the spaceborne lidar relative to the DEM elevation strips in the geographic coordinate system. , , , These are all registration parameters to be solved; , It is the row and column number after the registration parameter transformation of the ground altimeter strip of the spaceborne lidar.

4. The method according to claim 1, characterized in that, Step 4: Calculate the elevation projection length of the DEM elevation strip along the trajectory of the ground altimeter strip of the spaceborne lidar on the O-xz plane of the spatial rectangular coordinate system. Specifically, it includes: First, the terrain contour lines between adjacent ground points of the DEM are fitted with elevation change curves in the geographic coordinate system, as follows: ； In the formula, These are the row numbers of the DEM ground points in the geographic coordinate system. , , These are the fitting coefficients to be determined. Then, the derivative of the elevation change curve at each surface point is expressed as the slope change rate of each surface point along the strip direction in the DEM. Furthermore, using the coordinates, elevation, and slope change rate of the terrain contour line at both ends of the ground as observed values, an expression for solving the fitting coefficients using the least squares adjustment method is constructed: ； In the formula, It is residual error. It is the overall coefficient matrix. It is the fitting coefficient matrix. These are the observation error vectors, represented as follows: ； ； ； In the formula, , These are the row numbers of the terrain outline at the two ends of the ground. , It refers to the elevation of the terrain outline at the two ends of the ground. , It is the rate of change of slope of the two ground points along the DEM elevation strip direction; Then, the least squares estimation was used to fit the coefficient matrix. Perform the calculation: ； In the formula, This represents the weight matrix, which is an identity matrix with diagonal elements all being 1, indicating that the terrain information from the two ends of the ground has equal weight in fitting the function. Then based on the fitting coefficients , , To solve for the given terrain contour, calculate the projection lengths of the terrain contour onto the O-xz plane between adjacent ground points by integration: ； In the formula, The path representing the terrain outline, For the first i The projection length of the terrain contour line onto the O-xz plane; Finally, the total projected length of the terrain contour lines between all adjacent ground points in the DEM elevation strip is obtained by summing the projected lengths of the terrain contour lines in the DEM elevation strip. : 。 5. The method according to claim 1, characterized in that, Step 5 specifically includes: First, construct the error function: ； In the formula, Indicates taking The maximum value in; The radar elevation projection length calculated in step 2 As an observation, the DEM elevation projection length calculated by the registered model will be... Represented as a nonlinear function Then the error function is transformed into: ； Among them, the independent variable of the nonlinear function It is a height measurement point in the ground height measurement strip of the spaceborne lidar. i Initial row and column numbers in the geographic coordinate system; The registration parameters to be solved include scaling factors. Rotation angle and row and column offsets , ; Initial row and column numbers The corresponding elevation projection length error; The nonlinear least squares method is used to solve for the registration parameters in the error function by minimizing the sum of squared errors, thus obtaining the optimal solution for the registration parameters. : ； In the formula, It represents the number of altimeter points included in the ground altimeter strip of the spaceborne lidar.

6. The method according to claim 1, characterized in that, Step 6 includes the following specific steps: Step 6.1: For each ground altimeter strip of all spaceborne lidars in the study area, run steps 2 to 5 to obtain the registration parameter vectors for different strips. Step 6.2: Construct a polynomial fitting offset to correct the planar position of the entire DEM. The expression is as follows: ； In the formula, and These are the polynomial fitting coefficients. ; and The row and column numbers of the DEM in the geographic coordinate system. and This represents the offset of the DEM relative to the altimeter points of the spaceborne lidar in each ground altimeter strip; Step 6.3: Calculate the row and column numbers of each altimeter point in the ground altimeter strip from the spaceborne lidar in the geographic coordinate system. and the row and column offsets in the registration parameter vector obtained in step 6.

1. , Solve the polynomial fitting coefficients in equation (11) Given the fitting coefficient The expression (11) is the offset lookup table of DEM relative to the spaceborne lidar data.

7. An inversion device for correcting the planar position of an InSAR DEM based on lidar data, characterized in that, include: The preprocessing module is used to perform interferometric processing on the acquired InSAR data of the area to be studied to obtain the initial DEM of the image. It is also used to: calculate the elevation projection length along the orbital direction on the O-xz plane of the spatial rectangular coordinate system for each ground altimeter strip acquired by the spaceborne lidar in the area to be studied, denoted as the radar elevation projection length. ; The registration module is used to: construct a registration model to register the DEM and spaceborne lidar data of the area under study according to the principle of terrain contour consistency, and obtain the planar offset of the DEM relative to the spaceborne lidar data in the geographic coordinate system. The registration parameter calculation method for the registration model is as follows: Based on the geographical location of the radar ground altimetry strip, obtain the elevation strip corresponding to the DEM; then calculate the elevation projection length of the DEM elevation strip along the orbital direction of the spaceborne lidar ground altimetry strip on the O-xz plane of the spatial rectangular coordinate system. And then according to and The relative error between them is calculated and the registration parameters in the registration model are iteratively solved using a nonlinear least squares method; The lookup table construction module is used to: take the row and column numbers of all points in the ground altimetry strip of the spaceborne lidar in the geographic coordinate system and the offset of the DEM relative to them in the row and column directions as observations, and use the least squares adjustment method to solve the offset polynomials of the DEM relative to the ground altimetry strip of the spaceborne lidar in the row and column directions respectively, and obtain the offset lookup table of the DEM relative to the spaceborne lidar data. The DEM correction module is used to: use the offset lookup table obtained by the lookup table construction module to transform the initial DEM coordinates obtained by the preprocessing module to a new coordinate position, and obtain the DEM after planar position correction.

8. An electronic device, comprising a memory and a processor, characterized in that, The memory stores a computer program, and when the processor executes the computer program, the processor can implement the method as described in any one of claims 1 to 6.

9. A computer-readable storage medium, characterized in that, It stores a computer program that, when executed by a processor, implements the method as described in any one of claims 1 to 6.

Images