CN112700431A

CN112700431A - Flood coverage surface extraction method based on remote sensing image

Info

Publication number: CN112700431A
Application number: CN202110032720.6A
Authority: CN
Inventors: 黄文丽; 冯梅; 陈淑杰
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2021-01-11
Filing date: 2021-01-11
Publication date: 2021-04-23
Anticipated expiration: 2041-01-11
Also published as: CN112700431B

Abstract

The invention provides a flood coverage surface extraction method based on remote sensing images, which comprises the following steps: s1: acquiring optical remote sensing image data of a terrestrial satellite, calculating historical water body distribution probability of the optical remote sensing image data, and acquiring a perennial water area range distribution map by setting a distribution probability threshold; s2: acquiring satellite radar remote sensing image data, calculating statistical information of each pixel time domain within a set time range, and calculating a current water area distribution map by adopting a statistical threshold method; s3: and classifying the flood coverage types by combining the perennial water area and the current water area range, and extracting the current flood coverage ranges of different levels. The method provided by the invention can be used for rapidly and accurately acquiring the current flood coverage range by combining historical and current water body distribution.

Description

Flood coverage surface extraction method based on remote sensing image

Technical Field

The invention belongs to the technical field of remote sensing application, and particularly relates to a remote sensing image flood coverage surface extraction method based on water body distribution rate and a statistical threshold value.

Background

Flood is the most frequent natural disaster in the world and seriously threatens the life and property safety of people. Therefore, in order to rapidly analyze the disaster and provide rescue, rapid extraction of information about sudden flood is necessary. The satellite remote sensing image has great advantages in the aspect of extracting water as an important means for monitoring earth surface information in recent years.

Remote sensing data commonly used for large-scale flood assessment include optical and radar remote sensing data. The optical remote sensing data volume is large, the resolution ratio is high, and extensive research is carried out on the aspect of extracting the water area, but the water area is easy to be extracted and the cloud shadow image element can be mixed together. The radar remote sensing uses microwaves capable of penetrating through cloud layers, is not easily influenced by weather conditions, solves the problem of cloud interference in optical remote sensing, simultaneously has the characteristic of active remote sensing, enables the radar remote sensing to be imaged at night, and more importantly, in recent years, data of satellites such as Sentinel-1 and the like are freely disclosed to the world and are easy to obtain. Therefore, radar remote sensing data characterized by all-time and all-weather conditions has been increasingly used for evaluating flood disasters in these years.

A great deal of research and accumulation on the aspect of water area extraction of optical remote sensing data exist, and a mature threshold algorithm is used for extracting the water area range of an area. At present, various water body distribution rate products exist in the global range, and can be used as a basis for analyzing historical changes of the water body range. However, in the optical remote sensing image, cloud interference causes data to be easy to lose in a flood occurrence period, and key period distribution information cannot be obtained; the cloud shadow is easily confused with the water body or the wetland, and the accuracy of extracting the water body is influenced.

The radar remote sensing supports all-weather and all-weather information extraction, is not easily interfered by cloud layers in sudden flood areas, and has wide application prospect. The radar remote sensing image can reflect the reflection of the ground object to the electromagnetic wave and the heat radiation information of the ground object. For different ground objects, due to the difference of the internal structure composition and the physicochemical property, the reflection of electromagnetic waves and the heat radiation of the ground objects are different, so that the method can be used as a means for rapidly distinguishing the ground objects. However, due to the late development time, the radar satellite data is less in data volume accumulation, and the method is more suitable for range extraction of a water area during a sudden flood event.

It can be seen that a flood area with high accuracy cannot be obtained by using optical remote sensing data or radar remote sensing data alone, and therefore, how to more accurately extract the flood coverage area is the focus of current research.

Disclosure of Invention

The invention aims to provide a flood coverage area extraction method based on remote sensing images aiming at the defects of the prior art, the method combines history and current water body distribution to extract the flood coverage area range, and the current flood coverage area can be quickly and accurately obtained.

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

a flood coverage surface extraction method based on remote sensing images comprises the following steps:

s1: acquiring optical remote sensing image data of a terrestrial satellite, calculating historical water body distribution probability of the optical remote sensing image data, and acquiring a perennial water area range distribution map by setting a distribution probability threshold;

s2: acquiring satellite radar remote sensing image data, calculating statistical information of each pixel time domain within a set time range, and calculating a current water area distribution map by adopting a statistical threshold method;

s3: and classifying the flood coverage types by combining the perennial water area and the current water area range, and extracting the current flood coverage ranges of different levels.

