Method, device and equipment for correcting InSAR DEM (Digital Elevation Model) plane position based on laser radar data and medium

Abstract

The invention discloses a method, a device, equipment and a medium for correcting an InSAR DEM plane position based on laser radar data, and the method comprises the steps: obtaining InSAR data, carrying out the interference processing, and generating a DEM; acquiring satellite-borne laser radar data, combining with a ground point corresponding to the DEM, taking a satellite-borne laser radar common-rail strip as a unit, and performing registration based on a relative error of the projection lengths of the topographic contours of the satellite-borne laser radar data and the ground point to obtain an offset of the DEM relative to the satellite-borne laser radar data; constructing offset polynomial models of the DEM in longitude and latitude directions by using a least square adjustment method based on the registration plane offset; and geometric transformation is carried out on the DEM through an offset polynomial model, and finally plane position correction is realized. According to the method, the plane position of the DEM is corrected only through satellite-borne laser radar data, and the problem that an existing method depends on existing terrain products is solved.

Description

Technical Field

[0001] The present invention belongs to the field of geodetic measurement of synthetic aperture radar interferometry (InSAR), and in particular relates to a method, a device, a equipment and a medium for correcting the plane position of InSAR DEM based on laser radar data. Background Art

[0002] Digital Elevation Model (DEM) is one of the basic data sources for surveying and mapping and geographic information construction, and has important applications in infrastructure construction, disaster monitoring, agriculture and forestry, etc. Synthetic Aperture Radar (InSAR) technology has become the main technical means for DEM production due to its all-day, all-weather, large-scale, and high-resolution imaging characteristics.

[0003] At present, the commonly used InSAR technology to produce DEM is based on the DEM estimation of differential interferometry (D-InSAR), the core of which is to use the existing external DEM to simulate the phase, and then differentiate it from the InSAR interferometry phase, and then unwrap and transform the differential phase and superimpose it on the external DEM to generate a new DEM in the SAR coordinate system, and then rely on the external DEM for geocoding to the geographic coordinate system to generate the DEM. However, the geocoding process of converting DEM from the SAR coordinate system to the geographic coordinate system depends on the positioning and registration of the SAR image and the external DEM, which is usually affected by the SAR image orbit error and the DEM accuracy, and it is easy to transfer the trend plane error of the external DEM to the DEM estimated by InSAR.

[0004] On the other hand, satellite-borne lidar technology measures the height of the earth's surface by emitting laser pulses to the surface and receiving echo signals. It can densely collect height measurement points on a global scale and has high planar and elevation accuracy, and can be widely used in InSAR DEM control data.

[0005] In summary, the existing InSAR DEM plane positioning method is heavily dependent on the existing terrain products. Therefore, it is necessary to design a high-precision DEM plane position correction method that does not rely on the existing external DEM products and is only based on space-borne lidar data. Summary of the invention

[0006] The present invention provides a method, device, equipment and medium for correcting InSAR DEM plane position based on laser radar data, which can realize high-precision DEM plane position correction without relying on existing external DEM products.

[0007] In order to achieve the above technical objectives, the present invention adopts the following technical solutions: A method for estimating DEM based on InSAR and LiDAR data, comprising the following steps: Step 1: Obtain InSAR image data of the area to be studied and perform interferometric processing to generate InSAR DEM; Step 2: Obtain multiple ground altimetry data of the satellite-borne laser radar in the area to be studied, and each ground altimetry data in the same track is recorded as a radar ground altimetry strip; calculate the elevation projection length of each satellite-borne laser radar ground altimetry strip along the track direction on the O-xz plane of the spatial rectangular coordinate system, and record it as the satellite-borne laser radar elevation projection length ; Step 3, constructing a registration model between the satellite-borne laser radar ground altimetry strip and the DEM elevation strip; wherein the registration parameters of the registration model include the row and column offsets of the satellite-borne laser radar ground altimetry strip relative to the DEM in the geographic coordinate system; setting the initial registration parameters; Step 4: According to the registration model of the current registration parameters and the geographical location of the ground altimetry strip of the satellite-borne laser radar, 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 ground altimetry strip of the satellite-borne laser radar on the O-xz plane of the spatial rectangular coordinate system, which is recorded as the DEM elevation projection length calculated based on the current registration parameters. ; Step 5: Calculate the projection length of the satellite-borne laser radar elevation 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, according to the registration parameters calculated in step 5, the row and column numbers of all points in the satellite-borne laser radar ground altimetry strip in the geographic coordinate system and the offset of DEM relative to them in the row and column directions are used as observation quantities, and the least squares adjustment method is used to solve the offset polynomials of DEM relative to the satellite-borne laser radar ground altimetry strip in the row and column directions respectively, and obtain the offset lookup table of DEM relative to the satellite-borne laser radar data; Step 7: Using the offset lookup table obtained in step 6, the coordinates of the initial DEM obtained in step 1 are transformed to new coordinate positions to obtain a DEM after plane position correction.

