CN112036264A

CN112036264A - Automatic extraction method of superglacial moraine covering type glacier

Info

Publication number: CN112036264A
Application number: CN202010812775.4A
Authority: CN
Inventors: 杨成生; 李海; 惠文华; 张勤
Original assignee: Changan University
Current assignee: Changan University
Priority date: 2020-08-13
Filing date: 2020-08-13
Publication date: 2020-12-04
Anticipated expiration: 2040-08-13
Also published as: CN112036264B

Abstract

The invention discloses an automatic extraction method of superglacial moraine covering type glacier, which comprises the following steps: step 1, preprocessing an optical image of a research area to obtain multispectral wave bands and textural features of the research area; step 2, selecting proper terrain data to perform terrain analysis to obtain a slope, a plane curvature and a section curvature corresponding to a research area; step 3, calculating the flow displacement characteristics of the research area by using an optical offset method; and 4, step 4: and classifying by adopting a random forest algorithm and obtaining the superglacial moraine profile. The method fully utilizes the characteristic that the superglacial moraine has fluidity, uses an optical offset technology to obtain the flow displacement, and combines a random forest algorithm to extract the superglacial moraine, so that the automation degree and the extraction accuracy are high; the random forest algorithm solves the problems of difficult determination of characteristic threshold values, characteristic use sequence and combination, and compared with a decision tree, the classification accuracy is improved; the addition of the flow rate greatly improves the accuracy of the superglacial moraine boundary extraction.

Description

Automatic extraction method of superglacial moraine covering type glacier

Technical Field

The invention belongs to the field of researching modern glacier by using an optical remote sensing technology, and relates to an automatic extraction method of superglacial moraine covering type glacier.

Background

Glacier changes are considered to be sensitive indicators of climate change, as well as the most important fresh water resources on land. In recent years, with the background of global warming, the glaciers in himalaya regions called "asia watertowers" have attracted continuous attention as an important component of global mountain glaciers. Optical remote sensing technology has been playing an important role in studying glacier changes on a large scale. The superglacial moraine covering type glacier is used as an important component of the mountain glacier, and the automatic extraction of the superglacial moraine covering type glacier has important significance for researching remote sensing change of the glacier. However, due to the spectral similarity between the superglacial moraine covering glacier and the surrounding bare ground and the complexity of the superglacial moraine covering glacier, the automatic extraction of the superglacial moraine glacier is difficult to overcome.

In the conventional glacier change research, the automatic extraction of threshold classification is only for neat glacier, but for the extraction of superglacial glacier covering type glacier, a visual interpretation mode of an expert is frequently adopted, and the method is time-consuming and labor-consuming for the large-area glacier change research. Subsequently, the scholars gradually introduce the characteristics of gradient, temperature, surface curvature, texture, coherence and the like into the automatic extraction process of the superglacial moraine, so that the development of the automatic extraction of the superglacial moraine is promoted.

Haeber li (1986) suggested that it was also difficult to distinguish the boundaries of the superglacial moraine and the surrounding bare land, sometimes simply by visual interpretation. Frank Paul (2003) suggests that relying solely on multispectral information would make it difficult to find the boundary between the superglacial glacial ice and surrounding debris and gradually introduce topographical factors into the extraction of the superglacial ice.

Song Bo (2006) and the like semi-automatically extract superglacial moraine in Gomba by using 5 factors of NDSI, NDVI, gradient, curvature and thermal infrared band generated by Landsat images and DEM and adopting a method combining threshold classification and spatial analysis. However, the threshold classification method makes this method geographically limited and questionable for the choice of threshold values, for example the authors refer to regions with a gradient greater than 24 ° as non-superglacial regions.