Further, the specific method of step S1 is: extracting a wave band with DSWE model characteristics from optical remote sensing image data, calculating each pixel in the extracted wave band by adopting a DSWE model to obtain the distribution range of perennial waters with different levels, setting a distribution probability threshold, counting the probability that the perennial waters with different levels in a historical period are water bodies, and comparing the probability with the probability threshold to determine the range distribution of perennial non-waters and perennial waters.

Further, step S1 includes the following sub-steps:

s11: analyzing indexes used in the DSWE model and a calculation formula of each index;

s12: defining and calculating each index of each pixel in the extracted wave band by adopting a calculation formula of each index in the DSWE model;

s13, comparing whether each index obtained by calculating each pixel meets the condition according to the condition that each index in the DSWE model needs to meet, if so, taking 1 as the result, and if not, taking 0 as the result, and then respectively placing the results on ten million bits in sequence to enable each pixel to obtain a five-bit code;

s14, determining a classification standard, classifying the codes obtained by each pixel according to the corresponding classification standard to distinguish non-water bodies, water bodies with middle confidence level, water bodies with high confidence level and partial surface water bodies to obtain the classification result of the perennial water area;

s15: and counting the probability that each pixel is a water body in a historical time period, and setting a probability threshold value to determine the range distribution of perennial non-water areas and water areas.

Further, the probability threshold is set to 90%, that is, the probability that a certain pixel is a water body in a historical period is greater than 90%, and the pixel is defined as an perennial water body.

Further, the specific method of step S2 is:

s21, calculating backscattering coefficients of the satellite radar remote sensing image data in different polarization modes, orbit directions and acquisition modes;

s22, calculating an average backscattering coefficient and a standard deviation backscattering coefficient of a period of no flood according to historical hydrological conditions;

s23, substituting the calculation result into a Z-score threshold value method calculation formula to obtain Z of the SAR backscattering of the vertical emission vertical reception of the flood in the required period_VVScore value distribution plot and Z of vertical emission horizontal re-emission SAR backscatter_VHScore value profile.

Further, step S3 further includes the following sub-steps:

s31: given Z_VVScore threshold and Z_VHA score threshold, finding all pixels not in the perennial water area in the distribution map of the perennial water area obtained in step S1, and obtaining Z calculated in step S2 corresponding to all pixels_VVScore and Z_VHScore value distribution plot, Z for each pixel_VVScore and Z_VHScore value with given Z_VVScore threshold and Z_VHScore threshold is compared if Z of pel_VVScore lower than given Z_VVScore threshold and Z thereof_VHScore lower than given Z_VHScore threshold, then this pel is marked as a high confidence flood label; if Z is_VVScore and Z_VHOnly one of score is below a given threshold, this pel is marked as a medium confidence flood label; if Z is_VVScore and Z_VH-score is above a given threshold, respectively, then this pel is marked as non-flood;

s32: determining special frequent and sudden flood areas; in high-confidence flood, when the total flooding probability is more than 25%, defining the pixel as a frequent flood area; in the flood with medium confidence coefficient, when the total flooding probability is more than 25%, defining the pixel as the frequent flood; in non-flood, when the total flooding probability is greater than 25%, the pixel is defined as a sudden flood.

Further, given Z_VVScore threshold and Z_VHThe score threshold is determined by:

i. deriving the calculated Z_VVScore value distribution plot, Z_VHScore value distribution graph and radar remote sensing image of corresponding area:

ii. Randomly selecting some points in the radar remote sensing image and obtaining the point Z_VVScore value, Z_VH-finding the Z corresponding to it in the score value image_VVScore and Z_VH-score value and judging the category of the point, namely the point is a water body or a non-water body, and generating a scatter diagram;

iii, analyzing the scatter diagram to find Z for distinguishing the water body from the land_VVScore value and Z_VH-a range of score values;

iv, mixing Z_VVScore value and Z_VH-refinement of the range of score values and thinning of the Z_VVScore and Z_VH-score groups of different values are combined, calculated in relation to the precision verification image, all groups calculatedThe combined overall classification precision is found, and the combination which enables the overall classification precision to reach the maximum is used as the Z extracted in the current water area range_VVScore and Z_VH-score threshold.

Compared with the prior art, the invention has the beneficial effects that:

1) according to the method, the land satellite image is converted into historical water body distribution information, the current water body distribution is extracted based on the radar image by using a statistical threshold method, and the flood coverage area is extracted by combining the historical water body distribution and the current water body distribution, so that the whole algorithm idea is clear and intuitive, and the current flood coverage area can be rapidly and accurately obtained;

