Method and device for determining suitable vegetation growth area of inland river shore belt based on multi-source images

Abstract

The application provides a method and device for determining a suitable growth area of vegetation in an inland river shore belt based on multi-source images, and the method comprises the following steps: 1. Based on image data of a research river section, water and land area division is performed on each image, and a binary submerged image data set is constructed; 2. Shore beach long-duration average submergence frequency calculation based on pixel scale; 3. Shore vegetation coverage area is extracted from each image, and a binary coverage image set is constructed; 4. Long-duration inland shore beach vegetation appearance frequency calculation based on pixel scale; 5. The target river section area is divided into suitable growth areas of the inland shore belt vegetation; according to the grid data of the multi-year average submergence frequency FF and the vegetation appearance frequency VF of the target river section area, a scatter plot and a kernel density cloud map of FF~VF are obtained, and then the suitable growth area of the shore vegetation is divided, and the critical submergence frequency threshold of the suitable growth of the shore vegetation is obtained.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of ecological river management, and particularly relates to a method and device for determining a suitable area for vegetation growth in a bank zone of an inland river based on multi-source images. BACKGROUND

[0002] Bank zone greening is an important part of ecological river management. For example, the high bank line of the two banks of a city river section is often developed into a leisure landscape area covered with a large amount of vegetation, and ecological revetment and bank protection technology is often used to replace the traditional bank protection mode in river and channel regulation projects. In the design and construction of these projects, the bank zone vegetation configuration needs to be adjusted to local conditions, and the existing river topography, hydrology and other conditions need to be fully utilized and supplemented with appropriate artificial restoration to build a stable habitat condition suitable for plant growth. After the implementation of many ecological bank protection projects, a large area of vegetation has died, and some projects have even adopted the mode of "construction-construction effect feedback-modification and remediation", which has caused waste of resources and led to these phenomena. The reason is that the planning and design stage cannot clearly determine which areas of the bank zone are suitable for vegetation growth and what hydrological conditions the vegetation growth area needs to meet.

[0003] Unlike the simple case of a fixed water level in an artificial channel and a fixed location of a beach in a tidal zone, the topography of a natural river bank is complex and is subject to a large fluctuation in water level within a year and from year to year. The distribution of vegetation in the bank zone presents both a general characteristic of a step distribution along the river cross section and a strong spatial and temporal randomness, and it is difficult to determine the suitable area for vegetation growth by a quantitative method. In the past, some engineering construction has used a method combining field investigation and water power numerical simulation to determine the suitable area for vegetation growth, but this process requires long-term investigation data, has a high modeling difficulty and cost, and generally does not have the implementation conditions. SUMMARY

[0004] The present application is made to solve the above problems, and aims to provide a method and device for determining a suitable area for vegetation growth in a bank zone of an inland river based on multi-source images, which has a small amount of requirements for measured topographic data, strong data complementarity, a simple and quantitative operation process, fully considers the hydrological characteristics of a natural river with a large fluctuation in water level and strong randomness and the complex morphological characteristics of a bank with a high and low relief and a water-land interlaced complex form, and more accurate calculation of the submergence frequency, so as to provide a more scientific, efficient and reliable implementation approach for ecological river management.

[0005] To achieve the above object, the present application adopts the following scheme:

[0006] <Method>

[0007] The present application provides a method for determining a suitable area for vegetation growth in a bank zone of an inland river based on multi-source images, characterized in that it comprises the following steps:

[0008] Step 1. Based on the image data of the study reach, water-land area division is implemented on each image, and a binary inundation image dataset is constructed;

[0009] Step 2. The average inundation frequency of the beach length is calculated based on the pixel scale; specifically including the following sub-steps:

[0010] Step 2.1: According to whether the target river reach has long series water level observation data, select the calculation method of the beach inundation frequency, if it has water level observation data, go to step 2.2, otherwise go to step 2.10;

[0011] Step 2.2: Date as index, match the daily average water level of the date corresponding to each binary inundation image {B 11 ,.......,B ij}, recorded as {Z 11 ,.......,Z ij}, where Z ij represents the water level corresponding to the i-th image on the j-th date in the i-th year;

[0012] Step 2.3: According to the statistical requirements of regional inundation frequency, sort the daily average water level in T years from small to large, and divide it into m water level interval intervals at equal intervals, 25≤m≤35;

[0013] Step 2.4: Statistic the frequency F of water level {Z 11 ,.......,Z ij} falling in m water level interval, keep the water level interval with frequency F≥1, merge the water level interval with frequency F=0 and the adjacent interval, the number of water level interval after merging is recorded as m', 20≤m';