[0008] Furthermore, the elevation projection length of each satellite-borne laser radar ground height measurement strip along the track direction on the O-xz plane of the spatial rectangular coordinate system is calculated. , the calculation formula is:

[0009] In the formula, is the number of altimeter points contained in the ground altimeter strip of the satellite-borne laser radar, x is the row number in the geographic coordinate system, It is the elevation of the measuring point.

[0010] Furthermore, in step 3, the following registration model is constructed based on the plane position deviation between the ground altimetry strips of the satellite-borne lidar and the DEM elevation strips:

[0011] In the formula, , is the initial row and column number of the ground altimetry strip of the spaceborne laser radar in the geographic coordinate system, is the scaling factor, is the rotation angle, , It is the row and column offset of the satellite-borne lidar ground height measurement strip relative to the DEM height strip in the geographic coordinate system. , , , These are the registration parameters to be solved; , It is the row and column number of the ground altimetry strip of the satellite-borne laser radar after the registration parameters are transformed.

[0012] Furthermore, step 4 calculates the elevation projection length of the DEM elevation strip along the orbital direction of the satellite-borne lidar ground height measurement strip on the O-xz plane of the spatial rectangular coordinate system. , specifically including: First, the terrain contour lines between adjacent ground points of DEM are fitted with elevation change curves in the geographic coordinate system, which can be expressed as:

[0013] In the formula, It is the row number of the DEM ground point in the geographic coordinate system. , , is the desired fitting coefficient, 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, and then the coordinates, elevation and slope change rate of the ground points at both ends of the terrain contour line are used as observation values ​​to construct an expression for solving the fitting coefficient using the least squares adjustment method:

[0014] In the formula, is the residual error, is the row number matrix of the ground points at both ends, is the fitting coefficient matrix, is the observation error vector, which can be expressed as:

[0015]

[0016]

[0017] In the formula, , is the row number of the ground points at both ends of the terrain contour line, , is the elevation of the terrain contour line at the ground points at both ends. , It is the slope change rate of the ground points at both ends along the DEM elevation strip direction; Then use the least squares estimation to fit the coefficient matrix Solve:

[0018] In the formula, represents the weight matrix, whose value is the unit matrix with diagonal elements of 1, indicating that the weights of the terrain information of the ground points at both ends for the function fitting are equal; According to the fitting coefficient , , Solve for known terrain contour lines and calculate the projection line length of the terrain contour line between adjacent ground points on the O-xz plane by integration:

[0019] In the formula, represents the path of terrain contours, For the i The projection length of the topographic contour line on the O-xz plane; Finally, the projection lengths of the terrain contour lines between all adjacent ground points in the DEM elevation strip are summed to obtain the total projection length of the terrain contour lines in the DEM elevation strip.

[0020]

[0021] Furthermore, step 5 specifically includes: First, construct the error function:

[0022] In the formula, Indicates taking The maximum value in ; The radar elevation projection length calculated in step 2 As the observation, the DEM elevation projection length calculated by the registration model is It can be expressed as a nonlinear function: , then the error function is transformed into:

[0023] Among them, the nonlinear function independent variable It is a ground point in the ground altimetry strip of the satellite-borne laser radar. i The initial row and column number in the geographic coordinate system; are the registration parameters to be solved, including the scaling factor , rotation angle And row and column offsets , ; Initial row and column number The corresponding elevation projection length error; The nonlinear least squares method is used to solve the registration parameters in the error function by minimizing the sum of squared errors to obtain the optimal solution for the registration parameters. :

[0024] In the formula, It is the number of altimeter points contained in the ground altimeter strip of the spaceborne lidar.

[0025] The specific process of step 6 includes: Step 6.1: For each ground height measurement strip of all spaceborne lidars in the study area, run steps 2 to 5 to obtain the registration parameter vectors of different strips; Step 6.2, construct a polynomial fitting offset to correct the plane position of the entire scene DEM, and its expression is:

[0026] In the formula, and are the polynomial fitting coefficients, ; and is the row and column number of DEM in the geographic coordinate system, and is the offset of DEM in each ground altimetry strip relative to the altimetry point of the satellite-borne lidar; Step 6.3, the row and column number of the ground height measurement strip of the satellite-borne laser radar 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) , the known fitting coefficients The expression (11) is the offset lookup table of DEM relative to the space-borne lidar data.

[0027] An inversion device for correcting the plane position of InSAR DEM based on laser radar data, comprising: The preprocessing module is used to perform interferometric processing on the 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 of each ground height measurement strip of the satellite-borne laser radar in the area to be studied in the along-track direction on the O-xz plane of the spatial rectangular coordinate system, which is recorded as the satellite-borne laser radar elevation projection length ; The registration module is used to: construct a registration model to register the DEM and satellite-borne laser radar data of the study area according to the principle of terrain contour consistency, and obtain the plane offset of the DEM relative to the satellite-borne laser radar data in the geographic coordinate system; The registration parameter solution method of the registration model is as follows: obtain the elevation strip corresponding to the DEM according to the geographical location of the satellite-borne laser radar ground altimetry strip, and then calculate the elevation projection length of the DEM elevation strip along the orbit direction of the satellite-borne laser radar 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 by the nonlinear least square method; The lookup table construction module is used to: take the row and column numbers of all points in the satellite-borne laser radar ground altimetry strip in the geographic coordinate system and the offset of DEM relative to them in the row and column directions as observation quantities, use the least squares adjustment method to solve the offset polynomials of DEM relative to the satellite-borne laser radar ground altimetry strip in the row and column directions respectively, and obtain the offset lookup table of DEM relative to the satellite-borne laser radar data; The DEM correction module is used to: transform the initial DEM coordinates obtained in the preprocessing module to new coordinate positions using the offset lookup table obtained by the lookup table construction module to obtain the DEM after plane position correction.

[0028] An electronic device comprises a memory and a processor, wherein a computer program is stored in the memory, and when the processor executes the computer program, the processor can implement any of the above-mentioned methods for correcting the plane position of an InSAR DEM based on lidar data.

[0029] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements any of the above-mentioned methods for correcting the plane position of an InSAR DEM based on lidar data.

[0030] Beneficial Effects The present invention corrects InSAR based on satellite-borne laser radar data The method of DEM plane position realizes the correction of DEM plane position based only on satellite-borne lidar data: first, the InSAR data is processed by the commonly used conventional data processing method to obtain the initial DEM; then the satellite-borne lidar terrain altimetry data is used to calculate the projection length of its elevation terrain contour line in units of strips; then the initial registration parameters are set, and the satellite-borne lidar terrain altimetry data is moved according to the registration formula to obtain the corresponding DEM strip and fit the terrain contour of the adjacent point set of the DEM strip, and the fitting function coefficient is solved according to the elevation and slope change rate, and the projection length of the DEM strip terrain contour is solved by integration; then the relative error degree between the satellite-borne lidar terrain altimetry data and the corresponding DEM strip terrain contour projection length is calculated, and used as the error function to assist in iteratively solving the registration parameter vector; then, according to the registration parameters of several satellite-borne lidar terrain altimetry strips in the area to be studied, the offset polynomial of DEM relative to the satellite-borne lidar data is fitted, and an offset lookup table is constructed; finally, the DEM is geometrically transformed through the offset lookup table, its coordinates are transformed to the new geographical location in the geographic coordinate system, the final DEM result is estimated, and the DEM plane position correction is completed.

