A spatiotemporal retrospective analysis method of river and lake water encroachment based on remote sensing

Abstract

The invention discloses a spatiotemporal retrospective analysis method for river and lake water area encroachment based on remote sensing. Firstly, medium-resolution dense time series remote sensing image data and high-resolution remote sensing data of river and lake water area encroachment are collected respectively, and processed to obtain medium-resolution satellite remote sensing data set and high-resolution satellite remote sensing data set with consistent time series. Then, according to the high-resolution satellite remote sensing data set, various types of water body indexes and vegetation indexes are constructed to form characteristic bands of river and lake water bodies. River and lake water area and non-water area samples are selected, and an improved extreme gradient boosting algorithm is used to intelligently extract river and lake water areas to obtain the change range of river and lake water areas. Finally, according to the change range of river and lake water areas and the medium-resolution satellite remote sensing data set, a spatiotemporal retrospective evidence collection indicator system for river and lake water area encroachment is constructed, and quantitative extraction of long-time series retrospective evidence collection indicators for river and lake water area encroachment and key time point determination are carried out to provide support for cross-regional and long-time series spatiotemporal retrospective evidence collection for river and lake water area encroachment.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing image intelligent recognition, and specifically relates to the design of a spatiotemporal retrospective analysis method for river and lake water area encroachment based on remote sensing. Background Art

[0002] Rivers and lakes not only have economic and social functions such as water supply, aquaculture, transportation, shipping, sand mining, landscape, and entertainment, but also have ecosystem service functions such as providing biological habitats, maintaining biodiversity, absorbing pollutants, and regulating climate. They play a very important role in economic and social development and the construction of ecological civilization. The water area of ​​rivers and lakes is the core area of ​​river and lake supervision. Once encroached upon, it will have a serious impact on river flooding and ecological security. The encroachment of river and lake waters is mainly caused by the unreasonable use of river space by human society. It can be roughly divided into four categories: illegal occupation, illegal mining, illegal dumping, and illegal construction. For example, beach utilization within the water area, illegal sand mining or soil extraction, dumping, landfilling or stacking of solid waste, and construction of buildings and structures that hinder flood flow.

[0003] The encroachment of river and lake waters is large in scope, diverse in types and relatively hidden. Manual inspections require a lot of manpower and material resources, and the timeliness of inspections is low. Satellite remote sensing technology has the advantages of being macroscopic, accurate and real, and is an important means of conducting large-scale river and lake water surveys. With the rapid development of domestic high-resolution remote sensing technology, multiple high-resolution satellites such as Gaofen-1 to Gaofen-8 have been launched successively, forming a high-resolution multi-source satellite collaborative observation system, which provides rich and high-quality image data for large-scale and detailed retrospective analysis of river and lake water encroachment.

[0004] At present, relevant work has been carried out to monitor the encroachment of river and lake waters. Most of them use high-resolution satellite remote sensing data, and adopt time series comparison and manual visual interpretation to find the encroachment area of ​​river and lake waters. This method has the following shortcomings: (1) The manual visual interpretation method is inefficient, and the supervision of encroachment of river and lake waters in a large area is time-consuming and laborious. The timeliness of remote sensing monitoring of encroachment of river and lake waters is low. (2) There is a lack of spatiotemporal retrospective analysis of encroachment of river and lake waters, which makes it difficult to display the information of the whole process of encroachment of river and lake waters. The lack of evidence information makes it difficult to meet the needs of precise supervision. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that the existing remote sensing monitoring methods for river and lake water encroachment have low timeliness and lack of spatiotemporal retrospective analysis of river and lake water encroachment, making it difficult to meet the needs of precise supervision. A spatiotemporal retrospective analysis method for river and lake water encroachment based on remote sensing is proposed.

[0006] The technical solution of the present invention is: a spatiotemporal retrospective analysis method of river and lake water area occupation based on remote sensing, comprising the following steps:

[0007] S1. Collect medium-resolution dense time-series remote sensing image data of river and lake encroachment areas, and perform geometric and radiometric corrections to obtain a medium-resolution satellite remote sensing dataset with consistent time series.

[0008] S2. Collect high-resolution remote sensing data of river and lake encroachment areas and perform preprocessing to obtain a high-resolution satellite remote sensing data set with consistent time series.

[0009] S3. Intelligently extract river and lake water areas based on high-resolution satellite remote sensing data sets to obtain the range of changes in river and lake water areas.

[0010] S4. Conduct spatiotemporal retrospective forensic analysis of river and lake water encroachment based on the range of river and lake water changes and medium-resolution satellite remote sensing datasets.

[0011] Furthermore, the selection criteria for high-resolution remote sensing data of the river and lake encroachment area in step S2 are:

[0012] (1) Select clear sky images with no cloud coverage in river and lake areas;

[0013] (2) The image has clear surface layers and no obvious aerosol coverage;