[0014] Step 2.5: According to the daily average water level data year by year, statistic the number of days N ik (k=1,2,……,m') that the daily water level falls in m' water level interval in the i-th year (i=1,2,……,T), the number of days N ik is the time length of the k-th water level interval in the i-th year;

[0015] Step 2.6: According to the water level data {Z 11 ,.......,Z ij} corresponding to the time sequence image, statistic the number of images C that the water level Z ij of the j-th image in the i-th year falls in m' water level interval;

[0016] Step 2.7: According to step 2.5, match the representative image of each water level interval in the i-th year, coded as B ik , B ik represents the binary inundation image of the k-th water level interval in the i-th year;

[0017] Step 2.8: for all binary images B ik , the formula is operated on a pixel-by-pixel basis to obtain the submergence frequency of each pixel in the i-th year, and the submergence frequency values of all pixels constitute the grid data FF i of the spatial distribution of the submergence frequency of the target river reach region in the i-th year.

[0018] Step 2.9: on the basis of the submergence frequency of each year, the formula is operated on a pixel-by-pixel basis to obtain the spatial distribution grid data FF of the multi-year average submergence frequency within T years.

[0019] Step 2.10: if there is no water level observation data, the time coordinates of the binary image set {B 11 ,.......,B ij} are used as the basis for calculating the submergence duration, the time midpoint of the corresponding dates of the j-1th and jth images in the i-th year is denoted as t s , the time midpoint of the corresponding dates of the jth and j+1th images is denoted as t e , and N ij =t e -t s is taken as the representative period of the jth image, and the formula is operated on a pixel-by-pixel basis to obtain the submergence frequency of each pixel, which together form the spatial distribution grid data FF i ' of the submergence frequency in the i-th year, j=1,2,……,n i , and n i is the number of images collected in the i-th year; the calculation method of the spatial distribution grid data of the multi-year average submergence frequency is the same as step 2.8.

[0020] Step 3: the vegetation coverage area of the riparian zone is extracted from the images one by one, and a binary coverage image set is constructed.

[0021] Step 4: the vegetation appearance frequency of the riparian beach in a long duration is calculated on a pixel-by-pixel basis.

[0022] Step 5: the riparian zone in the target river reach region is divided into suitable zones; the scatter plot and kernel density cloud plot of FF~VF are obtained according to the grid data of the multi-year average submergence frequency FF and the vegetation appearance frequency VF of the target river reach region, and then the suitable zone for the growth of the riparian vegetation is divided, and the critical submergence frequency threshold for the suitable growth of the riparian vegetation is derived.

[0023] Preferably, the method for determining the suitable zone for the growth of the riparian vegetation in the inland river provided by the present application, in step 1, the binary submergence image data set {MB 11 ,.......,MB ij} and {SB11 ,.......,SB ij} into 30m pixel scale, then merge {MB 11 ,.......,MB ij} and {SB 11 ,.......,SB ij} into 30m precision binary inundation image dataset {B 11 ,.......,B ij}, where B ij represents the jth image of the ith year. Considering the spatial and temporal limitations of single image dataset in monitoring regional water body vegetation, synthetic aperture radar and other image datasets are introduced to make up for the problems such as insufficient data caused by cloud and fog interference of multispectral image, and the expansion of data is realized by using multisource remote sensing image dataset through differential water body index selection and pixel scale conversion, which provides sufficient data quantity for subsequent fine calculation of inundation frequency and vegetation appearance frequency.

[0024] Preferably, the method for determining the suitable vegetation growth area of the inland river shore belt based on multi-source images provided by the application comprises the following steps:

[0025] Step 1.1: Obtain image data of the target river section for T years from a data center (for example, Google Earth Engine), T≥10, wherein multispectral image data subjected to geometric correction, radiation calibration and atmospheric correction and other preprocessing methods are mainly used, such as Landsat SR secondary data products (30m precision), Sentinel-2 level 2A data (10m precision) and the like, and synthetic aperture radar data subjected to radiation positioning, correction and denoising are used as auxiliary, such as Sentinel-1 SAR GRD ground distance image data (10m precision). And filter multispectral image data that meet the quality and quantity conditions: (1) in the remote sensing image, the target river section area should meet the condition that no more than 5% of the cloud cover is blocked, and there is no damage and strip; (2) within each year, there should be available images in each season;

[0026] Step 1.2: use the QA band to mask the cloud layer for Landsat SR secondary data, and use the SLC band to mask the cloud layer for Sentinel-2 level 2A data; smooth and denoise the SAR GRD data; after removing the cloud, crop the target river section to form a time-series multispectral image dataset {MS 11 ,.......,MS ij}, and a time-series synthetic aperture radar dataset {SAR 11 ,.......,SAR ij}, where i represents the i-th year, j represents the j-th image in the i-th year, i = 1, 2, ..., T; j = 1, 2, ..., n i n i Let be the number of images collected in year i;

[0027] Step 1.3: Select {MS 11 MS ij Each image in MS ij ρ GREEN Green band and ρ SWIR The reflectivity data for the shortwave infrared band is calculated using the formula MNDWI = (ρ GREEN -ρ SWIR ) / (ρ GREEN +ρ SWIR Band operations were performed to obtain the improved Normalized Difference Water Index (MNDWI) image, thereby enabling the {MS} to be analyzed. 11 MS ij The dataset was transformed into an improved normalized differential water index MNDWI image dataset {M 11 ,.......,M ij}, where M ij This represents the j-th MNDWI image from year i; for the synthetic aperture radar dataset {SAR} 11 ,.......,SAR ij}, select each image SAR ij The dual-polarized data (VV and VH bands) from the Level 1 ground distance detection product after projection correction are used to calculate the Water Extraction Index (SDWI) according to the formula SDWI = ln(10 × VV × VH) - 8, thus forming the SDWI image dataset {S 11 ,.......,S ij}, where S ij This represents the j-th SDWI image of year i.

[0028] Step 1.4: Based on Step 1.3, for the image dataset {M} 11 ,.......,M ij} and {S 11 ,.......,S ij Output the grayscale histogram of each image in the target river section area. The horizontal axis represents the grayscale value (MNDWI or SDWI value), and the vertical axis represents the pixel frequency H under that grayscale value. The grayscale histogram should be bimodal.

[0029] Step 1.5: Based on the bimodal distribution characteristics of the gray histogram of each image, the OTSU maximum between-class variance method is used to determine the segmentation threshold of water and non-water body for each image. The implementation process is as follows: (1) classify the gray histogram of each image, first assume a classified gray threshold, then according to the formula 2 = ω0× ω1× (μ0- μ1) 2 Calculate the between-class variance σ 2 of the gray threshold, record the proportion of pixels less than the gray value in the map as ω0, and the average gray value as μ0; record the proportion of pixels greater than the gray value as ω1, and the average gray value as μ1; (2) constantly re- assume the classified gray threshold, repeat step (1) until the between-class variance σ 2 takes the maximum value, and the gray threshold threshold at this time is the segmentation threshold of water and non-water body.

[0030] Step 1.6: After calculating the water and non-water body segmentation threshold threshold for each image, the data sets {M 11 ,.......,M ij} and {S 11 ,.......,S ij} are binarized according to the segmentation threshold of each image, the water area pixels are assigned as 1, and the non-water body pixel area is assigned as 0, to obtain the binarized inundation data sets {MB 11 ,.......,MB ij} and {SB 11 ,.......,SB ij}, wherein MB ij represents the i-th year j-th binarized inundation image based on Landsat SR and Sentinel-2 2A data; SB ij represents the i-th year j-th binarized inundation image based on the synthetic aperture radar data set Sentinel-1 SAR data;

[0031] Step 1.7: For the 10m precision binarized images generated based on Sentinel-2 2A / Sentinel-1 SAR satellite in the binarized inundation image data sets {MB 11 ,.......,MB ij} and {SB 11 ,.......,SB ij}, the pixel scale is converted to 30m by using the resampling technology, and then {MB 11 ,.......,MB ij} and {SB 11 ,.......,SB ij} are merged into 30m precision binarized inundation image data set {B}.11 ,.......,B ij}, where B ij This represents the j-th image in the i-th year.

[0032] Preferably, in the method for determining suitable vegetation growth zones along inland riverbanks based on multi-source imagery provided by the present invention, in step 2.6, if the number of images C ≥ 2 in the k-th interval, then the image with the daily average water level closest to the median of that interval is selected as the representative image for that interval. For example, if the water level interval is divided into (32m~33m), then the image with the water level closest to 32.5m is selected as the representative image; if the number of images C = 0 in the k-th interval in a certain year, then the corresponding water level Z of the images in year T is traversed. ij Take the water level Z corresponding to the image. ij The image closest to the median of the water level interval and closest to the year is used as the replacement for the representative image of the k-th interval of that year.

[0033] Preferably, in the method for determining suitable vegetation growth zones along inland riverbanks based on multi-source imagery provided by the present invention, step 3 includes the following sub-steps:

[0034] Step 3.1: Based on the target river segment temporal multispectral image dataset {MS} obtained in Steps 1.1 and 1.2 11 MS ij}, take each image MS ij ρ RED Red light band and ρ NIR Infrared reflectance data is obtained using the formula NDVI = (ρ NIR -ρ RED ) / (ρ NIR +ρ RED Band operations were performed to obtain the Normalized Differential Vegetation Index (NDVI) imagery, and a long-term NDVI imagery dataset of the target river section was formed. 11 ,.......,V ij}, where V ij This represents the j-th NDVI image of year i;

[0035] Step 3.2: For the NDVI image dataset {V 11 ,.......,V ij Using a band calculator, pixel masks (NDVI≤0) are applied to the water area frame by frame. After masking, RNDVI images containing both light beaches and vegetation are generated, forming a dataset {RN}. 11 ,.......,RN ij};

[0036] Step 3.3: {RN 11 ,.......,RNij RN in} ij The image grayscale histogram shows a bimodal shape. The threshold' for segmenting the light beach and vegetation is calculated frame by frame using the OTSU maximum inter-class variance method in step 1.5.

[0037] Step 3.4: Divide the dataset {RN} according to the segmentation threshold ' for each image. 11 ,.......,RN ij Binarize each image individually, assigning a value of 1 to pixels in vegetated areas and a value of 0 to pixels in bare beach areas, to obtain the binarized cover dataset {MVB}. 11 ,.......,MVB ij}, where MVB ij This represents the j-th binary vegetation cover image of year i. By simplifying the classification of natural riverbed features, a relatively simple vegetation cover extraction method is used to achieve rapid identification and accurate extraction of the regional vegetation cover range.