[0031] The beneficial effects of this method are as follows: a method system for correcting the DEM plane position based only on space-borne laser radar data is constructed, overcoming the problem that the existing InSAR DEM geocoding is heavily dependent on existing terrain products; a method for aligning space-borne laser radar data with DEM based on the relative error of terrain contour is proposed, realizing the correction of DEM plane errors. In terms of accuracy, the present invention achieves an accuracy level and operating efficiency comparable to that of existing methods without relying on existing terrain products, and is an effective method for autonomous, controllable, efficient and high-precision DEM plane position correction. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 The present invention is a flowchart of the method described in the embodiment of the present invention.

[0033] Figure 2 The location of the verification area and the coverage of the DEM selected for the present invention.

[0034] Figure 3 It is a lookup table for the offset of 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 of DEM1 in the column direction; Figure 3(b) is the offset lookup table of DEM1 in the row direction; Figure 3 (c) is the offset lookup table of DEM2 in the column direction; Figure 3 (d) is the offset lookup table of DEM2 in the row direction.

[0035] Figure 4 These are the estimated DEM after plane position correction and the difference map between DEM and airborne DTM before and after plane position correction. Figure 4 (a) is the estimated DEM1 after plane position correction; Figure 4 (b) is the difference map between DEM1 and airborne DTM before plane position correction; Figure 4 (c) is the difference map between DEM1 and airborne DTM after plane position correction; Figure 4 (d) is the estimated DEM2 after plane position correction; Figure 4 (e) is the difference map between DEM2 and airborne DTM before plane position correction; Figure 4 (f) is the difference map between DEM2 and airborne DTM after plane position correction. DETAILED DESCRIPTION

[0036] The following is a detailed description of an embodiment of the present invention. This embodiment is based on the technical solution of the present invention, and provides a detailed implementation method and a specific operation process to further explain the technical solution of the present invention.

[0037] In order 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 Lu Tan 1 (LT-1) InSAR data and satellite-borne ICESat-2 ATL03 data; it is noted that the specific implementation described herein is only used to explain the present invention and is not intended to limit the present invention.

[0038] This embodiment provides a method for correcting the plane position of InSAR DEM based on satellite-borne laser radar data. Figure 1 As shown, the following steps are included: Step 1: According to the coordinate range of the selected experimental area, LT-1 InSAR data is acquired, the acquired data is interferometrically processed, and the InSAR DEM is generated by interferometric processing.

[0039] This implementation example selects 2 scenes of LT-1 data located in the experimental area in central Spain. The experimental area is located as follows: Figure 2As shown; the data is interferometrically processed to generate the InSAR DEM of the area to be studied. Interferometric processing technology is a prior art method, and reference may be made to the literature: 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 be described in detail.

[0040] Step 2: According to the coordinate range of the selected experimental area, the satellite-borne lidar data is obtained and projected from the spatial rectangular coordinate system to the same acquisition track (i.e., the same track strip) of the satellite flying over once. plane, and calculate the projection length of the terrain contour line :

[0041] In the formula, is the number of altimeter points contained in the ground altimeter strip of the satellite-borne laser radar, , is the row number of the height measurement point in the geographic coordinate system, It is the elevation of the measuring point.

[0042] Step 3: Construct a registration model. This model characterizes the plane position offset between the satellite-borne laser radar ground altimetry strip and the DEM elevation strip. It includes parameters such as scaling factor, rotation angle, and row and column offset between the satellite-borne laser radar ground altimetry strip and the DEM elevation strip in the geographic coordinate system. The formula is:

[0043] In the formula, , is the initial row and column number of the ground altimetry strip of the spaceborne laser radar in the geographic coordinate system, is the scaling factor, is the rotation angle, , It is the row and column offset of the satellite-borne laser radar ground altimetry strip relative to the DEM in the geographic coordinate system. , It is the row and column number of the ground altimetry strip of the satellite-borne laser radar after the registration parameters are transformed. , , , Together they constitute the registration parameters to be solved for the registration model.