The jiang religious (2012) and the like judge the boundary of the superglacial moraine according to the fact that the superglacial moraine has fluidity and shows low coherence in the SAR image interferogram. STEFAN LIPPL (2018), and a morphological filtering post-treatment process is added on the basis, so that the single superglacial moraine is automatically extracted. The coherence map generated using the InSAR technique can be affected by a number of factors. For example, the time-space baseline is long, the topography is large, and atmospheric errors cause serious phase loss phenomena, so that the extraction result depends on many manual processes. And due to the difference of the imaging principle, the result of the SAR image and the optical image introduces registration error.

Wu\281567adopts a method of combining an object-oriented waveband ratio method and a characteristic partition to classify glaciers in Bomi county. And remote sensing image characteristic classification divides the research area into a neat glacier area, a superglacial moraine area and a glacier area under a shadow under the illumination condition. And a plurality of characteristics of wave band ratio, vegetation index, water body index, gradient and texture are quoted in the period for threshold classification. The threshold is artificially determined by trial-and-error and induction methods. And indicates that the neat glacier classification accuracy is 100% and the superglacial moraine is 93%. However, multiple experiments show that the accuracy of the threshold classification method for classifying the large-area superglacial moraine is difficult to reach 93%, and the article does not give a comparison graph of the extraction result of the superglacial moraine, so that the validity of the method remains uncertain.

In conclusion, the studies on the automatic extraction of tilde covering glacier have been greatly progressed by the efforts of a large number of scholars. However, as for the result, the problems that the accuracy of the extraction result is low, and the threshold value of each feature and the sequence of using the features are difficult to determine still exist, so that the methods cannot be popularized or the automatic extraction of the large-area superglacial moraine cannot be completed.

Disclosure of Invention

Aiming at the defects in the prior art, the invention aims to provide an automatic extraction method of superglacial moraine covering type glacier, and solves the problems that in the prior art, the extraction result is not high in accuracy, the automation degree is not high, the threshold value is difficult to determine, and the method is difficult to be applied to large-area superglacial moraine extraction.

In order to solve the technical problems, the invention adopts the following technical scheme:

an automatic extraction method of tillite covering glacier comprises the following steps:

step 1: preprocessing an optical image of a research area to obtain a multispectral wave band; performing principal component analysis on the multispectral wave bands to obtain a first principal component with the most information content, and performing texture calculation by taking the first principal component as an input wave band to obtain a collaborative wave band serving as a texture feature of a research area; the optical image of the research area comprises a thermal infrared band;

step 2: selecting appropriate terrain data for terrain analysis: projecting the terrain data into a projection zone which is the same as the multispectral wave band of the image, and further obtaining the corresponding slope, plane curvature and section curvature of the research area through a 3D surface analysis tool;

and step 3: calculating the flow displacement characteristics of the investigation region by using an optical offset method:

step 3.1, respectively taking optical images with a difference of 1-2 years between the front time phase and the rear time phase of the research region as input data for calculating the optical offset, setting the size of a search window, the step length, the signal-to-noise ratio threshold value and the robust iteration value, and calculating to obtain preliminary EW and NS horizontal displacement images of the research region;

step 3.2, respectively carrying out error correction on the preliminary EW and NS direction horizontal displacement images obtained in the step 3.1 to obtain accurate EW and NS direction horizontal displacement images, and synthesizing the horizontal flow displacement images by a wave band calculation tool to obtain the flow displacement characteristics of the research area;

and 4, step 4: classifying by adopting a random forest algorithm and obtaining a superglacial moraine profile:

step 4.1, manually selecting samples on the image of the research area, wherein the samples comprise nude land, clean glacier, superglacial moraine and ice lake four types, and the samples are uniformly distributed in the image of the research area;

step 4.2, synthesizing the multispectral wave bands, the thermal infrared wave bands, the texture characteristics, the gradient, the plane curvature, the section curvature and the flow displacement characteristics obtained in the step 1, the step 2 and the step 3 into a multiband image;