[0014] (3) The image has no missing scan lines, stripes, noise, or abnormal pixels;

[0015] (4) The image is not covered with ice or snow.

[0016] Further, step S2 includes the following sub-steps:

[0017] S21. Collect 30m digital elevation model data and Sentinel-2 time-series high-resolution image data of river and lake encroachment areas, and use Sentinel-2 remote sensing data with a spatial resolution of 10m as reference image acquisition control points.

[0018] S22. Divide the time series high-resolution image into a grid of 10,000 × 10,000 pixels. In each grid, use the frequency phase matching method to collect image control points, remove control points with errors greater than 2 m, and retain 20 high-precision control points in each grid.

[0019] S23. Use high-precision control points to optimize the rational function coefficients of the panchromatic image in the time-series high-resolution image, build a high-precision geometric correction model to correct the panchromatic image, and achieve precise geometric correction of the multispectral image by aligning the panchromatic and multispectral images in the time-series high-resolution image.

[0020] S24. The Pansharp image fusion method is used to fuse the panchromatic and multispectral images in the time-series high-resolution images to obtain a high-resolution satellite remote sensing dataset with consistent time series.

[0021] Further, step S3 includes the following sub-steps:

[0022] S31. Divide the remote sensing images in the high-resolution satellite remote sensing data set into water and non-water areas, and use 10000×10000 pixel grids to spatially divide the remote sensing images. For each water and non-water area in each grid, randomly generate 500 sample points, and divide each sample point into water and non-water areas.

[0023] S32. Construct feature vectors of high-resolution remote sensing.

[0024] S33. Based on the characteristic vectors of high-resolution remote sensing and the water body or non-water body type of each sample point, the improved extreme gradient boosting algorithm is used to obtain the change range of river and lake water areas.

[0025] Furthermore, the feature vector V of the high-resolution remote sensing in step S32 ′ It is expressed as:

[0026] V ′ =V(B1,B2,B3,B4,B NDVI ,B NDWI ,B NDBSI ,B SLoPE )

[0027] Where V(·) represents the feature set, B1, B2, B3, B4 represent the pixel values ​​of the feature bands of band 1, band 2, band 3 and band 4 in the high-resolution remote sensing satellite image, respectively. NDVI Represents the pixel value of the normalized vegetation index characteristic band, B NDWI Represents the pixel value of the normalized water index characteristic band, B NDBSI Indicates the pixel value of the normalized bare soil index characteristic band, B SLOPE The pixel values ​​representing the slope index characteristic band.

[0028] The calculation formula of normalized difference vegetation index NDVI is:

[0029]

[0030] The calculation formula of normalized water index NDVI is:

[0031]

[0032] The calculation formula of normalized bare soil index NDBSI is:

[0033]

[0034] The calculation formula of slope index SLOPE is:

[0035]

[0036] where ρ nir Represents high-resolution remote sensing multispectral band 4, ρ red Represents high-resolution remote sensing multispectral band 3, ρ green represents high-resolution remote sensing multispectral band 2, ρ blue It represents high-resolution remote sensing multispectral band 1, dz / dx represents the ratio of the elevation z to the x direction, and dx / dy represents the ratio of the elevation z to the y direction.

[0037] Further, step S33 includes the following sub-steps:

[0038] S331, based on the feature vector V of high-resolution remote sensing ′ and the water or non-water type y of each sample point i , construct a set of K sample points {T1(V ′ ,y i ),T2(V ′ ,y i ),…,T K (V ′ ,y i )}, where y i =1 means the sample point is a water body, y i =0 means the sample point is not a water body, T K (V ′ ,y i ) represents the characteristics of the Kth sample point and the corresponding water body or non-water body type.

[0039] S332, by gradually adding regression trees to the extreme gradient boosting algorithm, the predicted value of water body or non-water body type for each sample point is obtained

[0040]

[0041] in represents the predicted value of the i-th water body or non-water body type generated for the t-th time, f k (·) represents the score of the kth regression tree, f t (·) represents the score of the tth regression tree.

[0042] S333, predicted value according to water body or non-water body type and the true value y i Construct the regularized objective function L(f t ):

[0043]

[0044] where l(·) represents the differential loss function, Ω(f t ) represents the regularization term, and the calculation formula is:

[0045]

[0046] Where γ and λ are constants that control the degree of regularization, T represents the number of leaves in the regression tree, and ω j Represents the weight of the j-th leaf node.

[0047] S334, calculate the regularization objective function L(f t )'s first-order statistic g i and the second-order statistic h i :

[0048]

[0049] in Express The partial derivative of .

[0050] S335, according to the first-order statistic g i and the second-order statistic h i For the regularized objective function L(f t ) is optimized to obtain the optimization objective function L (t) :

[0051]

[0052] Where n represents the number of water or non-water body types.

[0053] S336. Calculate the optimal solution of the optimization objective function

[0054]