[0044] Step 4: According to the registration model of the current registration parameters and the geographical location of the ground altimetry strip of the satellite-borne laser radar, 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 ground altimetry strip of the satellite-borne laser radar on the O-xz plane of the spatial rectangular coordinate system, which is recorded as the DEM elevation projection length calculated based on the current registration parameters. .

[0045] 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 geographic location of the satellite-borne laser radar strip after the offset to obtain the corresponding DEM elevation strip. Then, the elevation change curve of the adjacent points of the DEM elevation strip is fitted according to the terrain feature information such as the coordinates, elevation, slope and slope change rate of the DEM in the geographic coordinate system, and the elevation projection length of the DEM elevation strip along the orbit direction of the satellite-borne laser radar strip is calculated. Due to the intensive collection of satellite-borne laser radar data, the terrain contour between adjacent ground points corresponding to the ground point set of the DEM can be fitted as a quadratic function, expressed as:

[0046] It is the row number of the DEM ground point in the geographic coordinate system. , , is the fitting coefficient to be determined, which represents the coefficient of the quadratic function used to fit the terrain contour between adjacent ground points.

[0047] The length of the projection line of the terrain contour of adjacent ground points on the O-xz plane is the path integral of the function between the two endpoints, and the calculation formula is:

[0048] In the formula, Represents the path between adjacent points of the DEM ground point in the O-xz plane. It is the point concentration Segment path.

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

[0050] In the formula, is the residual error, is the overall coefficient matrix, is the fitting coefficient matrix, is the observation error vector. For a certain section of terrain contour line, the above matrices are expressed as:

[0051]

[0052]

[0053] In the formula, , is the row number of the ground points at both ends of the terrain contour line, , is the elevation of the terrain contour line at the ground points at both ends. , It is the slope change rate of the ground points at both ends along the DEM strip direction; On this basis, the least squares estimation is used to estimate the fitting coefficient matrix Solve:

[0054] In the formula, Represents the weight matrix, whose value is the identity matrix with diagonal elements of 1, indicating that the terrain information at the two endpoints has equal weight for function fitting.

[0055] Finally, after calculating the projection lengths of the terrain contour lines of adjacent point sets in the DEM strips one by one, the total projection length of the terrain contour lines of the DEM strips is obtained. :

[0056] Step 5: Calculate the projection length of the satellite-borne laser radar elevation obtained in step 2 And the DEM elevation projection length obtained in step 4 The relative error between them is used to construct the error function, and the registration parameters in the registration model are iteratively solved by the nonlinear least squares method.

[0057] When transforming the coordinates of the satellite-borne lidar data and its corresponding DEM ground points according to the registration parameter vector, is a fixed value, and the total length of the projection of the DEM strip terrain contour line will change. When the relative error between the two reaches the minimum, the registration is considered successful. Therefore, the relative error between the total length of the projection of the satellite-borne lidar terrain height measurement strip and the corresponding DEM strip terrain contour line under the current registration parameter conditions can be constructed and used as the error function:

[0058] In the formula, It means taking the maximum value between the projection length of the satellite-borne laser radar and the DEM strip terrain contour line; Introducing the error function into the registration formula, it can be characterized as a nonlinear model function:

[0059] In the formula, It is observations, expressed as the response variable of the initial row and column number of the satellite-borne laser radar ground altimetry strip in the geographic coordinate system before registration, that is, the length of the terrain contour line calculated after the satellite-borne laser radar altimetry data is transformed by the registration parameters, corresponding to , ; is the length of the terrain contour line corresponding to the DEM elevation strip obtained by the registration model function, that is, formula (2), where the independent variable is the row and column number of the i-th altimeter point in the ground altimeter strip of the satellite-borne laser radar in the geographic coordinate system before the registration transformation, is a parameter vector, including scaling factor, rotation angle, row and column offset, expressed as .

