Multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate areas

Abstract

The present invention provides a multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate zone, comprising the following steps: step 1: data collection and preprocessing; step 2: high-resolution multi-orbit dual-station system SAR-DEM generation; step 3: multi-orbit radar DEM fusion generation of high-resolution and high-precision DEM; step 4: multi-source high-resolution optical DEM generation; step 5: DEM differential processing; step 6: target area surface height / volume change monitoring. The present invention develops a surface height / volume monitoring method suitable for complex terrain and climate zone where field observation is difficult to carry out by fully taking into account the respective monitoring advantages of multi-source satellite data (optical, radar, laser), and the present invention will help promote the integrated use of multi-source satellite data and improve the ability of global change monitoring.

Description

Technical Field

[0001] The invention relates to the technical field of using multi-source remote sensing satellite data in monitoring surface changes in complex terrain and climate zones, and in particular to a multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones. Background Art

[0002] Under the background of global warming and increased human activities, the land surface is undergoing complex and profound changes, especially in complex terrain and climate zones characterized by glaciers, frozen soil and snow, where the surface changes are more significant. At the same time, these changes are often closely related to many natural disasters, such as ice avalanches, snow avalanches, rock avalanches, landslides, debris flows, glacial lake outburst floods, etc., which seriously threaten the safety of life and property of residents in nearby areas and the construction and operation safety of infrastructure and projects. In recent years, the frequency of natural disasters developed in complex terrain and climate zones has increased significantly. In addition, according to the current trend of climate change and human activities, the number of related disasters in the future will be more foreseeable, and the impact will be greater, especially in mountainous areas with large populations and high density. Therefore, how to quickly, efficiently and accurately obtain surface height / volume changes is of great significance for disaster assessment, disaster relief guidance and post-disaster reconstruction.

[0003] At present, the detection of changes in surface height / volume mainly relies on ground or near-ground observation methods, such as GPS, leveling, ground-based lidar, drone mapping, etc., on the one hand, and remote sensing satellite observations on the other. The former is mainly suitable for small-scale (such as engineering scale) monitoring, while the latter is used for large-scale surface height changes caused by natural processes, especially for remote mountainous areas that are difficult for people to approach or reach. Regarding satellite altimetry, it is currently divided into three main categories: optical stereo imaging satellites, interferometric SAR satellites (heavy orbit system mode and dual-station system mode), and laser altimetry satellites. Among them, optical stereo imaging satellites are the most numerous and the most widely used. The more representative satellites are Terra-ASTER, ALOS-PRISM, SPOT series, and my country's ZY-3 satellite. In recent years, with the launch of high-resolution and even ultra-high-resolution optical stereo imaging satellites, high-resolution and high-precision digital elevation models can be obtained quickly and on a large scale.

[0004] The capabilities of Digital Elevation Model (DEM) have been significantly improved. For example, the WorldView and GeoEye series satellites of the United States and the Pléiades series satellites of France all have sub-meter spatial resolution. Nevertheless, due to the inherent obstruction of clouds and fog by optical satellites, their application is greatly limited in some areas with complex climate and meteorological conditions. In contrast, synthetic aperture radar (SAR) can achieve all-weather and all-day observations without being disturbed by clouds and fog. The more representative heavy-orbit interferometric SAR satellites currently in orbit include Sentinel-1 (C band) of the European Space Agency and ALOS-2 (L band) of the Japanese Space Agency. For areas with more significant surface changes, heavy-orbit interferometric SAR is still not suitable for complex terrain and climate zones due to the limitation of temporal decoherence. For this reason, the newly developed dual-station interferometric SAR satellites in recent years can overcome the influence of temporal decoherence and realize the reconstruction of digital elevation models of global scale and all terrain elements. They have been successfully applied in the fields of global terrain mapping and global change monitoring. At present, the only satellites in orbit around the world are the German TerraSAR / TanDEM-X satellite system (X-band) and my country's Land Exploration No. 1 satellite system (L-band) to be launched in 2022. Although the dual-station SAR has great advantages, it is still affected by overlap and shadows in rugged mountainous areas with complex terrain due to the inherent limitations of radar sensor imaging geometry, resulting in blind spots or failure areas in monitoring, limiting the ability to monitor surface changes. In addition, for laser altimetry satellites (such as the US ICESat-1 / 2 and my country's Gaofen-7 satellite), data collection is currently mainly achieved in the form of observations of light spot footprints. Therefore, in mid- and low-latitude regions, the observation sampling points are relatively sparse due to orbital restrictions, especially in the cross-track direction. At the same time, laser sensors are also affected by clouds and fog, resulting in less available data in areas with complex climates.