4.3, inputting the samples in the step 4.1 and the multiband images in the step 4.2 into a random forest classification algorithm for learning and modeling; inputting the part which is not selected in the image of the research area into a random forest algorithm for automatic classification, and then outputting a classification result;

and 4.4, performing morphological classification post-treatment on the classification result to obtain the accurate superglacial moraine contour.

The invention also comprises the following technical characteristics:

specifically, the preprocessing in the step 1 includes atmospheric correction and image stitching;

the texture calculation in the step 1 adopts a texture calculation method based on second-order probability statistical filtering.

Specifically, the step 3.2 error correction includes long wavelength track error and streak artifact error;

correcting the long wavelength orbit error by a polynomial curve fitting method; the streak artifact error is corrected by the mean subtraction principle.

Specifically, the step 4.4 of performing morphological classification post-processing on the classification result comprises the following steps:

4.4.1, removing classification holes and deleting isolated points in the superglacial moraine region by using a Majority/least analysis algorithm for two times or more;

step 4.4.2, performing cluster analysis by using a column Classis algorithm;

4.4.3, carrying out vectorization and smoothing operation on the superglacial moraine to obtain a preliminary superglacial moraine contour vector;

and 4.4.4, screening through a certain area gradient to obtain the accurate superglacial moraine profile.

Compared with the prior art, the invention has the following technical effects:

the method fully utilizes the characteristic that the superglacial moraine has fluidity, uses an optical offset technology to obtain the fluidity displacement, and combines a random forest algorithm to extract the superglacial moraine; in the whole process, besides manual sample selection, the automation degree is high, and the extraction accuracy is high; the random forest algorithm solves the problems of difficult determination of characteristic threshold values, characteristic use sequence and combination, and compared with a decision tree, the classification accuracy is improved. The addition of the flow displacement greatly improves the accuracy of the superglacial moraine boundary extraction.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 shows a typical superglacial moraine on both sides of the test area (Landsat4, 3, 2 combination) of step 1 of the present invention;

FIG. 3 illustrates the first principal component band and texture features of step 1 of the present invention;

FIG. 4 shows topographic features of step 2 of the present invention, where A is DEM, B is slope, C is section curvature, and D is plane curvature;

FIG. 5 is a step 3 glacier flow displacement strip processing of the present invention; removing a pre-trend error image; removing a trend factor; c, removing the trend of the image; d, removing the strip to obtain an image;

FIG. 6 shows the results of the classification in step 4.2, where the red region is the classification result and the green line region is the manually demarcated superglacial moraine region;

FIG. 7 is the superglacial moraine extraction result after the area gradient screening and vectorization of step 4.3;

FIG. 8(a) is a profile obtained by a conventional classification method; (b) is a contour map obtained by random forest classification in the application;

FIG. 9(a) is a result diagram with no shift feature engaged; (b) is a result graph with the participation of displacement characteristics in the application;

FIG. 10 is a comparison of the superglacial acetic acid extraction area with the manually delineated area.

The present invention will be described in further detail with reference to the accompanying drawings and examples.

Detailed Description

The flowing speed of the mountain glacier in the Himalayan region is 2-50m/y, and the flowing speed of the superglacial glacier is about 0.5-10m/y, which belongs to larger displacement. Through multi-party verification, the precision of the optical offset technology is generally about 1/20 pixels, and the precision of the landsat image for measuring the glacier flow velocity is only about 0.75-1.5m, so that the optical offset technology is feasible for measuring the superglacial moraine flow velocity. The fluidity is an important factor for distinguishing the glacier from the surrounding bare land, so that the addition of the flow rate factor (namely, the flow displacement) can theoretically improve the extraction precision of the superglacial moraine.

A large number of researches show that the random forest algorithm has the characteristics of strong model generalization capability, low requirement on the number of samples and automatic determination of feature importance. The method solves the problems that the threshold is difficult to determine and the sequence and importance of the features are difficult to determine in the conventional classification method. The sample requirement is less, the generalization ability is strong, the universality of the classification method is enhanced, and the automatic extraction of the superglacial moraine can adapt to a wider area.