[0060] The parameter vector is solved using the nonlinear least squares method. , specifically: According to the constructed registration model, the residual sum of squares is calculated:

[0061] By minimizing the residual sum of squares, the optimal solution of the parameter vector is obtained:

[0062] In the formula, Indicates that when the residual sum of squares reaches the minimum, the parameter is solved The optimal solution of ; The optimal solution of the registration parameters calculated in step 5 As a reference, the plane position of the moving satellite-borne lidar terrain height measurement strip and the corresponding DEM strip are calculated to calculate the relative error of the length of the terrain contour projection line between the two, that is, When the relative error is less than a given threshold (e.g. 0.0001), the iteration stops and the optimal solution is obtained using the currently calculated registration parameters. As the output result, otherwise it is necessary to return to steps 3-5 to calculate the registration parameter vector again.

[0063] Step 6, construct an offset lookup table: according to the registration parameters calculated in step 5, the row and column numbers of all points in the satellite-borne laser radar ground altimetry strip in the geographic coordinate system and the offset of DEM relative to them in the row and column directions are used as observation quantities, and the least squares adjustment method is used to solve the offset polynomial of DEM relative to the satellite-borne laser radar ground altimetry strip in the row and column directions respectively, and the offset lookup table of DEM relative to the satellite-borne laser radar data is obtained, which is used to characterize the plane position offset degree of each pixel of DEM in the row and column directions. It includes the following sub-steps: Step 6.1: For each ground height measurement strip of the spaceborne lidar in the study area, calculate its registration parameters through steps 2 to 5; Step 6.2: The expression for correcting the plane position of the entire DEM using the polynomial fitting offset is:

[0064] In the formula, and are the polynomial fitting coefficients, ; and is the row and column number of DEM in the geographic coordinate system, and is the offset of DEM in each ground altimetry strip relative to the satellite-borne lidar altimetry point; for this polynomial, at least 5 satellite-borne lidar altimetry strips are required to solve its coefficients.

[0065] Step 6.3, use the row and column numbers of each altimeter point in the ground altimeter strip of the satellite-borne laser radar in the geographic coordinate system , and the row and column offsets in the registration parameter vector solved in step 6.1 , As the observed value, solve the polynomial coefficients in equation (12), and the fitting coefficients are known The expression (12) is used to construct the offset lookup table of DEM relative to the satellite-borne lidar data. Figure 3 The offset lookup table of the two DEM scenes in the area to be studied in the row and column directions is shown. The size of the lookup table is the same as the DEM size, and it stores the plane offset of each pixel of the DEM relative to the space-borne lidar data; Step 7, use the offset lookup table constructed in step 6 to perform geometric transformation on the DEM, transform its coordinates to the new coordinate position in the geographic coordinate system, output the final DEM, and complete the DEM plane position correction, such as Figure 4 (a) and Figure 4 (d) Using the airborne DTM as a reference, the InSAR DEM before and after the plane position shift is subtracted from the DTM. Figure 4 (b) and Figure 4(e) shows the difference between the two DEMs and the airborne DTM before the plane position shift. Figure 4 (c) and Figure 4 (f) shows the difference between the two DEMs and the airborne DTM after the plane position is offset.

[0066] Based on the above analysis, this embodiment uses the satellite-borne ICESat-2 ATL03 data as the satellite-borne laser radar ground altimetry data. It has dense sampling in the along-track direction and is suitable for registration in strips. The satellite-borne ICESat-2 ATL03 data is used to align with the DEM of the experimental area estimated by LT-1 InSAR data, and the plane position correction is completed to generate the final DEM of the experimental area, as shown in Figure 2. Figure 4 (a) and Figure 4 (d) At the same time, the airborne DTM is used as a reference to compare with the DEM before and after plane position correction. Figure 4 (b) and Figure 4 (e) shows the difference between the two DEMs and the airborne DTM before the plane position shift. Figure 4 (c) and Figure 4 (f) shows the difference between the two DEMs and the airborne DTM after the plane position offset. It can be seen that before the plane position correction, there is an obvious plane position deviation between the InSAR DEM and the airborne DTM. After the registration and offset correction, the plane position deviation between the two is well eliminated. In terms of accuracy, the present invention achieves an accuracy level comparable to that of existing methods without relying on existing terrain products. It is an effective method for autonomous, controllable, efficient and high-precision DEM plane position correction.