[0005] In summary, affected by climate change and human activities, the surface height / volume changes in complex terrain and climate zones are becoming more frequent, and the impacts caused by them are becoming more significant. However, for current observation techniques, optical satellites are susceptible to cloud and fog obstructions, radar satellites are susceptible to inherent monitoring blind spots such as overlap and shadows, and laser altimetry satellites have the problem of few observation sampling points, which greatly limit and hinder rapid and high-precision monitoring of surface height / volume changes. Therefore, how to comprehensively utilize multi-source, multi-platform, and multi-phase satellite data is a prerequisite for achieving high-precision surface height / volume change estimation in complex terrain and climate zones. However, there is currently no more systematic technical framework to solve this type of problem in a targeted manner. To this end, the present invention provides a technical framework for solving this type of problem. Summary of the invention

[0006] In view of the defects of the prior art, the present invention provides a multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate areas.

[0007] In order to achieve the above invention object, the technical solution adopted by the present invention is as follows:

[0008] The multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate areas includes the following steps:

[0009] Step 1: Data collection and preprocessing: First, collect multi-orbit dual-station SAR data covering the study area, publicly released low-resolution DEM product data, multi-platform multi-temporal high-resolution optical stereo image pairs, laser altimetry satellite data, etc. Then, preprocess the laser altimetry data, mainly including the conversion of coordinate reference system, data screening considering terrain information and removal of abnormal observations;

[0010] Step 2: High-resolution multi-orbit dual-station SAR-DEM generation: First, use the multi-orbit dual-station SAR data to generate interference patterns, then combine with external reference DEM data to obtain differential interference phase maps, then perform data filtering and phase unwrapping, further perform residual orbit error correction on the unwrapped phase map, and finally perform residual height conversion by using a conversion formula that takes into account the changes in imaging geometric parameters, and further combine with the original reference DEM to generate a high-resolution multi-orbit radar DEM of the target area, and finally use laser altimetry data to perform absolute positioning correction on the radar DEM;

[0011] Step 3: Multi-track radar DEM fusion to generate high-resolution and high-precision DEM: First, perform differential processing on radar DEM data from different tracks, check and remove the spatial trend term error and phase unwrapping error between DEMs based on the DEM difference map; then calculate the theoretical height error of each point in the radar DEM, and determine the weight based on it to fuse the DEM data, and finally generate a fused high-quality radar DEM;

[0012] Step 4: Multi-source high-resolution optical DEM generation: First, orthorectify the high-resolution optical image using the obtained high-resolution and high-precision radar DEM. Then, use the radar DEM and the orthorectified image as the vertical reference and horizontal reference to select control points and tie points. When the preliminary estimated aerial triangulation meets the accuracy requirements, perform dense image matching technology and estimate the height value based on the updated exterior orientation elements, and finally generate a high-resolution optical DEM.

[0013] Step 5: DEM differential processing: Using high-quality radar DEM as a reference, the generated multi-source optical DEM is aligned with the radar DEM in sequence using a three-dimensional registration method. Then, the images are arranged in sequence according to the time of image acquisition and differential processing is performed one by one to generate a height change map of the target area over multiple time periods.

[0014] Step 6: Target area surface height / volume change monitoring: Extract the boundaries of potential change areas, then calculate the target area surface height change based on the boundary information, and then estimate the volume of the target area surface change based on the area information.

[0015] Furthermore, the step 1 includes the following sub-steps:

[0016] Step 11: Query and collect archived SAR data of the dual-station system in the study area;

[0017] Step 12: Query and collect publicly released DEM product datasets in the study area;

