Automatic extraction method of aquatic vegetation and algal blooms in eutrophic lakes based on landsat images

Abstract

The application provides an automatic extraction method of eutrophic lake aquatic vegetation and algal bloom based on Landsat images, carries out a teanga transformation, NDVI index calculation and FAI index calculation on the Landsat image data respectively, extracts an AVI image based on the third component after the teanga transformation, and determines a classification threshold a of the AVI image for classifying vegetation and non-vegetation by using a spectral linear mixing model, determines a classification threshold b of the NDVI image and a classification threshold c of the FAI image respectively; a decision tree classification model is established according to the AVI, NDVI and FAI index values and the classification thresholds a, b and c, and an automatic extraction result map of different aquatic vegetation communities and algal bloom of the eutrophic lake is generated. The application can automatically, real-timely, large-scaledly and high-precisely identify and monitor the spatial distribution of different vegetation communities (floating leaf / emergent and submerged vegetation) and blue-green algal bloom of the lake, and has important practical significance for monitoring and tracing the mutual conversion dynamic process of the grass-type lake and the algal-type lake, lake ecological restoration management and grass and algae salvage.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing applications and relates to an automatic extraction method for aquatic vegetation and algal blooms in eutrophic lakes based on Landsat imagery. Background Technology

[0002] Aquatic vegetation and algal blooms are both primary producers in aquatic ecosystems, but their ecological roles in lake ecosystems are completely different. On the one hand, aquatic vegetation plays a positive role in maintaining the balance of lake ecosystems, material cycling, and water quality. On the other hand, cyanobacterial blooms are a result of severe eutrophication in lakes, posing a great threat to the ecosystem, damaging the lake landscape, seriously affecting the development of fisheries and tourism, and, more fatally, the toxic substances produced by algal blooms may enter the human body through the food chain, causing great harm to human health. Therefore, as important regulators and indicators of lake ecosystems, real-time and rapid monitoring of aquatic vegetation communities (floating-leaved / emergent and submerged vegetation) and algal blooms is of guiding significance for assessing the ecological functions of lake ecosystems. At the same time, the dynamic process of the transformation between grass-type and algal-type lakes, lake ecological restoration management, and the removal of grasses and algae are of important practical significance.

[0003] Traditional methods for monitoring aquatic vegetation and algal blooms include manual quadrat and transect surveys, which are time-consuming, labor-intensive, have limited monitoring range and sample size, and are difficult to monitor (many aquatic vegetation and algal bloom areas are difficult for people to reach), making it impossible to quickly, in real time, and over a large area obtain the types and distribution of aquatic vegetation and algal blooms. Remote sensing has the advantage of rapid, real-time, and large-area monitoring of surface information, but in eutrophic lakes, the difficulties and challenges in monitoring aquatic vegetation communities and algal blooms are: (1) Aquatic vegetation, especially floating-leaved vegetation, has a similar spectral curve to algal blooms, making it difficult to distinguish; (2) Submerged vegetation is located below the water surface, has a weak spectral signal, and is difficult to distinguish from open water bodies. Summary of the Invention

[0004] The purpose of this invention is to overcome the difficulties in remote sensing extraction of aquatic vegetation and algal blooms in eutrophic lakes in the prior art, and to provide an automatic extraction method for aquatic vegetation and algal blooms in eutrophic lakes based on Landsat imagery.

[0005] The above-mentioned technical objective of the present invention is achieved through the following technical solution:

[0006] An automatic extraction method for aquatic vegetation and algal blooms in eutrophic lakes based on Landsat imagery includes:

[0007] Landsat satellite surface reflectance imagery data covering the study area were acquired and processed as follows;

[0008] i. Perform tassel transformation on the image data, use the third component after tassel transformation as the aquatic vegetation index AVI, acquire AVI images, and use the spectral linear mixture model to determine the classification threshold a in the AVI images for classifying aquatic and non-aquatic vegetation.

[0009] ii. Obtain Normalized Difference Vegetation Index (NDVI) images based on image data to distinguish between floating / emergent vegetation and submerged vegetation in aquatic vegetation areas, and determine the classification threshold b in the NDVI images;

