A terrain modeling method based on collaborative multi-source satellite remote sensing data

Abstract

The invention discloses a terrain modeling method of multi-source satellite remote sensing data collaboration, which relates to the fields of satellite photogrammetry and geographic information technology. Firstly, multi-attribute criteria screening is performed using satellite-borne laser altimetry ICESat-2 data to obtain some high-precision data as high-quality control points, which are then introduced into the adjustment equation of a satellite stereo image pair to generate a dense point cloud, and then a large-area, high-precision DEM data is obtained through filtering operations. In some areas where DEM cannot be correctly produced due to cloud cover, a deep learning method and open-source DEM data are used to reconstruct the terrain. Finally, the DEM of a cloud-free area is fused with the DEM reconstructed by deep learning in a cloud area to obtain a large-scale, high-precision, full-coverage DEM, thereby concentrating the advantages of several multi-source satellite remote sensing data and providing relatively accurate data products for earth science surface modeling and analysis.

Description

Technical Field

[0001] The invention relates to the field of satellite photogrammetry and geographic information technology, and in particular to a terrain modeling method using multi-source satellite remote sensing data in collaboration. Background Art

[0002] Terrain modeling, as a hot topic in the field of geology, has important applications in people's livelihood and national defense. From the perspective of geology and geomorphology, terrain modeling is an accurate description of the shape of the earth's surface and the result of various landform evolution processes. It not only reflects the role of violent crustal movements such as geological changes, but also reflects minor changes such as soil erosion and soil loss. Therefore, terrain modeling is widely involved in scientific research fields such as environmental change, disaster prevention and control, hydrological analysis, and geological analysis, as well as engineering application scenarios such as military battlefields and game modeling, and is of great significance to engineering construction and national defense.

[0003] Satellite photogrammetry covers a wide area and can achieve sub-meter accuracy, but it is easy to cause obstruction in some high vegetation areas and areas with more clouds and other aerosols, thus affecting data quality. To meet the needs of large-scale, high-precision, full coverage, and low-cost terrain modeling, satellite mapping can only be used to produce regional terrain modeling. For cloud coverage and other aerosol obstruction problems, they are generally solved by adding other stereo pairs. However, since the cloud content of satellite photogrammetry stereo pairs is mostly between 1% and 10%, adding other stereo pairs to solve the problem will double the cost and data processing time, and the cost-effectiveness is low; the characteristic of deep learning technology is a large amount of prior data, so it is just right for cloud and other aerosol obstruction problems, and has important application prospects. Summary of the invention

[0004] In order to solve the above technical problems, the present invention provides a terrain modeling method based on multi-source satellite remote sensing data, comprising the following steps:

[0005] S1. Perform multi-attribute constraint screening on ICESat-2 satellite laser altimetry data in the terrain modeling sample area to obtain high-quality points after screening;

[0006] S2, relative orientation of the photogrammetric satellite stereo image pair of the area to be measured is performed by solving the rational polynomial function model, and then absolute orientation is performed using the screened high-quality points to generate a dense point cloud and DEM, and then the dense point cloud is filtered to generate a DEM;

[0007] S3, manually visually distinguish and extract the cloud-blocked area from the photogrammetric satellite stereo image pair of the area to be measured, and then perform mask processing on the DEM generated in step S2 to cut out the cloud-blocked area;

[0008] S4. Use hydrological analysis methods to extract valley lines, ridge lines, and gully lines from the DEM with the cloud and fog areas cut out;

[0009] S5, merging the valley lines, ridge lines and ditch lines extracted in step S4 into a three-band tif file, and cutting them according to the pixel size of 256×256 to produce a training data set, and then inputting the training data set into the conditional generative adversarial network CGAN for training to obtain a CGAN network with a mapping relationship between DEM and terrain skeleton;

[0010] S6, using the method in step S4 to extract the ridge lines and valley lines of the cloud-covered area from the open source 30m resolution ASTER GDEM, and resample them to 5m;

[0011] S7, using the method in step S4 to extract the gully line from the open source 1m resolution Google remote sensing image data and resample it to 5m, then merge the ridge line, valley line, and gully line extracted from the open source data into a three-band tif file and crop it according to 256×256 to form an open source terrain feature dataset;

[0012] S8, input the open source terrain feature dataset into the conditional generative adversarial network CGAN trained in step S5, and output the DEM reconstructed by deep learning in the cloud and fog area;