[0018] Step 13: Query and collect high-resolution optical stereo image pairs from multiple platforms and multiple phases in the study area;

[0019] Step 14: Query and collect laser altimetry satellite data of the study area;

[0020] Step 15: Screen the laser altimetry data and remove outliers.

[0021] Furthermore, the step 2 specifically includes the following sub-steps:

[0022] Step 21: Generate radar interferograms one by one using multi-orbit dual-station SAR images;

[0023] Step 22: Simulate terrain phase using low-resolution external DEM;

[0024] Step 23: Subtract the simulated terrain phase from the interference pattern to obtain a differential interference phase map, simulate the overlap and shadow areas of the radar image based on the external DEM data, and perform mask processing on the differential interference pattern accordingly;

[0025] Step 24: Adaptively filter the differential interference pattern after mask processing;

[0026] Step 25: Phase unwrapping of the filtered differential interferogram using the minimum cost flow method;

[0027] Step 26: Correct the residual orbit error using an empirical polynomial fitting method;

[0028] Step 27: Then, the corrected unwrapped phase is transformed to generate residual height, and it is added to the reference DEM to obtain a high-resolution multi-track radar DEM;

[0029] Step 28: Use laser altimetry data to perform absolute positioning correction on the generated DEM data.

[0030] Furthermore, the step 27 is specifically as follows: in the process of phase-to-height conversion, in order to take into account the spatial variation of imaging parameters (such as incident angle, spatial baseline and slant range, etc.) under the SAR field of view, an improved conversion formula is adopted, which can be expressed as formula (1):

[0031]

[0032] Where Δh represents the residual height difference, λ represents the radar wavelength, R represents the slant range, θ represents the incident angle, and B perp Represents the vertical baseline, represents the differential interference phase, and i represents the column number of data sampling along the distance direction.

[0033] Furthermore, the step 3 specifically includes the following sub-steps:

[0034] Step 31: First, perform data registration and differential processing on the DEMs of different tracks;

[0035] Step 32: Check and correct the spatial trend term error between DEMs;

[0036] Step 33: removing the unwrapping error in the local area;

[0037] Step 34: Calculate the theoretical height error of each single point in the DEM;

[0038] Step 35: Use weighted averaging to fuse the DEMs of different tracks and generate high-resolution and high-precision radar DEM.

[0039] Furthermore, the step 34 is specifically as follows: for the DEM generated by the dual-station SAR, the theoretical height error of a single point is related to the interferometric phase standard deviation and the height ambiguity, which can be expressed as formula (2);

[0040]

[0041] Among them, σ h Expressed as the height error of a single point, represents the standard deviation of the interferometric phase, which can be further calculated by formula (3). HoA represents the height ambiguity, which can be calculated by formula (4) for the dual-station SAR mode.

[0042]

[0043] Furthermore, the step 35 is specifically as follows: after obtaining the theoretical height error corresponding to each point in the DEM data, the weight factor (ω k ) can be expressed as formula (5). The fused elevation (h fusion ) can be further calculated by formula (6).

[0044]

[0045] Among them, k represents the number of radar DEMs in the same area. m represents the shadow and overlap binary image, 0 represents the overlap or shadow area, and 1 represents the effective pixel area.

[0046] Furthermore, the step 4 specifically includes the following sub-steps:

[0047] Step 41: Orthorectify the multi-source high-resolution images based on the generated high-resolution and high-precision DEM, and output them at the original resolution of the optical images;

[0048] Step 42: Using the orthorectified optical image as a horizontal reference and the high-resolution radar DEM as a vertical reference, the image matching technology is used to automatically generate control points and tie points, and the control points and tie points with low accuracy are checked and eliminated according to the aerial triangulation calculation results;

[0049] Step 43: Generate a multi-source high-resolution optical DEM based on photogrammetry using dense image matching techniques and updated exterior orientation elements.

[0050] Furthermore, the step 5 specifically includes the following sub-steps:

[0051] Step 51: Using the high-quality radar DEM as a reference, align all generated DEMs with it in turn using a three-dimensional analytical registration method, and check and correct the spatial trend item errors in turn;