[0010] iii. Based on image data, obtain phytoplankton index (FAI) images to distinguish algal blooms from water bodies in non-aquatic vegetation areas and determine the classification threshold c in the FAI images;

[0011] Based on the AVI, NDVI, and FAI index values ​​and their classification thresholds a, b, and c, the discrimination conditions are determined, a decision tree classification model is established, and automatic extraction results of different aquatic vegetation communities and algal blooms in eutrophic lakes are generated.

[0012] As a preferred implementation, Landsat satellite surface reflectance image data of the covered study area are obtained from the GEE cloud platform.

[0013] As a preferred embodiment, the method further includes preprocessing the Landsat satellite surface reflectance image data; the preprocessing includes image cloud removal and study area cropping.

[0014] In a preferred embodiment, the image declouding is performed by using the QA band and setting the values ​​of cloudShadowBitMask and cloudsBitMask to 0 before masking.

[0015] As a preferred implementation, the classification threshold a is determined based on the following formula:

[0016] a = p × AVI(w) + (1-p) × AVI(v)

[0017] Wherein, AVI(w) and AVI(v) represent the average AVI values ​​of pixels in pure water and pixels in densely growing floating leaf / emergent vegetation areas, respectively, and p is the proportion of pure water endmembers.

[0018] As a preferred implementation, AVI(w) is obtained based on the following steps:

[0019] 1) Obtain NDWI images of the study area and calculate their histograms;

[0020] 2) The NDWI value corresponding to the second highest peak in the NDWI histogram is used as the NDWI threshold for extracting pure water pixels;

[0021] 3) Extract all pixels with NDWI values ​​greater than the threshold as pure water pixels, and calculate the average AVI value of these pixels as the value of AVI(w).

[0022] As a preferred implementation, AVI(v) is obtained based on the following steps:

[0023] 1) By combining visual interpretation with measured data, select multiple pixels from images covering different lakes where floating leaves / emergent growth is dense;

[0024] 2) Calculate the average AVI value of these pixels and use this value as the value of AVI(v).

[0025] In a preferred embodiment, the classification threshold b and classification threshold c are determined based on the maximum gradient method.

[0026] As a preferred implementation method, the steps for determining the classification threshold b based on the maximum gradient method are as follows:

[0027] 1) Generate the gradient image of the NDVI image and define the gradient of a pixel as the difference between it and its neighboring pixels in a 3×3 window;

[0028] 2) Remove densely packed submerged and floating / emergent areas, retaining pixels in the boundary area between the two;

[0029] 3) Find the pixel with the largest gradient value among the remaining pixels, and use the NDVI value corresponding to that pixel as the NDVI threshold;

[0030] 4) Perform histogram statistics on the NDVI thresholds of multiple images, and use the average histogram minus twice the standard deviation as the classification threshold b of the NDVI map.

[0031] As a preferred implementation method, the steps for determining the classification threshold c based on the maximum gradient method are as follows:

[0032] 1) Generate the gradient image of the FAI image and define the gradient of a pixel as the difference between it and its neighboring pixels within a 3×3 window;

[0033] 2) Remove clean water bodies and areas with high concentrations of algal blooms, while retaining pixels in the boundary area between the two;

[0034] 3) Find the pixel with the largest gradient value among the remaining pixels, and use the FAI value corresponding to that pixel as the FAI threshold;

[0035] 4) Perform histogram statistics on the FAI threshold of multiple images, and use the average histogram minus twice the standard deviation as the classification threshold c of the FAI image.

[0036] As a preferred implementation, the decision tree classification model is as follows:

[0037] The study area is divided into aquatic vegetation and non-aquatic vegetation based on the AVI index. When a pixel in the study area satisfies AVI≥a, it is classified as aquatic vegetation; otherwise, it is classified as non-aquatic vegetation.

[0038] For aquatic vegetation, the NDVI index is used for classification. When NDVI≥b, it is classified as a floating-leaved / emergent plant; otherwise, it is classified as a submerged plant.

[0039] For non-aquatic vegetation, the FAI index is used for classification. When FAI ≥ c, it is classified as an algal bloom; otherwise, it is classified as a water body.