[0013] S9, merging the cloud-free area Gaofen-7 DEM generated in step S3 with the cloud-fog area DEM generated in step S8.

[0014] The technical solution further defined in the present invention is:

[0015] Furthermore, in step S1, the screening criteria for multi-attribute constraint screening of ICESat-2 satellite-borne laser altimetry data in the terrain modeling sample area are:

[0016] i. The cloud cover mark is less than 3;

[0017] ii. The city flag is equal to 0;

[0018] iii. The quantity mark is greater than 50;

[0019] iv. The difference between the laser point elevation and the ASTER GDEM elevation value is less than 75;

[0020] v. The difference between the average elevation of the laser point in each 100-meter statistical unit and the elevation value of the laser point interpolation surface is less than 75.

[0021] As described above, in the method for terrain modeling based on multi-source satellite remote sensing data collaboration, in step S2, the rational polynomial function model is as follows:

[0022]

[0023] in, represents the coordinates of the image point on the image plane, Indicates the image point at the corresponding ground coordinate, Represents the rational polynomial function model to be solved.

[0024] As described above, in a terrain modeling method using multi-source satellite remote sensing data, in step S2, absolute orientation is performed using the high-quality points obtained by screening as shown in the following formula:

[0025]

[0026] in, Indicates the 7 absolute orientation parameters that need to be solved; Represents the ground coordinate value of the model point; Represents the coordinates of the same model point in the auxiliary coordinate system of the image space; Represents the ground coordinate value of the model point; Indicates the model scaling factor; Indicates the photo posture angle The direction cosines are calculated as follows:

[0027]

[0028] in, Indicates the photo pose angle.

[0029] As described above, in a terrain modeling method based on multi-source satellite remote sensing data collaboration, in step S4, the valley line extraction method is specifically as follows: the runoff accumulation of each pixel is calculated based on the DEM with the cloud area cut out, and then the runoff accumulation threshold is used for segmentation, and pixels greater than the threshold are considered to be valley line pixels.

[0030] As described above, in a terrain modeling method based on collaborative multi-source satellite remote sensing data, in step S4, the ridge line extraction method is specifically as follows: the DEM with the cloud area cut out is inverted, and then the runoff accumulation of each pixel is calculated based on the inverted DEM, and then the runoff accumulation threshold is used for segmentation, and pixels greater than the threshold are considered to be ridge line pixels.

[0031] As described above, in a terrain modeling method using multi-source satellite remote sensing data, in step S4, the method for extracting the ditch line is specifically as follows: a 30*30 window analysis is performed on the DEM where the cloud area is cut out, and a statistical average is obtained to obtain an average DEM; then the DEM where the cloud area is cut out is subtracted from the average DEM to extract areas where the value is less than 0; then the boundary line of the area where the value is less than 0 is extracted to obtain the ditch line.

[0032] As described above, in a terrain modeling method for collaborative multi-source satellite remote sensing data, in steps S5 and S8, the conditional generative adversarial network CGAN includes a generator network and a discriminant network. The generator network is used to generate synthetic data from random noise and conditional information, and its network structure adopts a U-Net structure; the discriminant network is used to distinguish between data generated by the generator and real data, and its network structure only includes encoding operations, does not include decoding operations, and the activation function adopts tanh.

[0033] As described above in the method for terrain modeling based on collaborative multi-source satellite remote sensing data, in steps S6 and S7, the resampling method is set to the nearest neighbor sampling method.

[0034] As described above, a terrain modeling method based on the collaboration of multi-source satellite remote sensing data, in step S9, the specific method for fusing the GF-7 DEM in the cloud-free area with the DEM in the foggy area is as follows: a buffer zone with a width 5 times the resolution is generated at the junction of the GF-7 DEM in the cloud-free area and the DEM in the foggy area, the elevation value in the buffer zone uses the average value of the GF-7 DEM in the cloud-free area and the DEM in the foggy area, the elevation value of the GF-7 DEM is used in the cloud-free area outside the buffer zone, and the elevation value reconstructed by the conditional generative adversarial network (CGAN) is used in the foggy area outside the buffer zone.

[0035] The beneficial effects of the present invention are:

[0036] In the present invention, firstly, ICESat-2 satellite laser altimetry data is used to perform multi-attribute criterion screening to obtain some high-precision data as high-quality control points, and then the data is brought into the adjustment equation of the satellite stereo image pair to generate a dense point cloud, and then a large-area high-precision DEM data is obtained through filtering operation; in some areas where DEM cannot be correctly produced due to cloud cover, a deep learning method and open source DEM data are used to reconstruct the terrain; finally, the DEM of the cloud-free area is fused with the DEM reconstructed by deep learning in the cloud area to obtain a large-scale, high-precision, full-coverage DEM; when applied to areas with a single data shortage, the use of multi-source data can greatly improve the accuracy of terrain modeling, and the method of the present invention can also use low-resolution data to generate large-scale high-resolution terrain data, the generated model has excellent stability and flexibility, and the local detail adjustment of the generated terrain can be achieved by modifying the terrain features input by the model, thereby concentrating the advantages of several multi-source satellite remote sensing data, and providing more accurate data products for earth science surface modeling and analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a schematic diagram of the overall process of the present invention;

[0038] Figure 2 This is a comparison diagram before and after multi-attribute constraint screening of ICESat-2 satellite-borne laser altimetry data of a sample area in an embodiment of the present invention;

[0039] Figure 3 Schematic diagram of the Gaofen-7 satellite stereo image pair image of the sample area in the embodiment of the present invention;

[0040] Figure 4 is a comparison diagram before and after filtering in an embodiment of the present invention;

[0041] Figure 5 A schematic diagram of the process of generating DEM by deep learning in an embodiment of the present invention;

[0042] Figure 6 This is a schematic diagram of the fusion process of the cloud and fog area DEM generated by deep learning and the GF-7 DEM in the cloudless area in an embodiment of the present invention. DETAILED DESCRIPTION

[0043] This embodiment provides a terrain modeling method using multi-source satellite remote sensing data. Figure 1 As shown, the following steps are included:

[0044] S1. Figure 2 The ICESat-2 satellite laser altimetry data of the terrain modeling sample area shown in the figure on the left is subjected to multi-attribute constraint screening to obtain the following Figure 2 The figure on the right shows some high-quality points with relatively accurate numerical values ​​after screening. 309 points are retained after screening of the original 14,957 points. The screening criteria for multi-attribute constraint screening are:

[0045] i. The cloud cover mark is less than 3;

[0046] ii. The city flag is equal to 0;

[0047] iii. The quantity mark is greater than 50;

[0048] iv. The difference between the laser point elevation and the ASTER GDEM elevation value is less than 75;

[0049] v. The difference between the average elevation of the laser point in each 100-meter statistical unit and the elevation value of the laser point interpolation surface is less than 75.

[0050] S2, by solving the rational polynomial function model (Rational Function Model, RFM) Figure 3 The high-resolution 7 stereo image pair shown in the figure is relatively oriented, and then the high-quality points obtained by screening are used for absolute orientation to generate the following Figure 4 The dense point cloud and DEM shown in a and c are then filtered using cloth to retain only the ground points to generate a digital elevation model (DEM). Figure 4 As shown in b and d.

[0051] The rational polynomial function model is shown below:

[0052]

[0053] in, represents the coordinates of the image point on the image plane, Indicates the image point at the corresponding ground coordinate, Represents the rational polynomial function model to be solved.

[0054] The absolute orientation model is:

[0055]

[0056] in, Indicates the 7 absolute orientation parameters that need to be solved; Represents the ground coordinate value of the model point; Represents the coordinates of the same model point in the auxiliary coordinate system of the image space; Represents the ground coordinate value of the model point; Represents the model scaling factor.

[0057] Absolute directional model formula Indicates the photo posture angle The direction cosines are calculated as follows:

[0058]

[0059] in, Indicates the photo pose angle.

[0060] S3, manually visually distinguish and extract the cloud-blocked area from the photogrammetric satellite stereo image pair of the area to be measured, and then perform mask processing on the DEM generated in step S2 to cut out the cloud-blocked area.

[0061] S4. Use hydrological analysis methods to extract valley lines, ridge lines and gully lines from the DEM that cuts out the cloud and fog areas.

[0062] The specific method for extracting valley lines is as follows: the runoff accumulation of each pixel is calculated based on the DEM with the cloud area cut out, and then the runoff accumulation threshold is used for segmentation. The runoff accumulation threshold is selected in the range of 1000-2000, and pixels greater than the threshold are considered to be valley line pixels.