[0052] Step 52: Arrange all the registered DEM data in sequence according to the image acquisition time and perform differential processing to obtain the surface height change map of the study area in multiple time periods.

[0053] Furthermore, the step 6 specifically includes the following sub-steps:

[0054] Step 61: According to the surface height change map, determine the region where the surface height changes, and calculate the area of ​​the change region;

[0055] Step 62: Calculate the volume change of the target area surface according to the regional area and height change.

[0056] Compared with the prior art, the advantages of the present invention are:

[0057] 1. The multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones of the present invention fully takes into account the respective monitoring advantages of multi-source satellite data (optical, radar, laser), realizes the complementary advantages and data integration of various types of satellites, and can realize high-precision surface height / volume change estimation in complex terrain and climate zones, thereby enhancing the comprehensive collaborative application capability of multi-source satellite data and the ability to monitor surface height / volume changes.

[0058] 2. The present invention effectively solves the current problem of difficult on-site monitoring of surface height / volume changes in complex terrain and climate zones, and low or even ineffective single-source satellite observation accuracy. At the same time, the method and technical solution of the present invention is reasonable, highly feasible, and can be extended to any field where surface height / volume changes occur, thus helping to promote the integrated use of multi-source satellite data and improve the ability to monitor global changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a schematic diagram of a flow chart of an embodiment of the present invention;

[0060] Figure 2(a) shows the initial radar interferogram;

[0061] Figure 2(b) is the differential interferogram with the terrain phase removed;

[0062] Figure 2(c) is the estimated coherence graph;

[0063] Figure 2(d) is the differential interferogram after phase unwrapping;

[0064] Figure 2(e) is a phase image with the orbital trend removed and the overlap and shadow masking processed;

[0065] Figure 3(a) shows the digital elevation model obtained from the ascending SAR data;

[0066] Figure 3(b) shows the digital elevation model obtained from the descending SAR data;

[0067] Figure 3(c) is a high-resolution and high-precision digital elevation model obtained after data fusion;

[0068] Figure 4 The high-resolution optical DEM obtained by using ultra-high-resolution Pleiades optical stereo image pairs in the southwest of my country in the embodiment of the present invention;

[0069] Figure 5(a) to Figure 5(b) This is a height change map obtained by differential processing of high-resolution radar DEM and high-resolution optical DEM in the southwest of my country in the embodiment of the present invention;

[0070] Figure 5(a) shows the altitude variation of the glacier area and its surrounding areas;

[0071] Figure 5(b) shows the relationship between the change in glacier area height and altitude. DETAILED DESCRIPTION

[0072] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0073] like Figure 1 As shown, the multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate areas includes the following steps:

[0074] Step 1: Data collection and preprocessing. First, collect the dual-station SAR data, publicly released external DEM product data, laser altimeter satellite data, and multi-source high-resolution optical stereo image pair data required for the study area; then perform data screening and detection and removal of abnormal observations for the laser altimeter data considering terrain complexity;

[0075] Step 2: Generate multi-track dual-station SAR-DEM data. First, perform interferometric processing on the dual-station SAR image to generate a radar interferogram. Then, use the low-resolution external DEM as a reference to generate a simulated terrain phase. Further, perform differential processing to obtain a differential interferogram. Then, perform adaptive filtering on it. After that, use the minimum cost flow method to perform phase unwrapping and mask the overlapped and shadowed areas. Further, use the polynomial fitting method to correct the residual orbit error. The corrected unwrapped phase is converted to the height difference relative to the reference DEM using formula (1). Finally, it is added to the reference DEM to obtain a high-resolution radar DEM, and then laser altimetry data is used to perform absolute positioning correction on it.

[0076]

[0077] Where Δh represents the residual height difference, λ represents the radar wavelength, R represents the slant range, θ represents the incident angle, and B perp Represents the vertical baseline, represents the differential interference phase, and i represents the column number of data sampling along the distance direction.