[0040] As a preferred implementation method, the accuracy of aquatic vegetation and algal bloom extraction is evaluated using overall accuracy and Kappa coefficient.

[0041] This invention first divides the study area into aquatic and non-aquatic vegetation based on the third component of the tasseling transformation. It then designs a spectral linear mixture model to calculate the classification threshold. Subsequently, it further classifies aquatic and non-aquatic vegetation using the NDVI and FAI indices. This overcomes the problems of floating-leaved vegetation and algal blooms having similar spectral curves, making them difficult to distinguish, and submerged vegetation being below the water surface with weak spectral signals, making it difficult to distinguish from open water bodies. The method of this invention can automatically, in real-time, over a wide area, and with high precision identify and monitor the spatial distribution of different vegetation communities (floating-leaved / emergent and submerged vegetation) and cyanobacterial blooms in lakes. It has significant practical implications for monitoring and tracing the dynamic process of the transformation between grass-like and algal-like lakes, lake ecological restoration management, and the removal of grasses and algae. Attached Figure Description

[0042] Figure 1 This is a flowchart of the technical solution of the present invention.

[0043] Figure 2 This is a histogram of the NDVI and FAI thresholds in an embodiment of the method described in this invention.

[0044] Figure 3 The remote sensing extraction results of aquatic vegetation communities and algal blooms in Taihu Lake on September 27, 2021, are shown in the embodiments of the method described in this invention. Detailed Implementation

[0045] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0046] This paper presents an automatic extraction method for aquatic vegetation and algal blooms in eutrophic lakes based on Landsat imagery, using Taihu Lake, a typical large shallow lake in China (average depth 1.9 meters), as an example. Taihu Lake's latitude and longitude ranges from 119.55–120.34E and 30.55–31.32N. Surveys and data show that large areas of aquatic vegetation, mainly emergent / floating-leaved plants and submerged plants, are distributed in the eastern part of Taihu Lake, while large-scale algal blooms frequently occur in the central and northwestern parts of the lake. Currently, due to intense human activities in recent years, the aquatic vegetation and wetland environment of the lake have been severely disturbed and damaged. Real-time and rapid monitoring of the distribution area and size of aquatic vegetation and algal blooms in Taihu Lake using remote sensing technology is of significant guiding importance for lake ecological restoration and algal harvesting. The following is a case study using Taihu Lake as the research area, with the specific process as follows: Figure 1 As shown, it includes:

[0047] Step 1: Based on the acquisition time of the measured sample points, download the Landsat satellite surface reflectance image data (LANDSAT / LC08 / C01 / T1_SR) covering Taihu Lake, synchronized with the measured sample points, from Google Earth Engine (GEE). Perform image declouding on the GEE cloud platform, then crop the corresponding study area image from the image using the Taihu Lake vector boundary. Image declouding is performed using the QA band, setting the values ​​of cloudShadowBitMask and cloudsBitMask to 0 before masking.

[0048] Step 2: Perform tassel transformation on the study area image on September 27, 2021, and use the third component after tassel transformation as the aquatic vegetation index AVI to obtain the corresponding AVI image to distinguish between aquatic and non-aquatic vegetation.

[0049] Specifically, the tassel transformation is performed first using the following formula:

[0050] y = Cx + α

[0051] α represents a constant offset to avoid negative values ​​during the transformation process; x is the pixel vector of the Landsat multispectral image before transformation; y is the pixel vector of the Landsat multispectral image after transformation; C is the transformation matrix, as follows:

[0052]

[0053] Therefore, the formula for calculating the third component, Wetness (humidity index), is as follows:

[0054] Wetness = 0.1511R Blue +0.1973R Green+0.3283R Red +0.3407R NIR -0.7117R Swir1

[0055] -0.4559R Swir2

[0056] The formula for calculating AVI is:

[0057] AVI = Wetness