2) the method combines the advantages of the optical remote sensing data and the cloud layer, so that the obtained flood area range is more accurate; meanwhile, the research test is carried out based on a cloud computing platform, and data can be directly obtained and operated, so that the complex data downloading and processing operation can be greatly simplified, and the flood coverage range can be quickly extracted in response; in addition, the portability of the algorithm is high, a user does not need to download software, only an effective cloud computing platform number is needed, the execution code can be logged in a webpage, the operation result can be obtained, and the algorithm can also operate the code on a high-performance stand-alone platform after simple adjustment; in conclusion, the invention has stronger universality and is convenient to popularize and use.

Drawings

FIG. 1 is a flow chart of a method for flood coverage extraction according to an embodiment of the present invention;

FIG. 2 is a series of diagrams of a current water area extraction process according to an embodiment of the present invention; wherein, fig. 2(a) is a distribution diagram of the VV average backscattering coefficient, fig. 2(b) is a distribution diagram of the VV standard deviation backscattering coefficient, fig. 2(c) is a distribution diagram of the VV backscattering coefficient on the day of flooding of a certain city, and finally, fig. 2(d) is a distribution diagram of the calculated Z_VV-score profile;

FIG. 3 is a schematic diagram of flood classification according to an embodiment of the present invention;

FIG. 4 is a series of graphs illustrating the determination of threshold values according to an embodiment of the present invention; FIG. 4(a) shows a radar remote sensing image (for accuracy evaluation)(ii) a FIG. 4(b) is based on Z_VVScore and Z_VHDividing a water body and non-water body scatter diagram by the score value;

fig. 5 is a graph of flood extraction results according to an embodiment of the present invention.

Detailed Description

The technical solutions in the embodiments of the present invention will be described clearly and completely with reference to the following embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

It should be noted that the embodiments and features of the embodiments may be combined with each other without conflict.

The present invention is further illustrated by the following examples, which are not to be construed as limiting the invention.

The invention provides a flood coverage surface extraction method based on remote sensing images, which comprises the following steps:

s1: acquiring Landsat-8 land satellite optical remote sensing image data, calculating historical water body distribution probability of the data, and acquiring a perennial water area range distribution map by setting a distribution probability threshold;

in the step, terrestrial satellite optical remote sensing image data are obtained, a characteristic waveband required by calculation of a DSWE model is selected from the image data, and then the DSWE model is adopted for calculation. During calculation, firstly analyzing indexes used in the DSWE model and calculation formulas of the indexes, and then defining and judging the indexes by the calculation formulas of the indexes in the DSWE model for each pixel in the extracted wave band; the DSWE model is composed of calculation formulas of five indexes, wherein the formulas are as follows:

1)MNDWI＞123[scaled by 10000]

2)MBSRV＞MBSRN

3)AWEsh＞0

4)MNDWI＞-5000&SWIR1＜1000&NIR＜1500

5)MNDWI＞-5000&SWIR2＜1000&NIR＜2000

the MNDWI is a corrected normalized difference water index calculated according to the ground reflectivity; MBSRV is the sum of GREEN and RED, MBSRN is the sum of NIR and SWIR 1; AWEsh is the automatic water extraction index; SWIR1 is the shorter of the sensor short wave infrared band; SWIR2 is the longer of the sensor short wave infrared bands; NIR is the near infrared band; the calculation formula of the index is as follows:

MNDWI＝(GREEN-SWIR)/(GREEN+SWIR)*10000

AWEsh＝BLUE+2.5GREEN-1.5MBSRN-0.25SWIR2

MBSRV＝GREEN+RED

MBSRN＝NIR+SWIR1

in the above formula, GREEN is the GREEN band, BLUE is the BLUE band, and RED is the RED band.

Calculating each pixel according to the calculation formula of the 5 indexes, comparing the calculation result with the conditions which need to be met by the five indexes in the DSWE model, wherein the condition which is met is 1, the condition which is not met is 0, then respectively placing the results on ten million bits in sequence to enable each pixel to obtain a five-bit number, and then comparing the five-bit number with the division standard to obtain the attribute of the pixel: whether the water is high confidence water or middle confidence water or non-water or partial surface water, the division standard is as follows:

1100111111: 1 (high confidence water)

10111 10999：1

01111 01111：1

1100011000: 3 (partial surface water)

10000 10000：3

01000 01000：3

1001210110: 2 (middle confidence water)

10011 10011：2

10001 10010：2

01001 01110：2

00010 00111：2

0000000009: 0 (nonaqueous)