[0078] Step 3: Multi-track radar DEM data fusion. The multi-track radar DEM obtained in step (2) is first subjected to DEM registration and differential processing. The spatial trend error between DEMs is checked using the elevation difference map of the stable area and corrected using polynomial fitting. The unwrapping error is then checked and removed. The theoretical height error of each DEM is estimated using formulas (2), (3) and (4), and the DEM data is fused using weighted average using this as a weighted indicator. Finally, a high-resolution and high-precision radar DEM is generated. This process can be specifically expressed as formulas (5) and (6).

[0079]

[0080] Among them, σ h Expressed as the height error of a single point, represents the standard deviation of the interferometric phase, which can be further calculated by formula (3). HoA represents the height ambiguity, which can be calculated by formula (4) for the dual-station SAR mode.

[0081]

[0082] Among them, k represents the number of radar DEMs in the same area. m represents the shadow and overlap binary image, 0 represents the overlap or shadow area, and 1 represents the effective pixel area.

[0083] Step 4: Generate multi-source high-resolution optical DEM. First, use the high-resolution and high-precision radar DEM obtained after data fusion in step (3) to orthorectify the high-resolution optical image. Then, use the radar DEM and the orthorectified image as vertical reference and horizontal reference to select control points and connection points. When the preliminary estimated aerial triangulation meets the accuracy requirements, perform dense image matching technology and estimate the height value based on the exterior orientation elements, and finally generate a high-resolution optical DEM.

[0084] Step 5: DEM differential processing. First, the multi-source optical DEM obtained in step (4) is aligned with the radar DEM using the three-dimensional analytical registration method, and the spatial trend item errors are checked and corrected in turn. Then, the data are arranged according to the time of data acquisition and differential processing is performed pair by pair to generate a height change map of the target area for multiple time periods.

[0085] Step 6: Monitoring the surface height / volume changes in the target area: Based on the height change map obtained in step (5), the boundary of the potential change area is extracted, the area of ​​the change area and the average height change are calculated, and then the volume of the surface change in the target area is estimated.

[0086] In order to effectively solve the key problem in the current field of surface change monitoring - the estimation of height / volume changes in complex terrain and climate zones, the present invention provides a multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones by integrating multi-source satellite data.

[0087] The multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones described in the above-mentioned embodiments of the present invention can quickly and accurately realize large-scale surface height / volume change estimation by fully considering the respective monitoring advantages of multi-source satellite data (optical, radar, laser). The technical system has strong operability and scalability, and therefore will help to popularize and promote the application of multi-source satellite data in the field of surface change monitoring, and will also help to provide scientific guidance for disaster assessment and post-disaster emergency response.

[0088] When implementing the present invention, the southwestern region of my country was selected as the experimental test area. Affected by the southwest monsoon of the Indian Ocean, the region has abundant precipitation and is covered by clouds all year round. At the same time, the terrain is complex and the mountain range has a large drop, forming a natural landscape with vertical distribution of glaciers, snow and forests, which is a typical complex terrain and climate zone. This experiment mainly uses high-precision height / volume measurement in glacier areas to illustrate the effectiveness of the present invention, wherein the specific implementation steps of the present invention are as follows:

[0089] Step 1: TerraSAR / TanDEM dual-station SAR data (spatial resolution 2 meters) of the test area were collected in ascending and descending orbits. The external DEM used the latest NASA DEM (after data correction and interpolation, with a spatial resolution of 30 meters). The laser altimetry data used the ICESat-1 / 2 data released by the US Snow and Ice Science Data Center, and the high-resolution optical stereo images used the Pléiades data (spatial resolution of 0.5 meters). Then, the slope threshold was used to filter the ICESat-1 / 2 data, retaining the footprint point data in the flat area, and further remove the outliers based on the original data record file.

[0090] Step 2: Generate interferograms using the ascending and descending TerraSAR / TanDEM dual-station SAR data, as shown in Figure 2(a). Set the multi-view factor to 4*5 (range and azimuth) so that the output resolution is 10m. Then, take NASADEM as a reference and simulate the terrain phase with a resolution of 10m by upsampling. Then, perform differential processing to obtain the differential interferogram, and perform adaptive filtering on the corrected differential interferogram, as shown in Figure 2(b). Then, use the minimum cost flow method to phase unwrap the differential interferogram, as shown in Figure 2(d). Further, use polynomial fitting to remove the residual orbit error, as shown in Figure 2(e). Finally, use formula (1) to convert the unwrapped differential interferogram phase to the corresponding residual height, and add it to the upsampled reference DEM to obtain high-resolution ascending and descending radar DEMs. Then, use laser altimetry data to correct and adjust the absolute positioning of the DEM. The final results are shown in Figures 3(a) and 3(b).