The invention provides an automatic extraction method of superglacial moraine covering type glacier, which comprises the following steps:

step 1: preprocessing an optical image of a research area, including atmospheric correction and image splicing; the optical image of the research area comprises a thermal infrared band; the preprocessed image comprises multispectral wave bands; performing principal component analysis on the multispectral wave bands of the preprocessed image to obtain a first principal component PC1 with the largest information content, and performing texture calculation by taking the first principal component PC1 as an input wave band to obtain a collaborative wave band as a texture feature of a research region; the texture calculation in the step 1 adopts a texture calculation method based on second-order probability statistical filtering.

Specifically, the embodiment preprocesses the Landsat image in the 2018 research area, the Landsat image provided by the USGS is in the level of L2, and the image is subjected to orthorectification and geometric rectification, so that the preprocessing only includes atmospheric rectification and image stitching; then, performing principal component analysis on the multispectral wave bands to obtain a first principal component PC1 with the largest information content, performing texture analysis on the first principal component PC1, and selecting a texture calculation method based on second-order probability statistical filtering; selecting a synergistic wave band for texture analysis through comparison analysis, wherein the synergistic wave band can best distinguish the superglacial moraine; thus obtaining three groups of factors of multispectral 6 wave bands, thermal infrared wave bands and texture wave bands together; wherein the texture calculation is done in the corresponding Co-occurrence measures tool in ENVI 5.3. Fig. 3(a) is a first principal component image after principal component analysis is performed on multispectral bands of the Landsat image; fig. 3(b) shows a texture feature image calculated from the first principal component.

Step 2: selecting appropriate terrain data for terrain analysis: projecting the terrain data into a projection zone which is the same as the multispectral wave band of the image, and further obtaining a slope map, a plane curvature and a section curvature which correspond to a research area through a 3D surface analysis tool;

specifically, in order to obtain terrain data which is more in line with reality, ALOS-World-3DDSM data is adopted; the data is high-precision global digital earth model data freely provided by the Japanese space aviation research and development organization (JAXA)2015 in 5 months, the horizontal resolution is 30 meters, the elevation precision is 5 meters, and the data is one of the most accurate data in the world at present. The dsm is projected in arc-gis into the same projection band as the multi-spectral band WGS 84-45N. The corresponding slope map, plane curvature, section curvature were obtained using the 3D surface analysis tool in arc-gis. As shown in FIG. 4, A is a research area DEM; b, generating a gradient image by the DEM; c, generating a plane curvature image by the DEM; d, generating a section curvature image by the DEM.

And step 3: calculating the flow displacement of the research area by using an optical offset method;

step 3.1, setting the size and the step length of a search window according to the resolution of the optical image; for example, the method verifies that the resolution of the sentinel-2 image used in the experiment is 10m, and after multiple experiments, a search window is selected to be 64 multiplied by 64, and the step length is selected to be 4; the Landsat series image is selected to be 32 x 32, the step length is 4, and the effect is best; in addition, the signal-to-noise ratio threshold is selected to be 0.9, and the robust iteration value is selected to be 2; therefore, preliminary EW and NS direction horizontal displacement images of the research area can be obtained;

and 3.2, respectively carrying out error correction on the primary EW and NS direction horizontal displacement images obtained in the step 3.1 to obtain accurate EW and NS direction horizontal displacement images, and synthesizing the horizontal flow displacement images by a wave band calculation tool to obtain the flow displacement characteristics of the research area.