[0063] The specific method for extracting the ridge line is to invert the DEM with the cloud area cut off, then calculate the convergence accumulation of each pixel based on the inverted DEM, and then use the convergence accumulation threshold for segmentation, and consider the pixels greater than the threshold as ridge line pixels.

[0064] The specific method for extracting the ditch line is as follows: perform a 30*30 window analysis on the DEM where the cloud area is cut off and calculate the statistical average to obtain the average DEM; then use the original DEM to subtract the average DEM to extract the area with a value less than 0; then extract the boundary line of the area with a value less than 0 to obtain the ditch line.

[0065] S5. The valley lines, ridge lines and ditch lines extracted in step S4 are merged into a three-band tif file, and are cropped according to a pixel size of 256×256 to produce a training data set. The training data set is then input into a conditional generative adversarial network (CGAN) for training to obtain a CGAN network with a mapping relationship between the DEM and the terrain skeleton.

[0066] The conditional generative adversarial network CGAN consists of two parts: a generator network and a discriminator network. The generator network is used to generate synthetic data from random noise and conditional information, and its network structure adopts the U-Net structure; the discriminator network is used to distinguish between the data generated by the generator and the real data. Its network structure only includes encoding operations, does not include decoding operations, and the activation function uses tanh.

[0067] S6. Use the method in step S4 to extract the ridge lines and valley lines of the cloud-covered area from the open source 30m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM), and resample them to 5m using the nearest neighbor sampling method.

[0068] S7. Use the method in step S4 to extract the gully line from the open source 1m resolution Google remote sensing image data and resample it to 5m using the nearest neighbor sampling method. Then merge the ridge line, valley line, and gully line extracted from the open source data into a three-band tif file and crop it according to 256×256 to form an open source terrain feature dataset.

[0069] S8. Input the open source terrain feature dataset into the conditional generative adversarial network CGAN trained in step S5, and the output is as follows Figure 5 The DEM reconstructed by deep learning in the cloud area shown.

[0070] S9, the cloud-free area Gaofen-7 DEM generated in step S3 is merged with the cloud-fog area DEM generated in step S8, as shown in FIG. Figure 6As shown in the figure, the specific fusion method is: generate a buffer with a resolution 5 times wider at the junction of the GF-7 DEM in the cloudless area and the DEM in the foggy area, use the average value of the GF-7 DEM in the cloudless area and the DEM in the foggy area as the elevation value in the buffer, use the elevation value of the GF-7 DEM in the cloudless area outside the buffer, and use the elevation value reconstructed by the conditional generative adversarial network (CGAN) in the foggy area outside the buffer.

[0071] The method of this embodiment first uses the satellite-borne laser altimetry ICESat-2 data to perform multi-attribute criterion screening to obtain some high-precision data as high-quality control points, and then brings them into the adjustment equation of the satellite stereo image pair to generate a dense point cloud, and then obtains large-area, high-precision DEM data through filtering operations; in some areas where DEM cannot be correctly produced due to cloud cover, deep learning methods and open source DEM data are used to reconstruct the terrain; finally, the DEM of the cloud-free area is fused with the deep learning reconstructed DEM of the cloud area to obtain a large-scale, high-precision, full-coverage DEM; when applied to areas with a single data shortage, the use of multi-source data can greatly improve the accuracy of terrain modeling, and the method of the present invention can also use low-resolution data to generate a large range of high-resolution terrain data, the generated model has excellent stability and flexibility, and the local detail adjustment of the generated terrain can be achieved by modifying the terrain features input by the model, thereby concentrating the advantages of several multi-source satellite remote sensing data, and can provide more accurate data products for earth science surface modeling and analysis.

[0072] In addition to the above embodiments, the present invention may also have other implementation modes. Any technical solution formed by equivalent replacement or equivalent transformation falls within the protection scope required by the present invention.

Claims