[0055] Where q represents the sample, I j Represents a collection of water or non-water body types.

[0056] S337, the optimal solution As the leaf score of the regression tree, when the leaf score is less than the preset threshold, the regression tree stops growing the leaf depth and the classification results of water bodies and non-water bodies in the remote sensing image are obtained.

[0057] S338. Construct a river and lake water area set {A(T1,y a,b ),A(T2,y a,b ),…,A(T m ,y a,b )}, where A(T m ,y a,b ) indicates T m The type of river or lake water in the ath row and bth column of the remote sensing image at the moment, A(Tm ,y a,b )=1 means the river and lake water type is water body, A(T m ,y a,b )=0 means that the river or lake water type is non-water body.

[0058] S339. Obtain the river and lake water area change range A(a,b) according to the river and lake water area range set:

[0059]

[0060] Among them, A(a,b)=1 represents the area occupied by rivers and lakes, and A(a,b)=0 represents the area not occupied by rivers and lakes.

[0061] Further, step S4 includes the following sub-steps:

[0062] S41. Obtain retrospective evidence collection indicators for river and lake waters based on the river and lake waters management scope and river and lake waters encroachment area.

[0063] S42. Based on the scope of change of river and lake water areas and medium-resolution satellite remote sensing data sets, retrospective evidence collection is conducted on key time nodes of river and lake encroachment to obtain the duration of continuous encroachment of river and lake water areas.

[0064] S43. Based on the retrospective evidence collection indicators of river and lake waters and the continuous duration of river and lake waters encroachment, the spatiotemporal retrospective evidence collection information of river and lake waters encroachment is formed, and the spatiotemporal retrospective evidence collection analysis of river and lake waters encroachment is carried out.

[0065] Furthermore, the river and lake water area retrospective evidence indicators obtained in step S41 include the change area of ​​river and lake encroachment Length of natural shoreline or dike damage Changes in building or breeding area Interannual change rate of river and lake encroachment area Annual destruction length of natural coastline or levee and annual changes in construction or breeding areas

[0066] Changes in area occupied by rivers and lakes The calculation formula is:

[0067]

[0068] Where S T2 represents the area of ​​rivers and lakes at time T2, S T1 Represents the water area of ​​rivers and lakes at time point T1.

[0069] Length of natural shoreline or dike damage The calculation formula is:

[0070]

[0071] Where L T2 represents the length of natural shoreline or levee damage at time T2, L T1 Indicates the length of natural coastline or embankment damage at time point T1.

[0072] Changes in building or breeding area The calculation formula is:

[0073]

[0074] Where D T2 Indicates the building or breeding area at time T2, D T1 Indicates the building or breeding area at time point T1.

[0075] Interannual change rate of river and lake encroachment area The calculation formula is:

[0076]

[0077] Where S R2 represents the area of ​​rivers and lakes in year R2, S R1 It represents the water area of ​​rivers and lakes in the year R1, R1 represents the year when the encroachment of rivers and lakes begins, and R2 represents the year when the encroachment of rivers and lakes ends.

[0078] Annual destruction length of natural coastline or levee The calculation formula is:

[0079]

[0080] Where L R2 Indicates the length of natural shoreline or dike damage in year R2, L R1 Indicates the length of natural shoreline or embankment damage in year R1.

[0081] Annual changes in building or breeding area The calculation formula is:

[0082]

[0083] Where D R2 Indicates the building or breeding area in R2 year, D R1 Indicates the building or breeding area in year R1.

[0084] Further, step S42 includes the following sub-steps:

[0085] S421. According to the change range of river and lake water area A(i,j), determine the typical location of river and lake encroachment, and calculate the normalized water enhancement index eNDWI for each medium-resolution satellite remote sensing image of the typical location of river and lake encroachment:

[0086]

[0087] Where Green represents the reflectivity of the green light band, and MIR represents the reflectivity of the mid-infrared band.

[0088] S422. Draw a scatter plot with the date of the first medium-resolution satellite remote sensing image of the typical location of river and lake encroachment monitoring as the origin, the number of days from the date of each subsequent image to the first image as the x-axis, and the eNDWI value of the typical location of river and lake encroachment in each image as the y-axis, and construct a fitting function between x and y to fit the scatter plot:

[0089] y=Ax 3 +Bx 2 +Cx+D

[0090] Among them, A, B, C, and D are coefficients and constants obtained through polynomial fitting.

[0091] S423, determining the value of x when y=0 according to the fitting function, and selecting two adjacent image dates as the start time t1 and the end time t2 of the river and lake water area encroachment according to the value of x when y=0.

[0092] S424. Calculate the duration of river and lake water area occupation based on the start time t1 and the end time t2 of river and lake water area occupation:

[0093] Δt=t2-t1

[0094] Where Δt represents the duration of the occupation of river and lake waters.

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