Because of satellite orbit errors and imaging errors, the horizontal displacement image obtained in the previous step contains many errors, such as long-wavelength orbit errors and streak artifacts, which have a large influence surface, and therefore, it is necessary to correct errors in the horizontal displacement images in the two directions separately. Error correction includes long wavelength orbital error and streak artifact error; the long wavelength orbit error is corrected by a polynomial curve fitting method; the streak artifact error is corrected by the mean subtraction principle. Wherein the long wavelength orbit error can be well corrected by a polynomial curve fitting method; the method can be realized by an image debugging tool provided in COSI-CORR software; aiming at the streak artifact errors, the Destripe image tool provided by COSI-CORR software has great limitation, so that the method adopts the mean value subtraction principle proposed by Von light finance of the university of Zhongnan and is realized by matalab language programming, and the result has a good streak removal effect. And finally, synthesizing a horizontal flow dynamic displacement diagram by a waveband calculation tool in the ENVI, wherein the formula is as follows: sqrt (b1^2+ b2^ 2).

Specifically, the flow displacement of the research area superglacial moraine covering glacier in 2016 and 2018 is calculated by using an optical offset technology; selecting COSI-CORR software with better processing result for processing through comparison; the software is developed by California institute of technology and engineering, based on IDL language, and is embedded in ENVI 5.3; aiming at the banding problem in the offset calculation result, a mean value subtraction method is adopted and the banding problem is realized through matlab programming; and finally, a more accurate flow displacement graph is obtained. As shown in FIG. 5, A, the horizontal image of the NS direction obtained from the right side of the study area; removing error factors of long-wave band orbit errors by an image detrending tool; c: removing the long-wave band orbit error and then horizontally shifting the image in the NS direction; d: and removing the NS-direction horizontal displacement image after the stripe artifact is removed by a mean subtraction method.

It can be seen from the graph A that the horizontal shift image without error processing has severe streaking, and the graph D can be seen that the streaking is well corrected.

the random forest algorithm belongs to a supervised classification algorithm in the field of machine learning, and samples need to be selected in advance, provided for algorithm learning and modeled. Aiming at the ground feature characteristics of the region where the glacier is located and aiming at improving the accuracy of a classification algorithm, samples are selected and classified into bare land, clean glacier, superglacial tillite and glacier, the samples are uniformly distributed in a research area, and the area of the samples only accounts for about 0.1% of the total area.

The whole classification process is carried out in ENVI5.3, samples are selected in a roi form, finally classification is carried out by using a random forest algorithm provided by an ENVI expansion tool, 80-100 of decision tree selection is carried out, other parameters are determined by default, and the processing process is about 15 minutes.

And 4.4, performing morphological classification post-treatment on the classification result to obtain the accurate superglacial moraine contour. FIG. 7 shows the results of vectorization, smoothing, and area gradient screening of classification results; wherein the red line (line 1) is the superglacial moraine visual interpretation profile and the yellow line (line 2) is the superglacial moraine extraction profile according to the method of the present invention.

Specifically, a sample is selected in the research area, the sample is divided into bare land, neat glacier, superglacial moraine and ice lake, the sample is uniformly distributed in the research area, and the area of the sample only accounts for about 0.1% of the total area. And classifying by utilizing an ENVI expansion tool random forest algorithm, wherein the processing process is about 15 minutes. And then carrying out morphological classification post-processing operation on the classification result, completing loopholes, deleting isolated points and clustering, and carrying out vectorization, smoothing and area gradient deletion to obtain the superglacial moraine contour.

The following embodiments of the present invention are provided, and it should be noted that the present invention is not limited to the following embodiments, and all equivalent changes based on the technical solutions of the present invention are within the protection scope of the present invention.

Example 1:

in this example, the moraine of Xixia Ma peak and the staggered Lanma glacier lake group region in the middle section of the Himala mountain was extracted by the above method using the Landsat-8 image Sentiniel-2 and ALOS-World-3ddsm data. The images of the area in 2016 and 10 months and 2018 and 10 months are downloaded from the usgs website. However, due to cloud amount and image inventory problems, in the area 2016, the Sentinel-2 images of 2016, 12 and 2018, 10 are adopted to obtain the superglacial moraine displacement of 2016, 2018. The time difference is mainly two years because the superglacial moraine flows more slowly than the neat glacial moraine, and the flow amount for two years is adopted to highlight the flowing characteristic and help the classification. It is noted that a time span that is too long can cause decorrelation.