1. A terrain modeling method based on multi-source satellite remote sensing data, characterized in that: The following steps are involved: S1. Perform multi-attribute constraint screening on ICESat-2 satellite laser altimetry data in the terrain modeling sample area to obtain high-quality points after screening; S2, relative orientation of the photogrammetric satellite stereo image pair of the area to be measured is performed by solving the rational polynomial function model, and then absolute orientation is performed using the screened high-quality points to generate a dense point cloud and DEM, and then the dense point cloud is filtered to generate a DEM; S3, manually visually distinguish and extract the cloud-blocked area from the photogrammetric satellite stereo image pair of the area to be measured, and then perform mask processing on the DEM generated in step S2 to cut out the cloud-blocked area; S4. Use hydrological analysis methods to extract valley lines, ridge lines, and gully lines from the DEM with the cloud and fog areas cut out; S5, merging the valley lines, ridge lines and ditch lines extracted in step S4 into a three-band tif file, and cutting them according to the pixel size of 256×256 to produce a training data set, and then inputting the training data set into the conditional generative adversarial network CGAN for training to obtain a CGAN network with a mapping relationship between DEM and terrain skeleton; S6, extracting the ridge lines and valley lines of the cloud-covered area from the open-source 30m resolution ASTER GDEM using the method described in step S4, and resampling them to 5m; S7, using the method described in step S4 to extract the gully line from the open source 1m resolution Google remote sensing image data and resample it to 5m, then merge the ridge line, valley line, and gully line extracted from the open source data into a three-band tif file and crop it according to 256×256 to form an open source terrain feature dataset; S8, input the open source terrain feature dataset into the conditional generative adversarial network CGAN trained in step S5, and output the DEM reconstructed by deep learning in the cloud and fog area; S9, merging the cloud-free area Gaofen-7 DEM generated in step S3 with the cloud-fog area DEM generated in step S8.

2. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S1, the screening criteria for multi-attribute constraint screening of ICESat-2 satellite-borne laser altimetry data in the terrain modeling sample area are: i. The cloud cover mark is less than 3; ii. The city flag is equal to 0; iii. The quantity mark is greater than 50; iv. The difference between the laser point elevation and the ASTER GDEM elevation value is less than 75; v. The difference between the average elevation of the laser point in each 100-meter statistical unit and the elevation value of the laser point interpolation surface is less than 75.

3. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S2, the rational polynomial function model is as follows: in, represents the coordinates of the image point on the image plane, Indicates the image point at the corresponding ground coordinate, Represents the rational polynomial function model to be solved.

4. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S2, absolute orientation is performed using the high-quality points obtained by screening as shown in the following formula: in, Indicates the 7 absolute orientation parameters that need to be solved; Represents the ground coordinate value of the model point; Represents the coordinates of the same model point in the auxiliary coordinate system of the image space; Represents the ground coordinate value of the model point; Indicates the model scaling factor; Indicates the photo posture angle The direction cosines are calculated as follows: in, Indicates the photo pose angle.

5. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S4, the method for extracting the valley line is specifically as follows: the cumulative flow of each pixel is calculated based on the DEM with the cloud area cut out, and then the cumulative flow threshold is used for segmentation, and pixels larger than the threshold are considered to be valley line pixels.

6. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S4, the ridgeline extraction method is specifically as follows: inverting the DEM with the cloud area cut out, then calculating the flow accumulation of each pixel based on the inverted DEM, and then using the flow accumulation threshold for segmentation, and considering pixels greater than the threshold as ridgeline pixels.

7. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S4, the method for extracting the ditch line is specifically as follows: perform a 30*30 window analysis on the DEM where the fog area is cut out to obtain a statistical average value to obtain an average DEM; then use the DEM where the fog area is cut out to subtract the average DEM to extract the area where the value is less than 0; then extract the boundary line of the area where the value is less than 0 to obtain the ditch line.

8. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In steps S5 and S8, the conditional generative adversarial network CGAN includes a generator network and a discriminant network. The generator network is used to generate synthetic data from random noise and conditional information, and its network structure adopts a U-Net structure; the discriminant network is used to distinguish between data generated by the generator and real data, and its network structure only includes encoding operations, does not include decoding operations, and the activation function adopts tanh.

9. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In steps S6 and S7, the resampling method is set to the nearest neighbor sampling method.

10. The terrain modeling method based on multi-source satellite remote sensing data according to claim 1, characterized in that: In step S9, the specific method for fusing the GF-7 DEM in the cloudless area with the DEM in the foggy area is as follows: a buffer zone with a width of 5 times the resolution is generated at the junction of the GF-7 DEM in the cloudless area and the DEM in the foggy area, the elevation value in the buffer zone uses the average value of the GF-7 DEM in the cloudless area and the DEM in the foggy area, the cloudless area outside the buffer zone uses the elevation value of the GF-7 DEM, and the foggy area outside the buffer zone uses the elevation value reconstructed by the conditional generative adversarial network (CGAN).

Images