[0096] (1) The present invention uses multi-source satellite remote sensing data to construct a remote sensing spatiotemporal retrospective analysis method for river and lake water encroachment, breaking through the technical bottleneck of cross-regional and long-term retrospective evidence collection of river and lake water encroachment, and has high timeliness.

[0097] (2) The present invention utilizes time-series high-resolution remote sensing data and medium-resolution remote sensing data, combined with an improved extreme gradient boosting algorithm, to perform intelligent extraction of long-time series river and lake waters, which has high accuracy and strong universality.

[0098] (3) The present invention constructs a remote sensing retrospective evidence collection indicator system for river and lake encroachment, quantitatively analyzes the spatiotemporal process of river and lake encroachment, and provides an objective basis for the supervision of river and lake encroachment.

[0099] (4) The present invention combines the dense time series of Landsat series and the medium-resolution dense time series remote sensing image data of Sentinel-2 to realize the determination of key time nodes of river and lake encroachment, providing a clear basis for the spatiotemporal evidence collection of river and lake encroachment. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 Shown is a flow chart of a spatiotemporal retrospective analysis method for river and lake water area encroachment based on remote sensing provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0101] Now, exemplary embodiments of the present invention will be described in detail with reference to the accompanying drawings. It should be understood that the embodiments shown and described in the accompanying drawings are only exemplary and are intended to explain the principles and spirit of the present invention, rather than to limit the scope of the present invention.

[0102] The embodiment of the present invention provides a spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing, such as Figure 1 As shown, the following steps S1 to S4 are included:

[0103] S1. Collect medium-resolution dense time-series remote sensing image data of river and lake encroachment areas, and perform geometric and radiometric corrections to obtain a medium-resolution satellite remote sensing dataset with consistent time series.

[0104] In the embodiment of the present invention, in order to improve the determination accuracy of key time nodes of river and lake encroachment, long-time series Landsat5\7\8 and Sentinel-2 series remote sensing data of river and lake encroachment areas are obtained as medium-resolution dense time-series remote sensing image data.

[0105] S2. Collect high-resolution remote sensing data of river and lake encroachment areas and perform preprocessing to obtain a high-resolution satellite remote sensing data set with consistent time series.

[0106] In the embodiment of the present invention, in order to realize the accurate extraction of retrospective evidence indicators of encroachment on river and lake waters, appropriate high-resolution remote sensing data with a spatial resolution better than 2m is selected, and the data selection criteria are as follows:

[0107] (1) Select clear sky images with no cloud coverage in river and lake areas;

[0108] (2) The image has clear surface layers and no obvious aerosol coverage;

[0109] (3) The image has no missing scan lines, stripes, noise, or abnormal pixels;

[0110] (4) The image is not covered with ice or snow.

[0111] Step S2 includes the following sub-steps S21 to S24:

[0112] S21. Collect 30m digital elevation model (DEM) data and Sentinel-2 time-series high-resolution image data of river and lake encroachment areas, and use Sentinel-2 remote sensing data with a spatial resolution of 10m as reference image acquisition control points.

[0113] S22. Divide the time series high-resolution image into a grid of 10,000 × 10,000 pixels. In each grid, use the frequency phase matching method to collect image control points, remove control points with errors greater than 2 m, and retain 20 high-precision control points in each grid.

[0114] In the embodiment of the present invention, considering that the imaging width of high-resolution remote sensing data is large, a large number of redundant control points will be collected by collecting control points through panoramic images. Therefore, a control point collection method based on grid division is adopted to improve the collection efficiency and correction accuracy.

[0115] S23. Use high-precision control points to optimize the rational function coefficients of the panchromatic image in the time-series high-resolution image, build a high-precision geometric correction model to correct the panchromatic image, and achieve precise geometric correction of the multispectral image by aligning the panchromatic and multispectral images in the time-series high-resolution image.

[0116] S24. The Pansharp image fusion method is used to fuse the panchromatic and multispectral images in the time-series high-resolution images to obtain a high-resolution satellite remote sensing dataset with consistent time series. This method has a good effect in preserving image information, details and spectrum.

[0117] S3. Intelligently extract river and lake water areas based on high-resolution satellite remote sensing data sets to obtain the range of changes in river and lake water areas.

[0118] Step S3 includes the following sub-steps S31 to S33:

[0119] S31. Divide the remote sensing images in the high-resolution satellite remote sensing data set into water and non-water areas, and use 10000×10000 pixel grids to spatially divide the remote sensing images. For each water and non-water area in each grid, randomly generate 500 sample points, and divide each sample point into water and non-water areas.

[0120] In the embodiment of the present invention, Open Street Map is used as reference data for high-resolution satellite remote sensing data, and small water body and non-water body samples can be selected to ensure the richness of the samples. The remote sensing image is divided into water body and non-water body areas, where water bodies include rivers, lakes, rivers, pools, near-shore seawater, silt water bodies and aquaculture, and non-water body types include farmland, forests, cities, bare land, bushes, mountain shadows, cloud shadows and other features.