[0038] Step 3.5: Perform pixel scale transformation on the binarized overlay image from Sentinel-2 2A 10m resolution, and use resampling technology to convert its pixels from 10m to 30m, forming a 30m resolution binarized overlay dataset {P}. 11 ,.......,P ij}

[0039] Preferably, in the method for determining suitable vegetation growth zones along inland riverbanks based on multi-source imagery provided by the present invention, step 4 includes the following sub-steps:

[0040] Step 4.1: Use the binarized overlay image set {P} 11 ,.......,P ij The time coordinates are used as the basis for calculating the duration of vegetation appearance. The midpoint of the time between the corresponding dates of the (j-1)th and jth images in the i-th year is denoted as t. s The midpoint of the time corresponding to the date of the j-th image and the (j+1)-th image is denoted as t. e Take N ij =t e -t s As the representative time period of the j-th image; based on this, according to the formula The vegetation occurrence frequency of each pixel is calculated to obtain the spatial distribution raster data VF of vegetation occurrence frequency in year i. i ;

[0041] Step 4.2: Based on the annual vegetation occurrence frequency, each pixel is processed according to the formula. The calculations are performed to obtain raster data of the spatial distribution of the multi-year average vegetation occurrence frequency within year T.

[0042] The vegetation appearance frequency of the pixel fully considers the multi-year change process of the riparian zone vegetation, and the vegetation appearance frequency reflects not only the regional information of the riparian zone vegetation growth affected by the periodic flooding, but also the species of the vegetation. The region with low vegetation appearance frequency is the new pioneer plant of the riparian zone, and the region with high vegetation appearance frequency represents the perennial plant that has been stably grown.

[0043] Preferably, the method for determining the suitable growth region of the riparian zone vegetation in the inland river based on multi-source images provided by the application comprises the following sub-steps in step 5.

[0044] Step 5.1: According to the grid data of the multi-year average flooding frequency FF and the vegetation appearance frequency VF of the target river region, the flooding frequency and the vegetation appearance frequency are output pixel by pixel and stored in two arrays; and a scatter plot is drawn based on the FF-VF relationship of all pixels in the two arrays, wherein FF is the horizontal coordinate and VF is the vertical coordinate.

[0045] Step 5.2: The clustering analysis is performed on the FF-VF scatter plot, and the generated scatter plot is converted into a kernel density cloud chart; specifically, the KS density function of MATLAB is used to perform the kernel density estimation on the scatter plot data, the data is smoothed by using the Gaussian kernel function, the estimated density is calculated, the estimated probability density is normalized, and then visualized, and the FF-VF kernel density estimation cloud chart is drawn. The value in the cloud chart represents the proportion of the pixel under the condition of a certain flooding frequency and a certain vegetation appearance frequency combination to the total pixels.

[0046] Step 5.3: The contour lines are added to the FF-VF kernel density cloud chart, and three partitions are formed according to the closed condition of the 2% to 5% contour lines. The projection of the boundary between the two partitions on the vertical axis corresponds to two vegetation appearance frequency thresholds CVF1 and CVF2, according to which the VF grid data is divided into three categories, corresponding to three levels of vegetation appearance frequency, wherein the pixel with VF≥CVF2 corresponds to the stable vegetation region, the region with CVF2>VF>CVF1 is the transition vegetation region, and the region with CVF1≥VF is the unsuitable vegetation growth region; the visualization of these vegetation partitions is performed in combination with the geographic spatial information, and the riparian zone vegetation suitability partition chart is made.

[0047] Step 5.4: For the FF-VF scatter plot, FF is divided into 100 equidistant intervals according to the step value of 1%, the mean value of FF in each interval is FFA, and the mean value of VF in each interval is VFA. The logistic S-shaped curve is fitted according to FFA and VFA by using the origin software, and the general form is y=a+(b / (1+(x / c) d), where a, b, c, d are constants, a is the longitudinal offset, b is the upper limit of the curve, c is the center position, and d is the slope of the curve;

[0048] Step 5.5: Corresponding the two vegetation frequency CVF1, CVF2 defined in step 5.2 to the logistic S-shaped curve, the critical submergence frequencies CFF1, CFF2 between the vegetation suitable area, transition area and unsuitable area are obtained. In engineering practice, the above critical values can be used to implement appropriate modification of the shore zone terrain. By adjusting the terrain elevation, the submergence frequency of a certain position is greater than or less than the above critical value, so as to meet the ecological management needs.

[0049] <Device>

[0050] Further, the present application also provides a device for determining the vegetation growth suitable area of the inland shore zone based on multi-source images, which can automatically implement the above-mentioned method, and the device is characterized in that it comprises:

[0051] A data set construction unit, which performs water-land area division and constructs a binary submergence image data set based on the images of the research reach;

[0052] An average submergence frequency calculation unit, which calculates the average submergence frequency of the shore beach duration based on the pixel scale by using steps 2.1-2.10 as follows:

[0053] Step 2.1: Select the submergence frequency calculation method of the shore zone according to whether the target reach has long-term water level observation data. If the water level observation data is available, go to step 2.2, otherwise go to step 2.10;

[0054] Step 2.2: Match the daily average water level of the date corresponding to each binary submergence image {B 11 ,.......,B ij} with the date as the index, and record it as {Z 11 ,.......,Z ij}, where Z ij represents the water level corresponding to the i-th image in the j-th year;

[0055] Step 2.3: Sort the daily average water level in T years from small to large according to the statistical demand of the regional submergence frequency, and divide it into m water level interval areas at equal intervals according to the arithmetic interval, 25≤m≤35;

[0056] Step 2.4: Statistically count the frequency F of the water level {Z 11 ,.......,Z ij} falling into the m water level interval areas, retain the water level interval areas with frequency F≥1, and merge the water level interval areas with frequency F=0 with the adjacent interval areas. The number of water level interval areas after merging is recorded as m', 20≤m’;

[0057] Step 2.5: According to the daily average water level data, the number of days N in the i-th year (i = 1, 2, …, T) in which the daily water level falls within the m' water level intervals is counted ik (k = 1, 2, …, m'), the number of days N ik is the length of the time period of the k-th water level interval in the i-th year;

[0058] Step 2.6: According to the water level data {Z 11 ,.......,Z ij} corresponding to the time sequence images, the number of images C in which the water level Z ij corresponding to the j-th image in the i-th year falls within the m' water level intervals is counted;

[0059] Step 2.7: According to step 2.5, the representative image of each water level interval in the i-th year is matched one by one, and the encoding is denoted as B ik , B ik represents the binary inundation image of the k-th water level interval in the i-th year;

[0060] Step 2.8: For all binary images B ik , the formula is used to perform pixel-by-pixel operations to obtain the inundation frequency of each pixel in the i-th year. The values of the inundation frequency of all pixels constitute the raster data FF i of the spatial distribution of the inundation frequency of the target river section in the i-th year;

[0061] Step 2.9: On the basis of the annual inundation frequency, pixel-by-pixel operations are performed according to the formula to obtain the spatial distribution raster data FF of the multi-year average inundation frequency in T years;

[0062] Step 2.10: If there is no water level observation data, the time coordinates of the binary image set {B 11 ,.......,B ij} are used as the basis for calculating the inundation duration. The time midpoint of the corresponding dates of the j-1-th and j-th images in the i-th year is denoted as t s , and the time midpoint of the corresponding dates of the j-th and j+1-th images is denoted as t e . The N ij = t e -t s is taken as the representative period of the j-th image. On this basis, the formula is used to perform operations to obtain the inundation frequency of each pixel, which together form the spatial distribution raster data FF i ' of the inundation frequency in the i-th year, j = 1, 2, …, n i , n iThe number of images collected in the ith year; the calculation method of the multi-year average inundation frequency spatial distribution grid data is the same as step 2.8;