[0067] The above embodiments are preferred embodiments of the present application. Ordinary technicians in this field can also make various changes or improvements on this basis. Without departing from the overall concept of the present application, these changes or improvements should fall within the scope of protection required by the present application.

Claims

1. A method for correcting the plane position of InSAR DEM based on lidar data, characterized in that: The following steps are involved: Step 1: Obtain InSAR image data of the area to be studied and perform interferometric processing to generate InSAR DEM; Step 2: Obtain multiple ground altimetry data of the satellite-borne laser radar in the area to be studied, and each ground altimetry data in the same track is recorded as a radar ground altimetry strip; calculate the elevation projection length of each satellite-borne laser radar ground altimetry strip along the track direction on the O-xz plane of the spatial rectangular coordinate system, and record it as the satellite-borne laser radar elevation projection length ; Step 3, construct a registration model between the satellite-borne laser radar ground altimetry strip and the DEM elevation strip; wherein the registration parameters of the registration model include the row and column offsets, rotation angles and scaling factors of the satellite-borne laser radar ground altimetry strip relative to the DEM in the geographic coordinate system; set the initial registration parameters; Step 4: According to the registration model of the current registration parameters and the geographical location of the ground altimetry strip of the satellite-borne laser radar, 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 ground altimetry strip of the satellite-borne laser radar on the O-xz plane of the spatial rectangular coordinate system, which is recorded as the DEM elevation projection length calculated based on the current registration parameters. ; Step 5: Calculate the projection length of the satellite-borne laser radar elevation 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, according to the registration parameters calculated in step 5, the row and column numbers of all points in the satellite-borne laser radar ground altimetry strip in the geographic coordinate system and the offset of DEM relative to them in the row and column directions are used as observation quantities, and the least squares adjustment method is used to solve the offset polynomials of DEM relative to the satellite-borne laser radar ground altimetry strip in the row and column directions respectively, and obtain the offset lookup table of DEM relative to the satellite-borne laser radar data; Step 7: Using the offset lookup table obtained in step 6, the initial DEM coordinates obtained in step 1 are transformed to new coordinate positions to obtain a DEM after plane position correction.

2. The method according to claim 1, characterized in that: Calculate the elevation projection length of each satellite-borne lidar ground height measurement strip along the track direction on the O-xz plane of the spatial rectangular coordinate system , the calculation formula is: ； In the formula, is the number of altimeter points contained in the ground altimeter strip of the satellite-borne laser radar, is the row number of the height measurement point 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 plane position deviation between the ground altimetry strips of the satellite-borne lidar and the DEM elevation strips, the following registration model is constructed: ； In the formula, , is the initial row and column number of the ground altimetry strip of the spaceborne laser radar in the geographic coordinate system, is the scaling factor, is the rotation angle, , It is the row and column offset of the satellite-borne lidar ground height measurement strip relative to the DEM height strip in the geographic coordinate system. , , , These are the registration parameters to be solved; , It is the row and column number of the ground altimetry strip of the satellite-borne laser radar after the registration parameters are transformed.