[0121] S32. Construct feature vectors of high-resolution remote sensing.

[0122] In the embodiment of the present invention, the feature vector V of the high-resolution remote sensing ′ It is expressed as:

[0123] V ′=V(B1,B2,B3,B4,B NDVI ,B NDWI ,B NDBSI ,B SKOPE )

[0124] Where V(·) represents the feature set, B1, B2, B3, B4 represent the pixel values ​​of the feature bands of band 1, band 2, band 3 and band 4 in the high-resolution remote sensing satellite image, respectively. NDVI Represents the pixel value of the normalized vegetation index characteristic band, B NDWI Represents the pixel value of the normalized water index characteristic band, B NDBSI Indicates the pixel value of the normalized bare soil index characteristic band, B SLOPE Indicates the pixel value of the slope index characteristic band. In the embodiment of the present invention, a total of 8 characteristic bands are selected, including 4 image characteristic bands and 4 index characteristic bands.

[0125] The calculation formula of normalized difference vegetation index NDVI is:

[0126]

[0127] The calculation formula of normalized water index NDVI is:

[0128]

[0129] The calculation formula of normalized bare soil index NDBSI is:

[0130]

[0131] The calculation formula of slope index SLOPE is:

[0132]

[0133] where ρ nir Represents high-resolution remote sensing multispectral band 4, ρ red Represents high-resolution remote sensing multispectral band 3, ρ green represents high-resolution remote sensing multispectral band 2, ρ blue It represents high-resolution remote sensing multispectral band 1, dz / dx represents the ratio of the elevation z to the x direction, and dz / dy represents the ratio of the elevation z to the y direction.

[0134] S33. Based on the characteristic vectors of high-resolution remote sensing and the water body or non-water body type of each sample point, the improved extreme gradient boosting algorithm is used to obtain the change range of river and lake water areas.

[0135] Step S33 includes the following sub-steps S331 to S339:

[0136] S331, based on the feature vector V of high-resolution remote sensing′ and the water or non-water type y of each sample point i , construct a set of K sample points {T1(V ′ ,y i ),T2(V ′ ,y i ),…,T K (V ′ ,y i )}, where y i =1 means the sample point is a water body, y i =0 means the sample point is not a water body, T K (V ′ ,y i ) represents the characteristics of the Kth sample point and the corresponding water body or non-water body type.

[0137] S332, by gradually adding regression trees to the extreme gradient boosting algorithm, the predicted value of water body or non-water body type for each sample point is obtained

[0138]

[0139] in represents the predicted value of the i-th water body or non-water body type generated for the t-th time, f k (·) represents the score of the kth regression tree, f t (·) represents the score of the tth regression tree.

[0140] S333, predicted value according to water body or non-water body type and the true value y i Construct the regularized objective function L(f t ):

[0141]

[0142] where l(·) represents the differential loss function, Ω(f t ) represents the regularization term, and the calculation formula is:

[0143]

[0144] Where γ and λ are constants that control the degree of regularization, T represents the number of leaves in the regression tree, and ω j Represents the weight of the j-th leaf node.

[0145] S334, calculate the regularization objective function L(f t )'s first-order statistic g i and the second-order statistic h i :

[0146]

[0147] in Express The partial derivative of .

[0148] S335, according to the first-order statistic g i and the second-order statistic h i For the regularized objective function L(f t ) is optimized to obtain the optimization objective function L (t) :

[0149]

[0150] Where n represents the number of water or non-water body types.

[0151] S336. Calculate the optimal solution of the optimization objective function

[0152]

[0153] Where q represents the sample, I j Represents a collection of water or non-water body types.

[0154] S337, the optimal solution As the leaf score of the regression tree, when the leaf score is less than the preset threshold, the regression tree stops growing the leaf depth and the classification results of water bodies and non-water bodies in the remote sensing image are obtained.

[0155] S338. Construct a river and lake water area set {A(T1,y a,b ),A(T2,y a,b ),…,A(T m ,y a,b )}, where A(T m ,y a,b ) indicates T m The type of river or lake water in the ath row and bth column of the remote sensing image at the moment, A(T m ,y a,b )=1 means the river and lake water type is water body, A(T m ,y a,b )=0 means that the river or lake water type is non-water body.

[0156] S339. Obtain the river and lake water area change range A(a,b) according to the river and lake water area range set:

[0157]

[0158] Among them, A(a,b)=1 represents the area occupied by rivers and lakes, and A(a,b)=0 represents the area not occupied by rivers and lakes.

[0159] S4. Conduct spatiotemporal retrospective forensic analysis of river and lake water encroachment based on the range of river and lake water changes and medium-resolution satellite remote sensing datasets.

[0160] Step S4 includes the following sub-steps S41 to S43:

[0161] S41. Obtain retrospective evidence collection indicators for river and lake waters based on the river and lake waters management scope and river and lake waters encroachment area.

[0162] In the embodiment of the present invention, the river and lake water area retrospective evidence collection indicators include the change area of ​​river and lake encroachment Length of natural shoreline or dike damage Changes in building or breeding area Interannual change rate of river and lake encroachment area Annual destruction length of natural coastline or levee and annual changes in construction or breeding areas

[0163] Changes in area occupied by rivers and lakes The calculation formula is:

[0164]

[0165] Where S T2 represents the area of ​​rivers and lakes at time T2, S T1 Represents the water area of ​​rivers and lakes at time point T1.

[0166] Length of natural shoreline or dike damage The calculation formula is:

[0167]

[0168] Where L T2 represents the length of natural shoreline or levee damage at time T2, L T1 Indicates the length of natural coastline or embankment damage at time point T1.

[0169] Changes in building or breeding area The calculation formula is:

[0170]

[0171] Where D T2 Indicates the building or breeding area at time T2, D T1 Indicates the building or breeding area at time point T1.

[0172] Interannual change rate of river and lake encroachment area The calculation formula is:

[0173]

[0174] Where S R2 represents the area of ​​rivers and lakes in year R2, S R1 It represents the water area of ​​rivers and lakes in the year R1, R1 represents the year when the encroachment of rivers and lakes begins, and R2 represents the year when the encroachment of rivers and lakes ends.

[0175] Annual destruction length of natural coastline or levee The calculation formula is:

[0176]

[0177] Where L R2 Indicates the length of natural shoreline or dike damage in year R2, L R1 Indicates the length of natural shoreline or embankment damage in year R1.

[0178] Annual changes in building or breeding area The calculation formula is:

[0179]

[0180] Where D R2 Indicates the building or breeding area in R2 year, D R1 Indicates the building or breeding area in year R1.

[0181] S42. Based on the scope of change of river and lake water areas and medium-resolution satellite remote sensing data sets, retrospective evidence collection is conducted on key time nodes of river and lake encroachment to obtain the duration of continuous encroachment of river and lake water areas.

[0182] In the retrospective analysis of river and lake encroachment, it is extremely important to determine the key time nodes of river and lake encroachment. Since the high-resolution satellite remote sensing data has a long recurrence period, the embodiment of the present invention uses a medium-resolution satellite remote sensing data set, introduces the normalized water enhancement index eNDWI, and conducts intensive time series detection of rivers and lakes in typical areas to determine the key time nodes of river and lake encroachment.

[0183] Step S42 includes the following sub-steps S421 to S424:

[0184] S421. According to the change range of river and lake water area A(i,j), determine the typical location of river and lake encroachment, and calculate the normalized water enhancement index eNDWI for each medium-resolution satellite remote sensing image of the typical location of river and lake encroachment:

[0185]

[0186] Where Green represents the reflectivity of the green light band, and MIR represents the reflectivity of the mid-infrared band.

[0187] According to the calculated typical location eNDWI, the eNDWI index change graph is obtained. Since the value range of the eNDWI index is -1 to 1, the closer the eNDWI is to 1, the greater the probability of water body, and the closer the eNDWI is to -1, the greater the probability of non-water body. Therefore, the embodiment of the present invention uses eNDWI = 0 as the judgment threshold for changes in river and lake water area occupation, and the occupation time period is the time interval of each image before and after eNDWI = 0.

[0188] S422. Draw a scatter plot with the date of the first medium-resolution satellite remote sensing image of the typical location of river and lake encroachment monitoring as the origin, the number of days from the date of each subsequent image to the first image as the x-axis, and the eNDWI value of the typical location of river and lake encroachment in each image as the y-axis, and construct a fitting function between x and y to fit the scatter plot:

[0189] y=Ax 3 +Bx 2 +Cx+D

[0190] Among them, A, B, C, and D are coefficients and constants obtained through polynomial fitting.

[0191] S423, determining the value of x when y=0 according to the fitting function, and selecting two adjacent image dates as the start time t1 and the end time t2 of the river and lake water area encroachment according to the value of x when y=0.

[0192] S424. Calculate the duration of river and lake water area occupation based on the start time t1 and the end time t2 of river and lake water area occupation:

[0193] Δt=t2-t1

[0194] Where Δt represents the duration of the occupation of river and lake waters.

[0195] S43. Based on the retrospective evidence collection indicators of river and lake waters and the continuous time of encroachment on river and lake waters, the spatiotemporal retrospective evidence collection information on the encroachment on river and lake waters is formed, and the spatiotemporal retrospective evidence collection analysis on the encroachment on river and lake waters is carried out to provide timely and accurate information for the supervision and management of encroachment on river and lake waters.