[0063] The image set construction unit extracts the riparian vegetation coverage area from each image and constructs a binary coverage image set;

[0064] The frequency calculation unit calculates the riparian beach vegetation occurrence frequency in a long period based on the pixel scale;

[0065] The division determination unit divides the riparian vegetation in the target river section area into suitable zones and determines the critical inundation frequency threshold; according to the grid data of the multi-year average inundation frequency FF and the vegetation occurrence frequency VF of the target river section area, a scatter plot of FF~VF and a kernel density cloud plot are obtained, and then a suitable zone for the growth of riparian vegetation is divided, and the critical inundation frequency threshold for the suitable growth of riparian vegetation is derived;

[0066] The control unit is in communication with the data set construction unit, the average inundation frequency calculation unit, the image set construction unit, the frequency calculation unit, and the division determination unit, and controls the operation of them.

[0067] Preferably, the multi-source image-based determination device for the suitable growth zone of riparian vegetation in an inland river provided by the present application can further comprise an input display unit in communication with the control unit, allowing the operator to input operation instructions and displaying corresponding information according to the control instructions.

[0068] Preferably, in the multi-source image-based determination device for the suitable growth zone of riparian vegetation in an inland river provided by the present application, the input display unit can display the input, output data and processing process of the data set construction unit, the average inundation frequency calculation unit, the image set construction unit, the frequency calculation unit and the division determination unit in the form of a data table or an image in a static or dynamic manner according to the control instructions.

[0069] Effects of the invention

[0070] (1) Based on open multi-source remote sensing images, the data complement each other well, the time and spatial coverage are high, the image processing is based on a remote sensing cloud platform, the operation process is simple and quantitative, and the disadvantages of large field reconnaissance workload and strong experiential data processing are overcome.

[0071] (2) The calculation of inundation frequency requires less measured terrain data, the idea is clear, and the operation is simple. By making full use of the measured images and water level data year by year, considering the micro changes in riverbed topography caused by riverbed evolution between years in the scouring and silting balance state, avoiding the error caused by using a single terrain to calculate the inundation frequency of the riparian zone, greatly improving the calculation accuracy of the inundation frequency in the research period, and reflecting the real inundation frequency situation of the region more accurately, it is convenient for practical application.

[0072] (3) In the processing and statistical process of remote sensing image, the hydrological characteristics of natural river with large water level amplitude and strong randomness and the complex morphological characteristics of bank beach with ups and downs and water-land interlaced are fully considered, the calculation is more accurate and reliable, and has strong universality.

[0073] (4) Not only the spatial range of stable vegetation area, transition vegetation area and unsuitable growth area in the riparian zone can be obtained, but also the corresponding submergence frequency threshold of each area can be obtained, so that scientific basis can be provided for vegetation area configuration and cross section form reconstruction in ecological river management planning and design, and the needs of engineering practice are fully met. BRIEF DESCRIPTION OF DRAWINGS

[0074] Figure 1 The flow chart of the method for determining the vegetation growth suitable area of the inland riparian zone based on multi-source images according to the embodiment of the present application;

[0075] Figure 2 The comparison schematic diagram of actual riparian zone ground object distribution (a) and riparian zone ground object distribution (b) extracted based on water body index and vegetation index according to the embodiment of the present application;

[0076] Figure 3 The multi-year average submergence frequency distribution grid data graph (a) and the multi-year average vegetation appearance frequency distribution grid data graph (b) according to the embodiment of the present application;

[0077] Figure 4 The scatter point kernel density contour distribution schematic diagram of submergence frequency-vegetation appearance frequency (a) and the scatter point distribution and S-type fitting curve schematic diagram of submergence frequency-vegetation appearance frequency (b) according to the embodiment of the present application;

[0078] Figure 5 The riparian zone vegetation suitability zoning schematic diagram according to the embodiment of the present application. DETAILED DESCRIPTION

[0079] The specific implementation scheme of the method and device for determining the vegetation growth suitable area of the inland riparian zone based on multi-source images according to the present application will be described in detail below with reference to the drawings.

[0080] <Embodiment I>

[0081] This embodiment takes a certain river section as an example to illustrate the method of the present application. There is a large area of beach in the river section, and a hydrological station near the river section has long series of daily average water level observation data. The multi-spectral image data set is derived from the Landsat 5 / 7 Level-2 data product in the Google Earth Engine image database after preprocessing, and the time is from July 28, 1986 to January 1, 2004.

[0082] As Figure 1As shown, the method for determining suitable vegetation growth zones in inland river riparian zones based on multi-source imagery provided in this embodiment includes the following:

[0083] Step 1. Filter the image data and perform water and land area division on a frame-by-frame basis, constructing a binarized inundation dataset, including:

[0084] Step 1.1: Obtain multispectral images of the target river section from the Google Earth Engine data center over a period of 18 years, including 179 Landsat 5TM Level 2 SR images from 1986 to 2004 and 52 Landsat 7TM Level 2 SR images from 1999 to 2004. All image data have undergone preprocessing methods such as geometric correction, radiometric calibration, and atmospheric correction. The images were selected according to the following criteria: (1) the cloud cover in the target river section area of ​​the remote sensing image does not exceed 5%, and there are no defects or bands; (2) there are usable images for each year and each season. After screening the above-obtained images, a total of 163 Landsat images from 1986 to 2004 were found to be usable, of which 128 were Landsat 5 images and 35 were Landsat 7 images that could be identified.

[0085] Step 1.2: The Landsat Level 2 data is masked using the QA band, and after cropping and mosaicking, a regional multispectral image dataset {MS} is constructed. 11 MS 18j}, where i represents the i-th year and j represents the j-th image in the i-th year, for a total of 163 images.