[0058] Considering that AVI images are obtained by weighting the spectra of each band, multiplicative noise (such as atmospheric effects) cannot be eliminated, and they are highly sensitive to changes in the atmospheric or aquatic environment, making it impossible to obtain a universal (fixed) threshold using the traditional maximum gradient method. Therefore, this invention introduces a spectral linear mixture model to calculate the classification threshold 'a' of AVI images. This method can automatically obtain AVI thresholds for classifying aquatic and non-aquatic vegetation based on different AVI image features; that is, it adaptively obtains different classification thresholds 'a' based on different AVI images.

[0059] The calculation formula for the spectral linear mixing model is as follows:

[0060] S m,i (λ)=p×S W,i (λ)+(1-p)×S V (λ)

[0061] Among them, S m,i (λ) represents the mixed spectrum, S W,i (λ) and S V (λ) represents the spectrum of the endmembers of pure water and pure aquatic vegetation, and p is the proportion of the endmembers of pure water, which is 80%.

[0062] Based on the mixed spectrum, the threshold α of AVI is calculated using the following formula:

[0063] a = p × AVI(w) + (1-p) × AVI(v)

[0064] Where 'a' is the AVI threshold, and AVI(w) and AVI(v) represent the average AVI values ​​of pixels in pure water and pixels in areas with dense floating / emergent vegetation, respectively. The NDWI image of the scene is acquired, and the NDWI value corresponding to the second highest peak in the NDWI histogram is used as the threshold for extracting pure water; this value is calculated to be 0.754. All pixels with NDWI values ​​greater than 0.754 are extracted as pure water pixels.

[0065] Specifically, the NDWI calculation formula is as follows:

[0066] NDWI = (R GREEN -R NIR) / (R GREEN +R NIR )

[0067] Among them, R GREEN It is in the green band, with a wavelength range of 0.525-0.600 micrometers; R NIR It is in the near-infrared band, with a wavelength range of 0.845-0.885 micrometers.

[0068] Based on the measured data, 200 pixels with dense floating leaf / emergent growth were selected from images covering different lakes through visual interpretation.

[0069] By statistically analyzing the average AVI values ​​of pixels in pure water bodies and areas with dense growth of floating / emergent vegetation, the AVI(w) value was found to be -0.0488 and the AVI(v) value was found to be 0.02452. The two were then mixed according to the spectral linear mixing model at a ratio of 80%:20%, and the final AVI threshold a was found to be -0.03413.

[0070] Step 3: Obtain NDVI images of 60 scenes on Taihu Lake according to the NDVI calculation formula to distinguish between submerged vegetation and floating / emergent vegetation in aquatic vegetation areas. Use the maximum gradient method to determine the NDVI threshold of each scene.

[0071] Specifically, the NDVI calculation formula is as follows:

[0072] NDVI = (R NIR -R RED ) / (R NIR +R RED )

[0073] Among them, R RED It is in the red band, with a wavelength range of 0.63-0.68 micrometers; R NIR It is in the near-infrared band, with a wavelength range of 0.845-0.885 micrometers.

[0074] The classification threshold b is determined for the NDVI image using the maximum gradient method, and the steps are as follows:

[0075] 1) Generate the gradient image of the NDVI image and define the gradient of a pixel as the difference between it and its neighboring pixels in a 3×3 window;

[0076] 2) Remove densely packed submerged and floating / emergent areas, retaining pixels in the boundary area between the two;

[0077] 3) Find the pixel with the largest gradient value among the remaining pixels, and use the NDVI value corresponding to that pixel as the NDVI threshold;

[0078] 4) Perform histogram statistics on the NDVI thresholds of multiple images, and obtain the NDVI classification threshold b by subtracting twice the standard deviation from the mean of the histogram.

[0079] Histogram analysis of the NDVI thresholds of 60 images showed that the results followed a normal distribution, as follows: Figure 2 As shown. Subtracting twice the standard deviation from the mean of the NDVI threshold histogram yields an NDVI classification threshold b of 0.2.

[0080] Step 4: Obtain FAI images of 50 scenes on Taihu Lake according to the FAI calculation formula to distinguish water bodies and algal blooms in non-aquatic vegetation areas. Use the maximum gradient method to determine the FAI threshold for each scene.

[0081] Specifically, the FAI calculation formula is as follows:

[0082] FAI=R NIR -R′ NIR