[0091] Step 3: The ascending and descending DEMs obtained in step (2) are firstly registered and differentially processed, and the spatial trend term error is checked and corrected using the polynomial fitting method. Then, the phase unwrapping error is removed for pixels greater than 3 times the standard deviation. Then, the theoretical single-point error of the ascending and descending DEM data is calculated using formulas (2), (3) and (4), respectively. Then, the DEM data is fused using formulas (5) and (6). The fused high-resolution and high-precision DEM is shown in Figure 3(c).

[0092] Step 4: Use the high-quality radar DEM obtained in step (3) to orthorectify the Pléiades high-resolution optical image (panchromatic band data) and output it at the original resolution; then use the orthorectified optical image and radar DEM as the horizontal reference and vertical reference respectively, and automatically generate control points and connection points for the Pléiades stereo image pair through image matching technology. After manual inspection, only 12 control points and 24 connection points are retained, and the residual error of aerial triangulation is 0.5 pixel (less than the accuracy requirement of 1 pixel); then, based on the obtained imaging geometric parameters and using dense matching technology, the final optical DEM is generated, and the spatial resolution is adjusted to 10m through downsampling, as shown in Figure (4).

[0093] Step 5: The two DEMs obtained in steps (3) and (4) are registered with a three-dimensional analytical registration method based on terrain slope and aspect information (with a registration accuracy of 0.1 pixel). Then, the residual spatial trend term error is checked and corrected using a polynomial fitting method to obtain the final elevation difference map of the target area, as shown in Figure 5(a).

[0094] Step 6: According to step (5), the elevation difference map of the target area is obtained and the boundary of the target change area is delineated in combination with the optical image map after the surface change (in this experimental case, the glacier boundary data set is directly used) to obtain the area information of the change area, and further combined with the height information to obtain the volume change result.

[0095] Figure 5(a) to Figure 5(b)It is the result of the surface height change obtained by the embodiment of the present invention in the test area in southwest my country. From the height change diagram (the value range is ±15 meters, the black line is the glacier boundary, Figure 5 (a)), there is no residual spatial trend error in the stable area (non-glacier area), and the statistical histogram shows that it basically obeys the normal distribution, indicating that the height error is mainly random error and there is no systematic deviation. This also shows that the height change of the target area (glacier area) is approximately unbiased. From the height change of the glacier area, there is obviously a trend change along the rising altitude, as shown in Figure 5 (b). In the altitude band of 4900-5200m, the height change is almost 0; in the altitude band of 5200-5400m, the height change shows a linear increase trend, and then gradually stabilizes, with an average change of about 4 meters. It should be noted that the acquisition time of the radar DEM and optical DEM used in this embodiment is only 4 days different, and the final generated surface height change actually represents the penetration depth of the X-band in the ice and snow covered area. Judging from the final results, the method of the present invention is feasible and effective, and can fully realize high-precision surface height / volume change monitoring in complex terrain and climate zones.

[0096] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the implementation methods of the present invention, and should be understood that the protection scope of the present invention is not limited to such special statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the present invention.

Claims