[0196] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific 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 spatiotemporal retrospective analysis method for river and lake water encroachment based on remote sensing, characterized in that: The following steps are involved: S1. Collect medium-resolution dense time-series remote sensing image data of river and lake encroachment areas, and perform geometric and radiometric corrections to obtain a medium-resolution satellite remote sensing dataset with consistent time series. S2, collect high-resolution remote sensing data of river and lake encroachment areas and perform preprocessing to obtain a high-resolution satellite remote sensing data set with consistent time series; S3, intelligently extract river and lake water areas based on high-resolution satellite remote sensing data sets to obtain the range of river and lake water area changes; S4. Conduct temporal and spatial retrospective forensic analysis of river and lake water encroachment based on the range of river and lake water changes and medium-resolution satellite remote sensing data sets; The step S4 comprises the following sub-steps: S41. Obtain river and lake waters retrospective evidence collection indicators based on the river and lake waters management scope and river and lake waters encroachment area; S42. Based on the scope of change of river and lake waters and medium-resolution satellite remote sensing data sets, retrospectively collect evidence on key time nodes of river and lake encroachment to obtain the duration of continuous encroachment of river and lake waters; S43. Based on the retrospective evidence collection indicators of river and lake waters and the duration of continuous occupation of river and lake waters, form the temporal and spatial retrospective evidence collection information of river and lake waters occupation, and conduct temporal and spatial retrospective evidence collection analysis on river and lake waters occupation; The step S42 includes the following sub-steps: S421. According to the change range of river and lake water area A(a,b), determine the typical location of river and lake encroachment, and calculate the normalized water enhancement index eNDWI for each medium-resolution satellite remote sensing image of the typical location of river and lake encroachment: Where Green represents the reflectivity of the green light band, and MIR represents the reflectivity of the mid-infrared band; S422. Draw a scatter plot with the date of the first medium-resolution satellite remote sensing image of the typical location of river and lake encroachment monitoring as the origin, the number of days from the date of each subsequent image to the first image as the x-axis, and the eNDWI value of the typical location of river and lake encroachment in each image as the y-axis, and construct a fitting function between x and y to fit the scatter plot: y=Ax 3 +Bx 2 +Cx+D Among them, A, B, C, and D are coefficients and constants obtained through polynomial fitting; S423, determining the value of x when y=0 according to the fitting function, and selecting two adjacent image dates as the start time t1 of the river and lake water area encroachment and the end time t2 of the river and lake water area encroachment according to the value of x when y=0; S424. Calculate the duration of river and lake water area occupation based on the start time t1 and the end time t2 of river and lake water area occupation: Δt=t2-t1 Where Δt represents the duration of the encroachment on river and lake waters.

2. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 1 is characterized in that: The selection criteria of high-resolution remote sensing data of the river and lake encroachment area in step S2 are: (1) Select clear sky images with no cloud coverage in river and lake areas; (2) The image has clear surface layers and no obvious aerosol coverage; (3) The image has no missing scan lines, stripes, noise, or abnormal pixels; (4) The image is not covered with ice or snow.

3. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 1 is characterized in that: The step S2 comprises the following sub-steps: S21. Collect 30m digital elevation model data and Sentinel-2 time series high-resolution image data of river and lake encroachment areas, and use Sentinel-2 10m spatial resolution remote sensing data as reference image acquisition control points; S22, divide the time series high-resolution image into a grid of 10000×10000 pixels, collect image control points in each grid using the frequency phase matching method, remove control points with errors higher than 2m, and retain 20 high-precision control points in each grid; S23. Optimize the rational function coefficients of the panchromatic image in the time-series high-resolution image using high-precision control points, build a high-precision geometric correction model to correct the panchromatic image, and achieve precise geometric correction of the multispectral image by registering the panchromatic and multispectral images in the time-series high-resolution image; S24. The Pansharp image fusion method is used to fuse the panchromatic and multispectral images in the time-series high-resolution images to obtain a high-resolution satellite remote sensing dataset with consistent time series.

4. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 1 is characterized in that: The step S3 comprises the following sub-steps: S31. Divide the remote sensing images in the high-resolution satellite remote sensing data set into water and non-water areas, and use 10000×10000 pixel grids to spatially divide the remote sensing images. For the water and non-water areas in each grid, randomly generate 500 sample points, and divide each sample point into water and non-water areas. S32, construct feature vectors of high-resolution remote sensing; S33. Based on the characteristic vectors of high-resolution remote sensing and the water body or non-water body type of each sample point, the improved extreme gradient boosting algorithm is used to obtain the change range of river and lake water areas.

5. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 4 is characterized in that: The characteristic vector V′ of the high-resolution remote sensing in step S32 is expressed as: V′=V(B1,B2,B3,B4,B NDVI ,B NDWI ,B NDBSI ,B SLOPE ) Where V(·) represents the feature set, B1, B2, B3, B4 represent the pixel values ​​of the feature bands of band 1, band 2, band 3 and band 4 in the high-resolution remote sensing satellite image, respectively. NDVI Represents the pixel value of the normalized vegetation index characteristic band, B NDWI Represents the pixel value of the normalized water index characteristic band, B NDBSI Indicates the pixel value of the normalized bare soil index characteristic band, B SLOPE The pixel value representing the characteristic band of slope index; The calculation formula of the normalized difference vegetation index NDWI is: The calculation formula of the normalized water index NDWI is: The calculation formula of the normalized bare soil index NDBSI is: The calculation formula of the slope index SLOPE is: where ρ nir Represents high-resolution remote sensing multispectral band 4, ρ red Represents high-resolution remote sensing multispectral band 3, ρ green represents high-resolution remote sensing multispectral band 2, ρ blue It represents high-resolution remote sensing multispectral band 1, dz / dx represents the ratio of the elevation z to the x direction, and dz / dy represents the ratio of the elevation z to the y direction.

6. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 4 is characterized in that: The step S33 comprises the following sub-steps: S331, based on the feature vector V′ of high-resolution remote sensing and the water body or non-water body type y of each sample point i , construct a set of K sample points {T1(V′,y i ),T2(V′,y i ),…,T K (V′,y i )}, where y i =1 means the sample point is a water body, y i =0 means the sample point is not a water body, T K (V′,y i ) represents the characteristics of the Kth sample point and the corresponding water body or non-water body type; S332, by gradually adding regression trees to the extreme gradient boosting algorithm, the predicted value of water body or non-water body type for each sample point is obtained in represents the predicted value of the i-th water body or non-water body type generated for the t-th time, f k (·) represents the score of the kth regression tree, f t (·) represents the score of the tth regression tree; S333, predicted value according to water body or non-water body type and the true value y i Construct the regularized objective function L(f t ): where l(·) represents the differential loss function, Ω(f t ) represents the regularization term, and the calculation formula is: Where γ and λ are constants that control the degree of regularization, T represents the number of leaves in the regression tree, and ω j Represents the weight of the jth leaf node; S334, calculate the regularization objective function L(f t )'s first-order statistic g i and the second-order statistic h i : in Express The partial derivative of S335, according to the first-order statistic g i and the second-order statistic h i For the regularized objective function L(f t ) is optimized to obtain the optimization objective function L (t) : Where n represents the number of water or non-water body types; S336. Calculate the optimal solution of the optimization objective function Where q represents the sample, I j Represents a collection of water or non-water body types; S337, the optimal solution As the leaf score of the regression tree, when the leaf score is less than the preset threshold, the regression tree stops growing the leaf depth, and the classification results of water bodies and non-water bodies in the remote sensing image are obtained; S338. Construct a river and lake water area set {A(T1,y a,b ),A(T2,y a,b ),…,A(T m ,y a,b )}, where A(T m ,y a,b ) indicates T m The type of river or lake water in the ath row and bth column of the remote sensing image at the moment, A(T m ,y a,b )=1 means the river and lake water type is water body, A(T m ,y a,b )=0 means the river and lake water type is non-water body; S339. Obtain the river and lake water area change range A(a,b) according to the river and lake water area range set: Among them, A(a,b)=1 represents the area occupied by rivers and lakes, and A(a,b)=0 represents the area not occupied by rivers and lakes.

7. The spatiotemporal retrospective analysis method of river and lake water area encroachment based on remote sensing according to claim 1 is characterized in that: The river and lake water area retrospective evidence indicators obtained in step S41 include the change area of ​​river and lake encroachment Length of natural shoreline or dike damage Changes in building or breeding area Interannual change rate of river and lake encroachment area Annual destruction length of natural coastline or levee and annual changes in construction or breeding areas Changes in the area occupied by rivers and lakes The calculation formula is: Where S T2 represents the area of ​​rivers and lakes at time T2, S T1 It represents the water area of ​​rivers and lakes at time point T1; The length of natural coastline or levee damage The calculation formula is: Where L T2 represents the length of natural shoreline or levee damage at time T2, L T1 It indicates the length of natural shoreline or dike damage at time T1; Changes in building or breeding area The calculation formula is: Where D T2 Indicates the building or breeding area at time T2, D T1 Indicates the building or breeding area at time point T1; The interannual change rate of the area occupied by rivers and lakes The calculation formula is: Where S R2 represents the area of ​​rivers and lakes in year R2, S R1 represents the water area of ​​rivers and lakes in the year R1, R1 represents the year when the encroachment of rivers and lakes begins, and R2 represents the year when the encroachment of rivers and lakes ends; The annual destruction length of the natural coastline or levee The calculation formula is: Where L R2 Indicates the length of natural shoreline or dike damage in year R2, L R1 It indicates the length of natural shoreline or dike damage in year R1; The annual changes in the building or breeding area The calculation formula is: Where D R2 Indicates the building or breeding area in R2 year, D R1 Indicates the building or breeding area in year R1.

Images