[0083] R′ NIR =R RED +[R Swir1 -R RED ]×[(λ NIR -λ RED ) / (λ Swir1 -λ RED )]

[0084] Among them, R RED It is in the red band, with a wavelength range of 0.63-0.68 micrometers; R NIR It is in the near-infrared band, with a wavelength range of 0.845-0.885 micrometers; R Swir1 It is in the shortwave infrared band 1, with a wavelength range of 1.56-1.66 micrometers; λ RED λ is the center wavelength of the red band. NIR λ is the center wavelength of the near-infrared band. Swir1 It is the center wavelength of the shortwave infrared band 1.

[0085] The classification threshold c of the FAI image is determined using the maximum gradient method, and the steps are as follows:

[0086] 1) Generate the gradient image of the FAI image and define the gradient of a pixel as the difference between it and its neighboring pixels within a 3×3 window;

[0087] 2) Remove clean water bodies and areas with high concentrations of algal blooms, while retaining pixels in the boundary area between the two;

[0088] 3) Find the pixel with the largest gradient value among the remaining pixels, and use the FAI value corresponding to that pixel as the FAI threshold;

[0089] 4) Perform histogram statistics on the FAI threshold of multiple images, and obtain the FAI classification threshold c by subtracting twice the standard deviation from the mean of the histogram.

[0090] Histogram analysis of the FAI thresholds of 50 images showed that the results followed a normal distribution, as follows: Figure 2 As shown. Subtracting twice the standard deviation from the mean of the FAI threshold histogram yields a classification threshold c of 0.02 for FAI.

[0091] Step 5: Determine the discrimination criteria based on the thresholds of AVI, NDVI, and FAI, and build a decision tree classification model. The construction conditions are as follows:

[0092] Condition 1: When the pixels in the study area satisfy AVI≥-0.03413, they are classified as aquatic vegetation; when the pixels in the study area do not satisfy AVI≥-0.03413, they are classified as non-aquatic vegetation.

[0093] Condition 2: Among the pixels that meet Condition 1, if NDVI ≥ 0.2, it is identified as a floating-leaved / emergent plant; if NDVI < 0.2, it is identified as a submerged plant.

[0094] Condition 3: Among the pixels that do not meet condition 1, if the FAI ≥ 0.02, it is judged as an algal bloom; if the FAI < 0.02, it is judged as a water body.

[0095] Step 6: Graph the classification results in ArcGIS, as shown below. Figure 3 Based on the measured sample points and classification results, the overall classification accuracy and Kappa coefficient are calculated using the following formulas:

[0096] Overall classification accuracy:

[0097]

[0098] in, β represents the total number of samples that are correctly classified on the image, in addition to the actual samples.

[0099]

[0100] Where r is the total number of columns in the error matrix, i.e., the total number of categories; A ii It represents the number of pixels in the i-th row and j-th column of the error matrix, i.e., the number of correctly classified pixels; A i+ and A +i ...

[0101] The classification results were evaluated using the overall classification accuracy and the Kappa coefficient, as shown in Table 1.

[0102] Table 1 is an accuracy evaluation table of the automatic extraction results of aquatic vegetation and algal blooms in Taihu Lake obtained by the method described in the embodiments of the present invention. By comparing the monitoring results with the measured results, it can be seen that the monitoring accuracy of the method of the present invention for aquatic vegetation communities and algal blooms in eutrophic lakes (Taihu Lake) is higher than 80%, and the Kappa coefficient is greater than 0.75, which meets most application requirements. Figure 3 This is a schematic diagram illustrating the automatic extraction results of aquatic vegetation and algal blooms in Taihu Lake obtained by the method described in this embodiment of the invention. It can be seen that the advantages of this invention compared to existing traditional methods are: using limited measured sample points and remote sensing imagery, it can quickly and accurately obtain the distribution type and distribution area of ​​aquatic vegetation communities and algal blooms in the study area.

[0103] Table 1. Accuracy Evaluation of Extraction Results of Aquatic Vegetation and Algal Blooms in Taihu Lake

[0104]

Claims