4. The method according to claim 1, characterized in that: Step 4: Calculate the elevation projection length of the DEM elevation strip along the orbital direction of the satellite-borne lidar ground height measurement strip on the O-xz plane of the spatial rectangular coordinate system. , specifically including: First, the terrain contour lines between adjacent ground points of DEM are fitted with elevation change curves in the geographic coordinate system, which can be expressed as: ； In the formula, It is the row number of the DEM ground point in the geographic coordinate system. , , is the desired fitting coefficient, 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, and then the coordinates, elevation and slope change rate of the ground points at both ends of the terrain contour line are used as observation values ​​to construct an expression for solving the fitting coefficient using the least squares adjustment method: ； In the formula, is the residual error, is the overall coefficient matrix, is the fitting coefficient matrix, is the observation error vector, which can be expressed as: ； ； ； In the formula, , is the row number of the ground points at both ends of the terrain contour line, , is the elevation of the terrain contour line at the ground points at both ends. , It is the slope change rate of the ground points at both ends along the DEM elevation strip direction; Then use the least squares estimation to fit the coefficient matrix Solve: ； In the formula, represents the weight matrix, whose value is the unit matrix with diagonal elements of 1, indicating that the weights of the terrain information of the ground points at both ends for the function fitting are equal; According to the fitting coefficient , , Solve for known terrain contour lines and calculate the projection line length of the terrain contour line between adjacent ground points on the O-xz plane by integration: ； In the formula, represents the path of terrain contours, For the i The projection length of the topographic contour line on the O-xz plane; Finally, the projection lengths of the terrain contour lines between all adjacent ground points in the DEM elevation strip are summed to obtain the total projection length 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 the observation, the DEM elevation projection length calculated by the registration model is It can be expressed as a nonlinear function: , then the error function is transformed into: ； Among them, the nonlinear function independent variable It is a certain height measurement point in the ground height measurement strip of the satellite-borne laser radar i The initial row and column number in the geographic coordinate system; are the registration parameters to be solved, including the scaling factor , rotation angle And row and column offsets , ; Initial row and column number The corresponding elevation projection length error; The nonlinear least squares method is used to solve the registration parameters in the error function by minimizing the sum of squared errors to obtain the optimal solution for the registration parameters. : ； In the formula, It is the number of altimeter points contained in the ground altimeter strip of the spaceborne lidar.

6. The method according to claim 1, characterized in that: The specific process of step 6 includes: Step 6.1: For each ground height measurement strip of all spaceborne lidars in the study area, run steps 2 to 5 to obtain the registration parameter vectors of different strips; Step 6.2, construct a polynomial fitting offset to correct the plane position of the entire scene DEM, and its expression is: ； In the formula, and are the polynomial fitting coefficients, ; and is the row and column number of DEM in the geographic coordinate system, and is the offset of DEM in each ground altimetry strip relative to the altimetry point of the satellite-borne lidar; Step 6.3, by the row and column numbers of each altimeter point in the ground altimeter strip of the satellite-borne laser radar 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) , the known fitting coefficients The expression (11) is the offset lookup table of DEM relative to the space-borne lidar data.

7. An inversion device for correcting the plane position of InSAR DEM based on laser radar data, characterized in that: include: The preprocessing module is used to: perform interference 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 of each ground height measurement strip of the satellite-borne laser radar in the area to be studied along the track on the O-xz plane of the spatial rectangular coordinate system, which is recorded as the radar elevation projection length. ; The registration module is used to: construct a registration model to register the DEM and satellite-borne laser radar data of the study area according to the principle of terrain contour consistency, and obtain the plane offset of the DEM relative to the satellite-borne laser radar data in the geographic coordinate system; The registration parameter solution method of the registration model is as follows: obtain the elevation strip corresponding to the DEM according to the geographical location of the radar ground altimetry strip, and then calculate the elevation projection length of the DEM elevation strip along the orbit direction of the satellite-borne laser radar 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 by the nonlinear least square method; The lookup table construction module is used to: take the row and column numbers of all points in the satellite-borne laser radar ground altimetry strip in the geographic coordinate system and the offset of DEM relative to them in the row and column directions as observation quantities, use the least squares adjustment method to solve the offset polynomials of DEM relative to the satellite-borne laser radar ground altimetry strip in the row and column directions respectively, and obtain the offset lookup table of DEM relative to the satellite-borne laser radar data; The DEM correction module is used to: transform the initial DEM coordinates obtained in the preprocessing module to new coordinate positions using the offset lookup table obtained by the lookup table construction module to obtain the DEM after plane position correction.

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

9. A computer-readable storage medium, characterized in that: A computer program is stored thereon, and when the computer program is executed by a processor, the method according to any one of claims 1 to 6 is implemented.

Images