1. A multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones, characterized in that: The following steps are involved: Step 1: Data collection and preprocessing: First, collect multi-orbit dual-station SAR data covering the study area, publicly released low-resolution DEM product data, multi-platform multi-temporal high-resolution optical stereo image pairs, and laser altimetry satellite data; then preprocess the laser altimetry data, including the conversion of the coordinate reference system, data screening considering terrain information, and removal of abnormal observations; Step 2: Generate high-resolution multi-orbit dual-station SAR-DEM; It includes the following sub-steps: Step 21: Generate radar interferograms one by one using multi-orbit dual-station SAR images; Step 22: Simulate terrain phase using low-resolution external DEM; Step 23: Subtract the simulated terrain phase from the interference pattern to obtain a differential interference phase map, simulate the overlap and shadow areas of the radar image based on the external DEM data, and perform mask processing on the differential interference pattern accordingly; Step 24: Adaptively filter the differential interference pattern after mask processing; Step 25: Phase unwrapping of the filtered differential interferogram using the minimum cost flow method; Step 26: Correct the residual orbit error using an empirical polynomial fitting method; Step 27: Then, the corrected unwrapped phase is transformed to generate residual height, and it is added to the reference DEM to obtain a high-resolution multi-track radar DEM; Step 28: Use laser altimetry data to perform absolute positioning correction on the generated DEM data; Step 3: Multi-track radar DEM fusion generates high-resolution and high-precision DEM; Step 4: Generate multi-source high-resolution optical DEM; Step 5: DEM difference processing; Step 6: Monitoring of surface height / volume changes in the target area.

2. The multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones according to claim 1, characterized in that: The step 27 is specifically as follows: in the process of phase-to-height conversion, in order to take into account the spatial variation of imaging parameters under the SAR field of view, an improved conversion formula is adopted, which is expressed as formula (1): Where Δh represents the residual height difference, λ represents the radar wavelength, R represents the slant range, θ represents the incident angle, and B perp Represents the vertical baseline, represents the differential interference phase, and i represents the column number of data sampling along the distance direction.

3. The multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate zones according to claim 1, characterized in that: The step 3 includes the following sub-steps: Step 31: First, perform data registration and differential processing on the DEMs of different tracks; Step 32: Check and correct the spatial trend term error between DEMs; Step 33: removing the unwrapping error in the local area; Step 34: Calculate the theoretical height error of each single point in the DEM; Step 35: Use weighted averaging to fuse the DEMs of different tracks and generate high-resolution and high-precision radar DEM.

4. The multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones according to claim 3 is characterized by: The step 34 is specifically as follows: for the DEM generated by the dual-station SAR, the theoretical height error of a single point is related to the interferometric phase standard deviation and the height ambiguity, which is expressed as formula (2); Among them, σ h Expressed as the height error of a single point, represents the standard deviation of the interferometric phase, which is further calculated by formula (3). HoA represents the height ambiguity, which is calculated by formula (4) for the dual-station SAR mode: The step 35 is specifically as follows: after obtaining the theoretical height error corresponding to each point in the DEM data, the weight factor ω of the data fusion is k Expressed as formula (5); the fused elevation h fusion Further calculation by formula (6) yields: Among them, k represents the number of radar DEMs in the same area; m represents the shadow and overlap binary image, 0 represents the area is an overlap or shadow area, and 1 represents the area is a valid pixel area.

5. The multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate zones according to claim 1, characterized in that: The step 4 includes the following sub-steps: Step 41: Orthorectify the multi-source high-resolution images based on the generated high-resolution and high-precision DEM, and output them at the original resolution of the optical images; Step 42: Using the orthorectified optical image as a horizontal reference and the high-resolution radar DEM as a vertical reference, the image matching technology is used to automatically generate control points and tie points, and the control points and tie points with low accuracy are checked and eliminated according to the aerial triangulation calculation results; Step 43: Generate a multi-source high-resolution optical DEM based on photogrammetry using dense image matching techniques and updated exterior orientation elements.

6. The multi-source remote sensing monitoring method for surface height or volume change in complex terrain and climate zones according to claim 1, characterized in that: The step 5 comprises the following sub-steps: Step 51: Using the high-quality radar DEM as a reference, align all generated DEMs with it in turn using a three-dimensional analytical registration method, and check and correct the spatial trend item errors in turn; Step 52: Arrange all the registered DEM data in sequence according to the image acquisition time and perform differential processing to obtain the surface height change map of the study area in multiple time periods.

7. The multi-source remote sensing monitoring method for surface height or volume changes in complex terrain and climate zones according to claim 1, characterized in that: The step 6 comprises the following sub-steps: Step 61: According to the surface height change map, determine the region where the surface height changes, and calculate the area of ​​the change region; Step 62: Calculate the volume change of the target area surface according to the regional area and height change.

Images