1. An automatic extraction method of eutrophic lake aquatic vegetation and algal bloom based on Landsat images, characterized in that, The method comprises the following steps: i. Perform a tasseled cap transformation on the image data, and take the third component of the tasseled cap transformation as an aquatic vegetation index AVI to obtain an AVI image, and determine a classification threshold a for classifying aquatic vegetation and non-aquatic vegetation in the AVI image based on the following formula: wherein AVI(w) and AVI(v) represent the average AVI values of pure water body pixels and pixels in a floating leaf / emergent vegetation dense area respectively, and p is the proportion of pure water body end members; ； ii. Obtain a normalized difference vegetation index NDVI image based on the image data, which is used to distinguish floating leaf / emergent vegetation and submerged vegetation in the aquatic vegetation area, and determine a classification threshold b in the NDVI image based on the following steps: 1) generate a gradient image of the NDVI image, and define the gradient of a pixel as the difference from adjacent pixels in a 3*3 window; 2) remove submerged and floating leaf / emergent dense areas, and retain the pixels in the boundary area between the two; 3) find the pixel with the maximum gradient value in the remaining pixels, and take the NDVI value corresponding to the pixel as the NDVI threshold; 4) perform histogram statistics on the NDVI threshold of multi-scene images, and take the value obtained by subtracting twice the standard deviation from the average value of the histogram as the classification threshold b of the NDVI image; iii. Obtain a phytoplankton algae index FAI image based on the image data, which is used to distinguish algal blooms and water bodies in the non-aquatic vegetation area, and determine a classification threshold c in the FAI image based on the following steps: 1) generate a gradient image of the FAI image, and define the gradient of a pixel as the difference from adjacent pixels in a 3*3 window; 2) remove clean water bodies and high-concentration algal bloom areas, and retain the pixels in the boundary area between the two; 3) find the pixel with the maximum gradient value in the remaining pixels, and take the FAI value corresponding to the pixel as the FAI threshold; 4) perform histogram statistics on the FAI threshold of multi-scene images, and take the value obtained by subtracting twice the standard deviation from the average value of the histogram as the classification threshold c of the FAI image; Determine the discrimination conditions according to the AVI, NDVI and FAI index values and their classification thresholds a, b and c, establish a decision tree classification model, and generate an automatic extraction result image of different aquatic vegetation communities and algal blooms in the eutrophic lake. Obtain the Landsat satellite surface reflectance image data covering the study area from the GEE cloud platform.

2. The method of claim 1, wherein, Further comprising pre-processing the Landsat satellite surface reflectance image data; the pre-processing comprises image cloud removal and study area cropping; the image cloud removal is completed by using the QA band and setting the values of cloudShadowBitMask and cloudsBitMask to 0, and then masking.

3. The method of claim 1, wherein, AVI(w) is obtained based on the following steps:

4. The method of claim 1, wherein, 1) obtain the NDWI image of the image in the study area, and perform histogram statistics; 2) take the NDWI value corresponding to the second highest peak in the NDWI histogram as the NDWI threshold for extracting pure water pixels; 3) extract all the pixels with an NDWI value greater than the threshold as pure water pixels, and take the average AVI value of these pixels as the value of AVI(w). ​ 5. The method of claim 1, wherein, AVI(v) is obtained based on the following steps: 1) Select several pixels in dense floating leaf / emergent growth areas from images covering different lakes by visual interpretation combined with measured data; 2) Calculate the average AVI value of these pixels and take the value as the value of AVI(v).

6. The method of claim 1, wherein, The decision tree classification model is as follows: Based on AVI index, the study area is divided into aquatic vegetation and non-aquatic vegetation. When the pixel of the study area satisfies AVI a, it is judged as aquatic vegetation, otherwise it is judged as non-aquatic vegetation; For aquatic vegetation, the NDVI index is used for classification, when NDVI b then it is classified as floating / emergent, otherwise it is classified as submerged; For non-aquatic vegetation, the FAI index is used for classification, when FAI c, then it is classified as algal bloom, otherwise it is classified as water body.

7. The method of claim 1, wherein, The extraction accuracy of aquatic vegetation and algal blooms is evaluated by using overall accuracy and Kappa coefficient.

Images