In order to embody the superiority of the method, the same characteristics are adopted at the same time, decision tree threshold classification is carried out, each threshold is confirmed for multiple times through a trial and error method, the optimal classification result is basically achieved, and the classification result is shown in fig. 8 (a); meanwhile, in order to show the importance of the displacement features, random forest classification without the displacement features is performed under the same conditions, and the classification result is shown in fig. 9 (a).

FIG. 8: a is a decision tree threshold classification result (red line, i.e., line 3, is the superglacial moraine visual interpretation contour, green line, i.e., line 4, is the contour obtained by decision tree threshold classification), and b is a random forest classification result (red line, i.e., line 5, is the superglacial moraine visual interpretation contour, and yellow line, i.e., line 6, is the contour obtained by random forest classification); in the graph, the traditional decision tree threshold classification can be seen, even if the classification result is still poor under the assistance of displacement characteristics, the random forest classification can obtain a better classification result; the method has the advantage of high accuracy of classification results.

FIG. 9: a is a random forest preliminary classification result without participation of flow displacement characteristics; and b is the result of the preliminary classification of the random forest with the participation of the flow displacement characteristics (the method of the invention).

The comparison shows that the accuracy of the random forest classification result without the participation of the flow displacement features is improved relative to the threshold classification result of fig. 8, but a plurality of 'bugs' and large-area misclassifications exist. In the classification result with the flow displacement participating in the classification method, the phenomena of 'bugs' and wrong classification are greatly improved, so that the accuracy of the classification result is further improved.

From comparison in fig. 9, it is found that the automatic extraction accuracy of the random forest is greatly improved compared with the threshold classification result. The addition of the displacement characteristics greatly improves the adverse conditions of large and numerous gaps, inaccurate edge extraction and the like in the superglacial moraine extraction. The main reason is that the displacement has a certain continuity in the plane with respect to other factors. The problem of more gaps in the extraction result is well solved. Further, the fluidity is significantly different from that of the ambient stable bare table as an important feature of the superglacial moraine. Therefore, from the other areas, the addition of the displacement features excludes the false extraction and omission of some edge areas. Thereby further improving the extraction accuracy.

Fig. 10 is for verification of extraction accuracy. The range of the superglacial moraine in the area is manually delineated by combining the cataloguing data as a standard result. The large 15 glaciers were selected for comparison and found to have an average area accuracy of 89.7%. And a large amount of small superglacial moraine is found in the classification result, and the regions really have the characteristics of spectrum and fluidity of the superglacial moraine through manual careful interpretation and extraction of corresponding flow rate, so that the small regions are really superglacial moraine. The risk of missing small glaciers is proved to be possessed by artificially delineating the superglacial moraine glacier, and the method provided by the text can avoid most of the missing problems. Therefore, the method better solves the problems that the threshold value is difficult to determine, the feature combination is not good, the large-area applicability is not good and the automation degree is low in the automatic superglacial moraine extraction at present.

Claims

1. An automatic extraction method of superglacial moraine covering type glacier is characterized by comprising the following steps:

2. The automated superglacial moraine covering-type glacier extraction method according to claim 1, wherein the pretreatment in the step 1 comprises atmospheric correction and image stitching;

3. The automated superglacial moraine covering-type glacier extraction method of claim 1, wherein the step 3.2 error correction comprises long wavelength orbital error and streak artifact error;

4. The automated superglacial moraine covering-type glacier extraction method according to claim 3, wherein the step 4.4 of performing morphological classification post-processing on the classification result comprises the following steps:

step 4.4.2, performing cluster analysis by using a column Classis algorithm;