Extracting perennial water areas and perennial flood areas, specifically extracting 16 years of land satellite optical remote sensing images from 2000 to 2016, calculating by using the DSWE method, and defining any pixel as a perennial water area if the probability that the pixel is a water body (including high-confidence water, middle-confidence water and partial surface water) within 16 years is more than 90%; and if the probability of the pixel being a water body (including high-confidence water, middle-confidence water and partial surface water) is more than 25%, defining the pixel as a frequent flood area, and obtaining images of the perennial water area and the frequent flood.

S2: acquiring remote sensing image data of a Sentinel-1 satellite radar, calculating statistical information of each pixel time domain within a set time range, and calculating a current water distribution map by adopting a statistical threshold method; the step further comprises the substeps of:

Previous studies have shown that water has a strong absorption of incident energy and generally exhibits a weak reflectivity over most of the wavelength range of remote sensors. Thus, for bodies of water and land, the VH, VV backscattering values may differ, i.e. flooding may have an effect on the Z-score values of the vertical transmit vertical receive (VV) and vertical transmit horizontal re-transmit (VH) SAR backscattering, which are calculated as:

wherein, p is a polarization mode, m is a sensor acquisition mode, d is a track direction, the acquisition mode comprises IW and SM, and the track directionThe direction of the movement includes the rising and the falling,

represents the average backscattering coefficient when the polarization mode is p, the sensor acquisition mode is m and the track direction is d,

representing the backscatter coefficients for a polarisation mode p, a sensor acquisition mode m and a track direction d at time t,

the standard deviation backscattering coefficient is shown when the polarization mode is p, the sensor acquisition mode is m and the track direction is d. The Z-score value of SAR backscattering with polarization mode of vertical transmission and vertical reception is expressed as Z_VV-score, Z-score value of SAR backscattering with polarization mode being vertical emission horizontal re-emission expressed as Z_VH-score。

According to hydrological data of a certain city in historical period, the probability of flood generation in 3 to 5 months in the past year is low. By consulting hydrological data of the water and utility bureau of a certain city in 2019, it is verified that no flood occurs between 3 months and 1 days, and therefore 3-5 months are selected as the limit of the historical baseline of the certain city in the research.

The average backscattering coefficient for different polarization, orbit, and acquisition modes was then calculated using all the Sentinel-1 images acquired between these dates

The calculation result is shown in FIG. 2(a), the standard deviation backscattering coefficient

The calculation result of (c) is shown in fig. 2 (b). In 2019, 6 and 14 days, flooding a certain city, and calculating the VV backscattering coefficient of the current day

The calculation results are shown in fig. 2 (c). Finally, according to the calculation formula, the calculation can be carried outZ_VVThe results of the calculation of-score are shown in FIG. 2 (d). In the same way, we can also calculate Z_VH-score value.

S3: classifying the flood coverage types by combining the perennial water area and the current water area range, and extracting the current flood coverage ranges of different levels; the step further comprises the substeps of:

In this embodiment, we can use the GEE code to implement the above classification process, and assign different values to different types of areas, and use different colors to represent, for example, 0 represents land, 1, 2 represents medium confidence flood, and 3 represents high confidence flood (burst); 10 represents the usual flood, but the monitoring extraction is non-flood, 11, 12 represents the mid confidence flood, and the usual flood before, 13 represents the usual high confidence flood (usual); 20 represents perennial waters. To facilitate understanding of the partitioning process, the partitioning manner is summarized into a two-dimensional matrix legend, wherein the flooding confidence based on SAR constitutes one axis, and the historical flooding information based on Landsat constitutes another axis (lower left corner of fig. 3).

And in step S31, Z is given_VVScore threshold and Z_VHThe score threshold is determined by:

i. deriving the calculated Z_VVScore value distribution plot, Z_VH-score value distribution map and radar remote sensing image of corresponding area;

iii, analyzing the scatter diagram to find Z for distinguishing the water body from the land_VVScore and Z_VH-a range of score;

iv, mixing Z_VVScore value and Z_VH-refinement of the range of score values and thinning of the Z_VVScore value and Z_VHCombining a plurality of groups of different values of the score value, calculating a contact precision verification image, calculating the overall classification precision of all the combinations, and finding the combination which enables the overall classification precision to be maximum to be used as Z extracted in the current water area range_VVScore and Z_VH-score threshold.

Specifically, in the present embodiment, fig. 4 is a series of threshold determination charts, and Z is a series of threshold determination charts in the process of extracting flood water_VVScore and Z_VHAccuracy of score threshold selection has a large impact on accuracy, and to select the highest discrimination threshold, the study downloaded Sentinel-2 images from 2019, 6, 14 days (fig. 4(a)), and extracted Z for that day_VVScore and Z_VHComparison of the score images allows to find the Z distinguishing between bodies of water and non-bodies of water_VVScore and Z_VHApproximate range of score, the procedure is as follows:

with the export function of GEE, the calculation will be performedZ obtained_VV-score and_ZVH-the score image is downloaded, together with the Sentinel-2 satellite data with the corresponding area resolution of 10 m.

Opening the three pictures in ENVI at the same time, randomly selecting 40 points by using the random sampling and cursor positioning functions of the three pictures, and recording the Z points in an excel table_VVScore and Z_VHScore value and visually judge the category to which the point belongs (water or non-water, water assigned 1, non-water assigned 0), and generate a scatter plot, as shown in fig. 4 (b).

Analyzing the scatter diagram, and visually finding Z which can better distinguish water from land_VVScore value and Z_VH-score value in the approximate range (-2 to-3).

Will Z_VVScore and Z_VH-score threshold is refined and divided into-2, -2.5, -3 sections respectively to form nine combinations, different threshold combinations are selected, flood extraction results are obtained by running in GEE, accuracy verification images are connected, 250 points are randomly generated, the actual type (visual judgment whether the pixel is water or land) and the extraction type (judgment according to the pixel value) of each point are obtained by using an ENVI cursor positioning function, the results are recorded in an excel table, the overall classification accuracy is calculated to obtain overall classification accuracy tables with different thresholds, see table 1, and the threshold (in the embodiment, when Z is used as the threshold) which enables the overall classification accuracy to be maximum is found out_VVScore and Z_VHWhen scores are all-2 and all-2, the obtained overall classification accuracy is maximum) as Z extracted in the current water area range_VVScore and Z_VH-score threshold.

TABLE 1 is Z_VVScore and Z_VHOverall classification accuracy from different combinations of score

FIG. 5 is a flood extraction graph of Jian City, according to the method described above, using a determined threshold (Z)_VVScore and Z_VHScore is-2), flood conditions of 2019, 9 and 14 days are extracted from GEE, and the extraction result is shown in figure 5.

While the invention has been described with reference to a preferred embodiment, it will be understood by those skilled in the art that various changes in form and detail may be made therein without departing from the spirit and scope of the invention.

Claims

1. A flood coverage surface extraction method based on remote sensing images is characterized by comprising the following steps:

2. The method for extracting flood coverage based on remote sensing images according to claim 1, wherein the specific method in step S1 is as follows: extracting a wave band with DSWE model characteristics from optical remote sensing image data, calculating each pixel in the extracted wave band by adopting a DSWE model to obtain the distribution range of perennial waters with different levels, setting a distribution probability threshold, counting the probability that the perennial waters with different levels in a historical period are water bodies, and comparing the probability with the probability threshold to determine the range distribution of perennial non-waters and perennial waters.

3. The method for extracting flood coverage based on remote sensing images as claimed in claim 2, wherein the step S1 further comprises the following sub-steps:

4. The method for extracting flood coverage based on remote sensing images as claimed in claim 3, wherein the probability threshold is set to 90%, that is, the probability that a certain pixel is a water body in a historical period is greater than 90%, and the pixel is defined as a perennial water body.

5. The method for extracting flood coverage based on remote sensing images according to claim 1, wherein the specific method in step S2 is as follows:

6. The method for extracting flood coverage based on remote sensing images as claimed in claim 5, wherein the step S3 further comprises the following sub-steps:

s32: determining special frequent and sudden flood areas; in high-confidence flood, when the total flooding probability is more than 25%, defining the pixel as a frequent flood area; in the flood with medium confidence coefficient, when the total flooding probability is more than 25%, defining the pixel as the frequent flood; in non-flood, when the total flooding probability > 25%, then the pixel is defined as a sudden flood.

7. The method of claim 6, wherein the given Z is the flood coverage area_VVScore threshold and Z_VHThe score threshold is determined by:

i, deriving the calculated Z_VVScore value distribution plot, Z_VH-score value distribution map and radar remote sensing image of corresponding area;

ii, randomly selecting some points in the radar remote sensing image and enabling the points to be located in Z_VVScore value, Z_VH-finding the Z corresponding to it in the score value image_VVScore and Z_VH-score value and judging the category of the point, namely the point is a water body or a non-water body, and generating a scatter diagram;

iii, analyzing the scatter diagram to find Z for distinguishing water body from land_VVScore value and Z_VH-a range of score values;

iv, reacting Z_VVScore value and Z_VH-refinement of the range of score values and thinning of the Z_VVScore and Z_VHCombining a plurality of groups of different values of score, calculating a contact precision verification image, calculating the overall classification precision of all combinations, and finding the combination which enables the overall classification precision to be maximum as Z extracted in the current water area range_VVScore and Z_VH-score threshold.