[0086] Step 1.3 Select {MS 11 MS 18j Each image MS ij ρ GREEN Green band and ρ SWIR Shortwave infrared surface reflectance data are obtained using the formula MNDWI = (ρ GREEN -ρ SWIR ) / (ρ GREEN +ρ SWIR Band operations were performed to obtain the improved normalized difference water index (MNDWI) image, and {MS} was then used to perform band operations. 11 MS 18j The dataset was transformed into an improved normalized differential water index MNDWI image dataset {M 11 ,.......,M 18j}, where M ij This represents the j-th MNDWI image of year i.

[0087] Step 1.4 Based on Step 1.3, for the image dataset {M11 ,.......,M 18j Output the grayscale histogram of each image in the target river section area. The horizontal axis is the grayscale value (MNDWI value), and the vertical axis is the pixel frequency H under that grayscale value. The grayscale histogram should be bimodal.

[0088] Step 1.5: Based on the bimodal distribution characteristics of the gray-level histogram of each image, the Otsu's maximum inter-class variance method is used to determine the segmentation threshold of water bodies and non-water bodies for each image. The implementation process is as follows: (1) Perform classification trial calculations for each gray-level histogram. First, assume a gray-level threshold for classification, and then calculate according to the formula σ 2 =ω0×ω1×(μ0-μ1) 2 Calculate the inter-class variance σ under this grayscale threshold. 2 (1) The proportion of pixels in the image with a gray value less than the gray value is denoted as ω0, and its average gray value is μ0; the proportion of pixels with a gray value greater than the gray value is denoted as ω1, and its average gray value is μ1; (2) Repeat step (1) by continuously re-assuming the classification gray value threshold until the inter-class variance σ is found. 2 The grayscale threshold at which the maximum value is taken is the segmentation threshold between water and non-water bodies. In this embodiment, the segmentation threshold is determined using the above method, and a set of segmentation thresholds for water and non-water bodies is as follows: 0.25184, 0.31019, 0.22061, 0.27788, 0.18011, 0.26064, 0.28566, 0.20819...

[0089] Step 1.6: Apply the above threshold to the dataset {M} 11 ,.......,M 18j The images in the dataset are binarized frame by frame, with pixels in water areas assigned a value of 1 and pixels in non-water areas assigned a value of 0, resulting in a binarized flooded dataset {MB}. 11 ,.......,MB 18j}, where MB ij This represents the j-th multispectral binarized flooding image of the i-th year based on Landsat SR (30m resolution) data, totaling 163 binarized flooding images. For example... Figure 2 The image shown is a comparison between the binarized water body extraction results and the true-color image of this river section on April 14, 2003.

[0090] Step 1.7: In this embodiment, no corresponding Sentinel data is selected from the historical data, so there is no need to unify the pixel scale. The submerged dataset {MB} is binarized. 11 ,.......,MB 18j} Directly construct a 30m precision time-series binarized flooding dataset {B 11 ,.......,B 18j}

[0091] Step 2. The pixel-based coastal beach length duration average inundation frequency calculation includes:

[0092] Step 2.1: In this embodiment, 18 years of daily average water level observation data are obtained from a hydrological station near the river section, and the data is complete, which is transferred to step 2.2.

[0093] Step 2.2: Date is used as an index to match the daily average water level corresponding to the date of each binary inundation image {B 11 ,.......,B ij}, which is recorded as {Z 11 ,.......,Z 18j}, where Z ij represents the water level corresponding to the i-th image in the j-th year.

[0094] Step 2.3: In this embodiment, according to the regional inundation frequency statistical requirements, the daily average water level of 18 years is sorted from small to large, and 0.5m is used as the division interval to divide 34(25≤m≤35) water level interval.

[0095] Step 2.4: The frequency F of the water level {Z 11 ,.......,Z 18j} falling into the 34 water level intervals is counted, and the water level interval with frequency F≥1 is retained, and the water level interval with frequency F=0 is combined. In this embodiment, the higher water level interval(42.0-44.5m) and the lower water level interval(27.5-28.5m) are combined according to the F=0 condition of the image frequency, and the number of water level intervals after combination is 29, and 29(20≤m') water level intervals have high-quality image data.

[0096] Step 2.5: According to the annual daily average water level data, the number of days N ik (k=1,2,……,29) of the daily average water level falling into the 29 water level intervals in the i-th year(i=1,2,……,18) is counted, and the number of days N ik is the time period length of the k-th water level interval in the i-th year. In this example, the measured water level data of the hydrological station in 2001 is taken as an example, the minimum water level in 2001 is 29.2m, and the maximum water level is 38.9m, and the water level interval in the year only has distribution in the 29.0-39.0 water level interval, and the number of days N 18k of the water level appearing in each water level interval is counted: (29.0-29.5) days 19, (29.5-30.0) days 57, (30.0-30.5) days 30,……,(37.0-37.5) days 23, (37.5-38.0) days 13, (38.0-38.5) days 13, (38.5-39.0) days 3.

[0097] Step 2.6: Based on the time-series image water level data {Z 11 ,.......,Z 18j}, Statistical analysis of the daily average water level Z of the j-th image within the i-th year. ij The number of images C falling within the 29 water level intervals. If the number of images C ≥ 2 in the k-th interval of the year, then the image whose daily average water level is closest to the median water level of that interval is selected as the representative image for that interval. For example, if the water level interval is divided into (32m~33m), then the image with the water level closest to 32.5m is selected as the representative image. If the number of images C = 0 in the k-th interval of the year, then the daily average water level Z of the images throughout the year T is used. ij Take the average daily water level Z from the image. ij The image closest to the median water level of the given water level interval and closest in time to that year is used as the representative image for the k-th interval of that year. In this example, taking 2001 as an example, its representative image time and corresponding water level, and the number of days in the representative period are: 2001 / 1 / 10 (30.7) days 34, 2001 / 1 / 18 (30.3) days 30, 2001 / 5 / 10 (33.4) days 11, 2001 / 6 / 11 (35.4) days 10, 2001 / 7 / 5 (38.3) days Number 13, 2001 / 7 / 21 (34.2) days 23, replacement image 2000 / 3 / 12 (29.9) days 53, substitute 2002 / 6 / 22 (36.8) days 23, ..., 2001 / 8 / 22 (37.2) days 14, 2001 / 10 / 1 (37.6) days 13, 2001 / 9 / 7 (38.6) days 3.

[0098] Step 2.7: Following Step 2.6, match each water level interval of year i with its representative image, and encode it as B. ik B ik This represents the binarized inundation image of the k-th water level interval in the i-th year.

[0099] Step 2.8: For all binarized images B ik According to the formula Raster operations are performed pixel-by-pixel to obtain the inundation frequency of each pixel in year i. The inundation frequency values ​​of all pixels constitute the raster data FF of the spatial distribution of inundation frequency in the target river segment area in year i. i In this embodiment, the spatial distribution raster data of the annual inundation frequency of the target river section for a total of 18 years was obtained sequentially.

[0100] Step 2.9 Based on the annual flooding frequency over a total of 18 years, each pixel is calculated according to the formula. Calculations were performed to obtain the multi-year average inundation frequency spatial distribution raster data (FF) for this river section. Figure 3 .

[0101] Step 3. Extract the coastal vegetation cover area from each image and construct a binary vegetation cover image set.

[0102] Step 3.1: Based on the target river segment temporal multispectral image dataset {MS} obtained in Steps 1.1 and 1.2 11 MS 18j A total of 163 scenes were captured, and each image was taken as a MS. ij ρ RED Red light band and ρ NIR Infrared reflectance data is obtained using the formula NDVI = (ρ NIR -ρ RED ) / (ρ NIR +ρ RED Band operations were performed to obtain the Normalized Differential Vegetation Index (NDVI) imagery, and a long-term NDVI imagery dataset of the target river section was formed. 11 ,.......,V 18j}, where V ij This represents the j-th NDVI image in year i.

[0103] Step 3.2: For the NDVI image dataset {V 11 ,.......,V 18j Using a band calculator, pixel masks (NDVI≤0) are applied to the water area frame by frame. After masking, RNDVI images containing both light beaches and vegetation are generated, forming a dataset {RN}. 11 ,.......,RN 18j}

[0104] Step 3.3: {RN 11 ,.......,RN 18j 163 scenes in} ij The image grayscale histogram exhibits a bimodal shape. The threshold ' for segmenting the light beach and vegetation is calculated frame-by-frame using the OTSU maximum inter-class variance method described in step 1.5. In this embodiment, a set of RNDVI segmentation thresholds for light beach-vegetation are as follows: 0.43095, 0.40349, 0.37662, 0.32654, 0.23620, 0.39247, 0.49846, 0.38018, 0.42798, 0.53307, ...

[0105] Step 3.4 Segment the dataset {RN} according to the segmentation threshold ' for each image. 11 ,.......,RN 18j} are binarized and binary classification results are automatically generated, with the vegetation area pixel assigned as 1 and the beach pixel area assigned as 0. The binarized vegetation data set {MVB 11 ,.......,MVB 18j} is obtained, wherein MVB ij

[0106] represents the 30m-precision binarized vegetation image of the i-th year and the j-th image generated based on the multispectral data set {MS 11 ,.......,MS 18j} of Landsat data. In this embodiment, the accuracy evaluation of the extraction results of the vegetation water body and the beach is performed by using the confusion matrix, as shown in Table 1, the overall accuracy of the classification reaches 99.3%, and the kappa coefficient reaches 98.3%, which proves that the vegetation coverage extraction method provided by the present application has high accuracy and reliability in the shore zone feature classification. Figure 2 The binarized water body extraction result of the shore zone on April 14, 2003 is given.

[0107] Table 1 Accuracy evaluation of the confusion matrix

[0108]

[0109] Step 3.5 In this embodiment, the binarized vegetation data set {MVB 11 ,.......,MVB 18j} is directly constituted into the 30m-precision time-series binarized vegetation data set {P 11 ,.......,P 18j} without uniform pixel scale during the period without corresponding Sentinel data.

[0110] Step 4. Shore beach vegetation frequency calculation in a long period based on the pixel scale, including:

[0111] Step 4.1: The time coordinates of the binarized vegetation image set {P 11 ,.......,P 18j} are used as the basis for calculating the vegetation occurrence period, the time midpoint of the corresponding date of the i-th year and the j-th image is recorded as t s , the time midpoint of the corresponding date of the j-th and j+1-th image is recorded as t e , and N ij =t e -t s is taken as the representative period of the j-th image. On this basis, the vegetation occurrence frequency of each pixel is calculated according to the formula , wherein (j=1, 2, …, n i ), ni The number of images collected in the ith year, the vegetation frequency spatial distribution grid data VF of the ith year is obtained i In this embodiment, the vegetation frequency of 2001 is taken as an example, the vegetation coverage image time and its representative period N ij are: 2001 / 1 / 10, day 13, 2001 / 1 / 18, day 28, 2001 / 3 / 7, day 36, 2001 / 3 / 31, day 24, 2001 / 4 / 24, day 16, …, 2001 / 7 / 21, day 12, 2001 / 7 / 29, day 8, 2001 / 8 / 26, day 8, …, 2001 / 10 / 25, day 24, 2001 / 11 / 8, day 16, 2001 / 11 / 26, day 40. The annual vegetation frequency spatial distribution grid data VF of the 18-year riparian zone is obtained by sequentially calculating i .

[0112] Step 4.2: On the basis of the total of 18-year annual vegetation frequency VF i , the formula is operated pixel by pixel to obtain the spatial distribution of the average vegetation frequency of the river section for many years VF see Figure 3 .

[0113] Step 5. The division of the suitable area of the riparian zone in the target river section area includes:

[0114] Step 5.1: According to the grid data of the average submergence frequency FF and the vegetation frequency VF of the target river section area for many years, the submergence frequency and the vegetation frequency are outputted pixel by pixel, and stored in two arrays. The FF~VF relationship of all pixels in the two arrays is plotted, and a scatter plot is drawn, in which FF is the horizontal coordinate and VF is the vertical coordinate. The FF~VF scatter plot is shown in FIG. 5. Figure 4 .

[0115] Step 5.2: Cluster analysis is performed on the FF~VF scatter plot, and the generated scatter plot is converted into a probability density plot. The specific steps are as follows: the KS density function of MATLAB is used to perform kernel density estimation on the scatter plot data, a Gaussian kernel function is used to smooth the data, and the estimated density is calculated. After normalizing the estimated probability density, visualization is performed, and an FF~VF kernel density estimation cloud chart is drawn. The value in the cloud chart represents the proportion of pixels under the condition of a certain submergence frequency and a certain vegetation frequency combination to the total pixels. The FF~VF kernel density estimation cloud chart is shown in FIG. 6. Figure 4 .

[0116] Step 5.3: Add isograms to the FF~VF kernel density cloud chart, generally three partitions can be formed according to the closed condition of the 2%~5% isograms, and the projection of the boundary between any two of the partitions on the longitudinal axis corresponds to two vegetation frequency thresholds CVF1, CVF2. In this embodiment, according to the closed condition of the 2% isogram, CVF1=20%, CVF2=60%, and accordingly the VF grid data can be divided into three categories, corresponding to three vegetation frequency levels, i.e. high, medium and low, wherein the pixels with VF≥60% correspond to the stable vegetation area; the area with 20%<VF<60% corresponds to the transition vegetation area, and the area with VF≤20% corresponds to the vegetation unsuitable growth area. The visualization of these vegetation partitions in combination with geographic spatial information can produce a vegetation suitability partition map of the riparian zone, which is shown in Fig. 2. Figure 5 .

[0117] Step 5.4: For the FF~VF scatter chart, FF is divided into 100 equidistant intervals with a step value of 1%, the mean value of FF in each interval is FFA, and the mean value of VF in each interval is VFA. The logistic S-shaped curve is fitted according to FFA and VFA by using the origin software, and in this embodiment, the fitting result of the logistic S-shaped curve is y=0.33598+90.72094 / (1+(x / 14.74212) 2.20445 ), wherein the center position of the flood frequency value is 14.74212. The scatter chart of FF~VF and the fitted logistic S-shaped curve are shown in Fig. 3. Figure 4 .

[0118] Step 5.5: The two vegetation frequency thresholds CVF1=20% and CVF2=60% defined in step 5.2 are mapped to the logistic S-shaped curve, and the critical flood frequency thresholds CFF1=11% and CFF2=26% between the suitable area, the transition area and the unsuitable area of the riparian zone vegetation in the target river reach are obtained. In engineering practice, the riparian terrain can be appropriately modified according to these critical values, and the flood frequency of a certain position can be made greater or smaller than the above critical values by changing the terrain elevation, so as to meet the ecological management requirements.

[0119] <Embodiment Two>

[0120] The embodiment two provides a device for determining the vegetation growth suitable area of the riparian zone in an inland river based on the method of the present application, which comprises a data set construction unit, an average flood frequency calculation unit, an image set construction unit, an occurrence frequency calculation unit, a division determination unit, an input display unit and a control unit.

[0121] The data set construction unit can perform the content described in step 1 above, and based on the images of the research river reach, the water-land area division is implemented and the binary flood image data set is constructed.

[0122] The average inundation frequency calculation unit can perform the content described in step 2 above, and calculate the average inundation frequency of the beach length of the riparian zone based on the pixel scale.

[0123] The image set construction unit can perform the content described in step 3 above, and extract the vegetation coverage area of the riparian zone from each image and construct a binary coverage image set.

[0124] The occurrence frequency calculation unit can perform the content described in step 4 above, and calculate the vegetation occurrence frequency of the riparian zone beach in a long period based on the pixel scale.

[0125] The division determination unit can perform the content described in step 5 above, and divide the riparian vegetation in the target river section area and determine the critical inundation frequency threshold; according to the grid data of the average inundation frequency FF and the vegetation occurrence frequency VF of the target river section area in many years, a scatter plot of FF~VF and a kernel density cloud map are obtained, and then the suitable area of the riparian vegetation growth is divided, and the critical inundation frequency threshold of the suitable growth of the riparian vegetation is derived.

[0126] The input display unit is connected in communication with the control unit, so that the operator can input operation instructions, and display corresponding information according to the control instructions. Specifically, the input display unit can display the input and output data and processing process of the data set construction unit, the average inundation frequency calculation unit, the image set construction unit, the occurrence frequency calculation unit, and the division determination unit in the form of a data table or an image statically or dynamically according to the control instructions.

[0127] The control unit is connected in communication with the data set construction unit, the average inundation frequency calculation unit, the image set construction unit, the occurrence frequency calculation unit, the division determination unit, and the input display unit, and controls the operation of them.

[0128] The above embodiments are only examples of the technical solutions of the present application. The method and device for determining the suitable growth area of the riparian vegetation in the inland river based on multi-source images are not limited to the content described in the above embodiments, but are limited to the scope defined in the claims. Any modification, supplement or equivalent replacement made by the person skilled in the art on the basis of the embodiments is within the scope claimed by the claims of the present application.

Claims

1. A method for determining a suitable area for vegetation growth in an inland river riparian zone based on multi-source images, characterized in that, The method comprises the following steps: Step 1: based on the image data of the study reach, water-land area division is implemented frame by frame, and a binary inundation image data set is constructed; Step 2: based on pixel scale, the long-term average inundation frequency of the beach of the riparian zone is calculated; specifically comprising the following sub-steps: Step 2.1: according to whether the target river reach has long series water level observation data, the calculation method of the inundation frequency of the riparian zone is selected, if the water level observation data is available, step 2.2 is entered, otherwise step 2.10 is entered; Step 2.2: Match each binary inundation image with the date the average daily water level of the corresponding date, denoted as where is the water level of the corresponding date of the th image in the th year. Step 2.3: According to the statistical requirements of regional flood frequency, the multi-year daily average water levels in T years are sorted from small to large, and then divided into equal intervals by arithmetic intervals water level interval, ; Step 2.4: Statistics of water level Fall in The frequency of water level interval , the frequency of water level interval, and the adjacent interval are combined, and the number of water level intervals after combination is recorded as , , ;​ Step 2.5: According to the daily average water level data of each year, the number of days in each year in which the daily water level falls in each water level interval is counted year , the number of days in the period of the year in which the daily water level falls in the water level interval​​​​​ Step 2.6: Corresponding water level data of time-series images , the number of images corresponding to water level data falling in each water level interval in the year ;​​​​ Step 2.7: Follow step 2.5 to... Each water level interval for a given year was matched with its representative image, and coded as follows: , Indicates the first Year Binarized flooding images of each water level range; Step 2.8: for all binary images , the formula is calculated pixel by pixel to obtain the submergence frequency of each pixel in the year , and the values of the submergence frequency of all pixels constitute the raster data of the spatial distribution of the submergence frequency of the target river section in the year ;​ Step 2.9: On the basis of the annual inundation frequency, the pixel-by-pixel operation is performed according to the formula to obtain the multi-year average inundation frequency spatial distribution grid data in T years ; Step 2.10: If water level observation data is unavailable, use a binarized image set. The time coordinates are used as the basis for calculating the flooding duration, and the first... Year Page and the first The midpoint of the time corresponding to the date of the image is denoted as , No. Page and the first The midpoint of the time corresponding to the date of the image is denoted as ,Pick As the first The representative time period of each image, based on this formula The submersion frequencies of each pixel are obtained through calculation, and they together form the first... Raster data of spatial distribution of annual inundation frequency , j =1,2,......, n i , n i For the first i The number of images collected annually; the calculation method for the multi-year average flooding frequency spatial distribution raster data is the same as step 2.8; Step 3: the vegetation coverage area of the riparian zone is extracted frame by frame from the image, and a binary coverage image set is constructed; Step 4: based on pixel scale, the long-term vegetation appearance frequency of the beach of the riparian zone is calculated; Step 5: Divide the riparian vegetation within the target river section into suitable zones; based on the multi-year average inundation frequency of the target river section. and frequency of vegetation raster data obtained The scatter plot and kernel density cloud map were used to delineate suitable areas for riparian vegetation growth and to deduce the critical submergence frequency threshold for suitable riparian vegetation growth.

2. The method for determining the vegetation growth suitable area of the riparian zone in the inland river according to claim 1, characterized in that: wherein, In step 1, the 10m-accuracy binary inundation image dataset and is converted to pixel scale of 30m using resampling technique, followed by merging and into 30m-accuracy binary inundation image dataset where denotes the th image of the th year.

3. The method for determining the vegetation growth suitable area of the riparian zone in the inland river according to claim 1, characterized in that: wherein, Step 1 comprises the following sub-steps: Step 1.1: from the data center, image data of the target river reach for T years is obtained, T≥10, wherein multi-spectral image data subjected to geometric correction, radiation calibration and atmospheric correction preprocessing is mainly used, and synthetic aperture radar data subjected to radiation positioning, correction and denoising is used as a supplement, and multi-spectral image data meeting the quality and quantity conditions is selected: (1) in the remote sensing image, the target river reach area should meet the following conditions: no more than 5% of the cloud cover is blocked, no damage or strip; (2) within each year, there should be available images in each season; Step 1.2: Mask the cloud layer using the QA band for Landsat SR secondary data, and the SLC band for Sentinel-2 Level 2A data; smooth and denoise the SAR GRD data; after removing the cloud, crop the target river section to form the target river section time-series multi-spectral image dataset , time-series synthetic aperture radar dataset , wherein represents the year , represents the image of the th th image of the year ; , is the number of images collected in the year i ; Step 1.3: Select Each image of Green band and The reflectivity data for the shortwave infrared band is calculated according to the formula. By performing band calculations, the improved normalized differential water index is obtained. Images, thus making Dataset transformed into improved normalized differential water index Image dataset ,in Indicates the first The year's first frame Imagery; for synthetic aperture radar datasets Select each image Dual-polarization data from Level 1 ground distance detection products that have undergone projection correction are analyzed according to the water extraction index formula. The water extraction index is obtained by calculating the band data. ,form Image dataset ,in Indicates the first The year's first frame image; Step 1.4: Based on the image data set of step 1.3, output the gray level histogram of each image of the target river section area, with the horizontal coordinate being the gray level value and the vertical coordinate being the pixel frequency of the gray level value and The gray level histogram should be bimodal ​ Step 1.5: based on the bimodal distribution characteristics of the gray level histogram of each image, the OTSU maximum between-class variance method is used to determine the segmentation threshold of water and non-water body frame by frame; Step 1.6: Calculate water and non-water segmentation threshold for each image After that, segment the dataset according to the segmentation threshold of each image and Binaryzation for each image, water region pixel is assigned as 1, non-water region pixel is assigned as 0, get binaryzation inundation dataset and where represents the binaryzation inundation image of the th image in the th year based on Landsat SR and Sentinel-22A data; represents the binaryzation inundation image of the th image in the th year based on Sentinel-1 SAR data. Step 1.7: Binary inundation image dataset and 10m accuracy binary images generated from Sentinel-2 2A / Sentinel-1 SAR satellites, the pixel scale is converted to 30m using resampling techniques, followed by and merged into a 30m accuracy binary inundation image dataset where represents the th image of the th year.

4. The method for determining the vegetation growth suitable area of the riparian zone in the inland river according to claim 1, characterized in that: wherein In step 2.6, if the first Number of images within the interval Then, the image whose daily average water level is closest to the median of that interval is taken as the representative image of that interval; if the first day of a certain year... Number of images within the interval Then iterate through the images corresponding to the water levels within year T. Take the water level corresponding to the image The image closest to the median of that water level range and closest in time to that year will be used as the first image of that year. The interval represents the replacement of the image.

5. The method for determining the vegetation growth suitable area of the riparian zone in the inland river according to claim 1, characterized in that: wherein Step 3 comprises the following sub-steps: Step 3.1: Obtain the time-series multi-spectral image dataset of the target river reach based on the image datasets obtained in steps 1.1 and 1.2 , take the reflectance data of each image of the red light band and infrared band , and perform band operation according to the formula to obtain the normalized difference vegetation index image, and form a long time-series vegetation index image dataset of the target river reach region , wherein represents the image of the th year Step 3.2: For each image in the dataset image dataset , mask water region pixels using band calculator , generate image containing beach + vegetation after masking image, form dataset ; Step 3.3: In the image gray histogram presents bimodal, using step 1.5 in the maximum inter-class variance method calculated out of the threshold of the light beach and vegetation image gray histogram presents bimodal, using step 1.5 in the maximum inter-class variance method calculated out of the threshold of the light beach and vegetation ; Step 3.4: Segmentation threshold for each image Data set Each image is binarized, with the vegetation region pixels assigned a value of 1 and the mudflat region pixels assigned a value of 0, to obtain a binary vegetation data set wherein represents the year the binary vegetation image; Step 3.5: Pixel scale conversion is performed on the 10m precision binary canopy cover imagery using resampling techniques to convert its pixels from 10m to 30m, forming a 30m precision binary canopy cover dataset .

6. The multi-source imagery based inland river riparian zone vegetation growth suitability zone determination method of claim 1, characterized in that: wherein, Step 4 comprises the following sub-steps: Step 4.1: Use the binarized overlay image set The time coordinates are used as the basis for calculating the duration of vegetation appearance, and the first... Year Page and the first The midpoint of the time corresponding to the date of the image is denoted as , No. Page and the first The midpoint of the time corresponding to the date of the image is denoted as ,Pick As the first The representative time period of each image; Based on this, according to the formula The frequency of vegetation occurrence in each pixel is calculated to obtain the first... Annual vegetation occurrence frequency spatial distribution raster data ; Step 4.2: Based on the frequency of vegetation occurrence year by year, the formula is used to calculate the spatial distribution of the average frequency of vegetation occurrence in T years. Step 4.2: Based on the frequency of vegetation occurrence year by year, the formula is used to calculate the spatial distribution of the average frequency of vegetation occurrence in T years.

7. The method for determining the vegetation growth suitable area of the riparian zone in the inland river according to claim 1, characterized in that: wherein Step 5 comprises the following sub-steps: Step 5.1: Based on the multi-year average inundation frequency of the target river section. and frequency of vegetation The raster data is processed pixel-by-pixel, outputting the flooding frequency and vegetation occurrence frequency, and stored in two arrays; based on the data of all pixels in the two arrays... Plot the relationships using a scatter plot; Step 5.2: For Cluster analysis is performed on the scatter plot, and the resulting scatter plot is converted into a kernel density cloud plot; Step 5.3: For Contour lines were added to the kernel density cloud map, and three zones were formed based on the closure of the 2%~5% contour lines. The projections of the boundaries between each pair of the three zones onto the vertical axis corresponded to two vegetation occurrence frequency thresholds. , Based on this, The raster data is divided into three categories, corresponding to high, medium, and low vegetation occurrence frequencies. The pixel corresponds to a stable vegetation zone. The area is a transitional vegetation zone. These areas are designated as unsuitable for vegetation growth. By combining geospatial information, these vegetation zones are visualized to create a riparian vegetation suitability zoning map. Step 5.4: For each of the 100 bins, calculate the mean value of the scatter plot, the bins are divided into 100 equidistant intervals with a step value of 1%, the mean value of each interval is , the mean value of each interval is , and a logistic S-shaped curve is fitted according to and ;​​ Step 5.5: Two vegetation occurrence frequencies are calculated for each cell , Corresponding to the logistic S-curve, the critical submergence frequencies between the suitable, transition and unsuitable zones for vegetation are derived , .

8. The device for determining the suitable area of vegetation growth in the inland river shore zone based on multi-source images, characterized in that, including: The data set construction unit, based on the image of the study reach, implements water-land area division frame by frame, and constructs a binary inundation image data set; The average inundation frequency calculation unit, based on pixel scale, calculates the long-term average inundation frequency of the beach of the riparian zone by using the following steps 2.1-2.10; Step 2.1: according to whether the target river reach has long series water level observation data, the calculation method of the inundation frequency of the riparian zone is selected, if the water level observation data is available, step 2.2 is entered, otherwise step 2.10 is entered; Step 2.2: Match each binary inundation image with the date the average daily water level of the corresponding date, denoted as where is the water level of the corresponding date of the th image in the th year. Step 2.3: According to the statistical requirements of regional flood frequency, the multi-year daily average water levels in T years are sorted from small to large, and then divided into equal intervals by arithmetic intervals water level interval, ; Step 2.4: Statistics of water level Fall in The frequency of water level interval , the water level interval of the frequency , and the adjacent interval are combined, and the number of water level intervals after combination is recorded as , ; ; Step 2.5: According to the daily average water level data of each year, the number of days in each year in which the daily water level falls in each water level interval is counted. year ​​​​​​ Step 2.6: Corresponding water level data of time-series images , the number of images corresponding to water level data falling in each water level interval in the year ;​​​​ Step 2.7: Match the representative image of each water level interval in the year 2008, and encode it as , , represents the binary image of the flood in the year 2008 in the water level interval of the 1th water level interval. ​​ Step 2.8: For all binarized images According to the formula Perform operations pixel by pixel to obtain the first... The flooding frequency of each pixel in a given year, and the flooding frequency values ​​of all pixels constitute the first... Raster data on the spatial distribution of inundation frequency in the target river section area for the year. ; Step 2.9: On the basis of the annual inundation frequency, the pixel-by-pixel operation is performed according to the formula to obtain the multi-year average inundation frequency spatial distribution grid data in T years ; Step 2.10: If there is no water level observation data, use the binary image set The time coordinate of the image as the basis for calculating the duration of submergence, the time midpoint of the corresponding date of the first and the second image is recorded as , and the time midpoint of the corresponding date of the first and the second image is recorded as , and is taken as the representative period of the first image. On this basis, the submergence frequency of each pixel is calculated according to the formula , and they collectively form the submergence frequency spatial distribution grid data of the first year, , j =1,2,......, n i , n i is the number of images collected in the first i year; the calculation method of multi-year average submergence frequency spatial distribution grid data is the same as step 2.8; The image set construction unit, the vegetation coverage area of the riparian zone is extracted frame by frame from the image, and a binary coverage image set is constructed; The appearance frequency calculation unit, based on pixel scale, calculates the long-term vegetation appearance frequency of the beach of the riparian zone; The delineation and determination process involves dividing the riparian vegetation within the target river section into suitable zones and determining the critical inundation frequency threshold; based on the multi-year average inundation frequency of the target river section... and frequency of vegetation raster data obtained The scatter plot and kernel density cloud map were used to delineate suitable areas for riparian vegetation growth and to deduce the critical inundation frequency threshold for suitable riparian vegetation growth. The control unit, in communication with the data set construction unit, the average inundation frequency calculation unit, the image set construction unit, the appearance frequency calculation unit and the division determination unit, controls the operation of them.

9. The multi-source imagery based inland river riparian zone vegetation growth suitability zone determination apparatus of claim 8, wherein, Also included are: The input display unit is connected to the control unit for inputting operation instructions from the operator and displaying corresponding information according to the control instructions. 10.The inland river riparian zone vegetation growth suitability area determination device based on multi-source images according to claim 9, characterized in that: wherein, The input display unit can statically or dynamically display the input, output data and processing process of the data set construction unit, the average submergence frequency calculation unit, the image set construction unit, the occurrence frequency calculation unit and the division determination unit in the form of a data table or an image according to the control instructions.

Images