A eutrophication index inversion method and system based on Resource-1 hyperspectral data

Abstract

The present invention provides a eutrophication index inversion method based on Resource No. 1 hyperspectral data, comprising the following steps: Step 1: Acquire remote sensing data and preprocess it, including radiometric calibration, FLAASH atmospheric correction, and image cropping; Step 2: Resample the measured spectral reflectance; Step 3: Construct a sensitive band combination; Step 4: Train a CatBoost model; Step 5: Calculate the eutrophication index. The application of this technical solution can accurately reflect COD, DIN, PO4 3‑ The multiplicative synergistic effect of the two models can be avoided to avoid the evaluation bias of single parameter or linear model.

Description

Technical Field

[0001] The present invention relates to the field of environmental remote sensing technology, and in particular to a eutrophication index inversion method and system based on Resource No. 1 hyperspectral data. Background Art

[0002] Existing technologies commonly use empirical formulas or semi-empirical models (such as band ratio methods and multivariate regression) to invert single water quality parameters. For example, using Gaofen-1 satellite data, a soft classification model is constructed based on the nonlinear relationship between band reflectance and turbidity to invert total suspended solids concentration. Alternatively, based on Sentinel-2 data, an inversion formula is established using the exponential relationship between band combinations (such as B4 / B3) and chlorophyll a concentration. These technologies artificially select features, lack generalizability to complex Class II water bodies (containing a mixture of multiple pollutants), and fail to reflect the comprehensive state of eutrophication.

[0003] Model bias caused by single parameter: only relying on chlorophyll a concentration as the eutrophication indicator, without considering chemical oxygen demand (COD), inorganic nitrogen (DIN), active phosphate (PO4 3- Because eutrophication is essentially a comprehensive effect caused by excessive nutrient input, a single parameter cannot characterize the nonlinear coupling relationship between pollutants, resulting in systematic deviations between the evaluation results and the actual eutrophication status.

[0004] Insufficient band adaptability: The Sentinel-2 spectral range (400-2200 nm) is not optimized for the high spectral characteristics of Class II coastal waters, especially for PO4 3- The spectral resolution of the sensitive short-wave infrared band (1200-1300 nm) is insufficient, which makes it impossible to effectively extract PO4 3- Concentration-dependent spectral features.

[0005] Limitations of linear models: The use of linear regression models (such as the linear relationship between NDCI and chlorophyll a) does not take into account the nonlinear interference of suspended matter and colored dissolved organic matter (CDOM) in Class II water bodies on the optical signal, resulting in inversion errors being transmitted to the eutrophication assessment results through the single-parameter model. Summary of the Invention

[0006] In view of this, the purpose of the present invention is to provide a eutrophication index inversion method and system based on Resource No. 1 hyperspectral data to accurately reflect the COD, DIN, PO4 3- The multiplicative synergistic effect of the two models can be avoided to avoid the evaluation bias of single parameter or linear model.

[0007] To achieve the above object, the present invention adopts the following technical solution: a eutrophication index inversion method based on Resource No. 1 hyperspectral data, comprising the following steps:

[0008] Step 1: Acquire remote sensing data and preprocess it, including radiometric calibration, FLAASH atmospheric correction, and image cropping.

[0009] Step 2: Resample the measured spectral reflectance;

[0010] Step 3: Construct sensitive band combination;

[0011] Step 4: Train the CatBoost model.

[0012] Step 5: Calculate the eutrophication index.

[0013] In a preferred embodiment, the radiation calibration in step 1 specifically includes the following processing steps:

[0014] For the radiometric calibration processing of the Resource No. 1 AHSI image in the present invention, you can directly click the Apply FLAASH Settings button in the Radiometric Calibration parameter setting panel of the remote sensing image processing platform ENVI. The remaining parameters will be automatically adjusted to the data input format required by FLAASH atmospheric correction. Click OK to perform radiometric calibration processing.

[0015] In a preferred embodiment, the FLAASH atmospheric correction in step 1 specifically includes the following processing steps:

[0016] In the remote sensing image processing platform ENVI, select the FLAASH atmospheric correction function. In the pop-up interface, set the satellite altitude Sensor Altitude to 778km, the regional average ground elevation to 0.023, the pixel size to 30m, the atmospheric model Atmospheric Model to Tropical, the water vapor inversion WaterRetrieval to Yes, and select the water vapor absorption band to 940nm, the Aerosol Model to Martime, the spectral smoothing to Yes, and the width to 3. After setting the parameters, click the Apply button to perform the FLAASH atmospheric correction processing.

[0017] In a preferred embodiment, the image cropping in step 1 specifically includes the following processing steps:

[0018] The image remote sensing reflectance approximate estimation method is used to remove the influence of sky light on the surface reflectance image after FLAASH atmospheric correction. The calculation formula is shown in formula (1):

[0019] (1)

[0020] Where λ is the band of the remote sensing image, that is, the original spectral value, Represents the remote sensing reflectance estimated approximately on the image, represents the surface reflectivity, represents the smaller value of the near-infrared and short-wave infrared surface reflectance, and π is the circumference of the circle. The atmospheric correction image is clipped using NDWI to identify the water area of ​​the study area. This image is used to invert the concentration of water quality parameters and is input into the constructed CatBoost model. The NDWI calculation formula is as follows:

[0021]

[0022] Where b20 is the 20th band of the Ziyuan-1 AHSI, and b41 is the 41st band of the Ziyuan-1 AHSI. The areas with pixel values ​​greater than 0.5 in the processed image are extracted to obtain the water area of ​​the study area.

[0023] In a preferred embodiment, in step 2, the inversion of water quality parameter concentrations also requires real water body reflectance; using the real water body reflectance as a medium, a nonlinear relationship model between satellite imagery and measured water quality parameter concentrations is established by resampling the real water body reflectance to the Resource-1 AHSI remote sensing reflectance; using the 1 nm resolution spectral curve measured by the ASD spectrometer, the measured spectral reflectance is resampled to the 10 nm spectral interval of the Resource-1 AHSI using a linear interpolation method;

[0024] That is, for each sample j, the original spectrum is defined as: λ=[λ1,λ2,...,λ n ], corresponding to the spectral value S j (λ)=[S j1 ,S j2 ,...,S jn ], define the target wavelength sequence as λ * =[λ1 * ,λ2 * ,...,λ n * ], then after interpolation, at each The spectral value at is:

[0025]

[0026] Where, is the target interpolation wavelength, The jth sample after interpolation at the target interpolation wavelength The estimated spectral value on 、 The closest Two adjacent wavelengths, 、 The jth sample at wavelength 、 Finally, the chemical oxygen demand concentration COD, inorganic nitrogen concentration DIN, active phosphate concentration PO4 3- Take the logarithm with base 10 for processing, i.e. lg(COD), lg(DIN), lg(PO4 3- ).

[0027] In a preferred embodiment, in step 3, a combination of ratio, normalization, three-band and ratio, three-band inverse difference product satellite band is constructed, and the Pearson correlation coefficient is used to calculate the satellite band. Calculate the band combination and lg(COD), lg(DIN), lg(PO4 3- ) correlation coefficient; select combinations with an absolute value of correlation coefficient > 0.6 and exclude combinations with a band interval < 20 nm:

[0028] (2)

[0029] Where, is the covariance, is the product of the standard deviations of the two variables; the detailed calculation formula is:

[0030]

[0031] Among them, X is the combined result of the reflectivity of each band, R COD 、R DIN 、R PO4 3- Represents COD, DIN, PO4 respectively 3- Correlation with the combined results of the band reflectances.

[0032] In a preferred embodiment, in step 4, the resampled AHSI reflectance and its corresponding water quality parameters are divided into a training set and a test set in a ratio of 8:2; the CatBoost library is imported into Python, and the selected parameters are used to invert lg(COD), lg(DIN), lg(PO4 3- )concentration;

[0033] Suppose there are two hyperparameters Parme1 and Parme2, and the parameter values ​​corresponding to the two hyperparameters are a1, a2 and b1, b2, b3, namely:

[0034]

[0035] Traverse all combinations , all parameter combinations are used to train and analyze the model, and finally the model parameters are selected.

[0036] In a preferred embodiment, in step 5, the CatBoost inversion model obtained in step 4 is used to input the atmospherically corrected AHSI remote sensing image, and the lg(COD), lg(DIN), lg(PO4 3- ) value, perform mathematical operations on it, convert it into a logarithmic value, and use the eutrophication index calculation formula as shown in formula (3):

[0037] (3)

[0038] Where, E is the eutrophication index, dimensionless; COD is the chemical oxygen demand concentration obtained by inversion; PO4 3- is the active phosphate concentration obtained by inversion; DIN is the inorganic nitrogen concentration obtained by inversion.

[0039] The present invention also provides a eutrophication index inversion system based on Resource No. 1 hyperspectral data, which runs the above-mentioned eutrophication index inversion method based on Resource No. 1 hyperspectral data.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] 1. Realize multi-parameter coordinated nonlinear eutrophication assessment

[0042] Accurately reflect COD, DIN, PO4 3- The multiplicative synergistic effect of the two models can be avoided to avoid the evaluation bias of single parameter or linear model.

[0043] 2. Accurate feature extraction of hyperspectral band combinations

[0044] Based on the 330-band data of Resource No. 1 AHSI, four types of combinations are constructed: band ratio (R), normalization (ND), sum ratio (SR), and reciprocal difference product (RDP). 3- Improve the inversion accuracy and reduce spectral noise interference.

[0045] 3. Dynamic monitoring of eutrophication with high spectral resolution

[0046] Taking advantage of the wide width (60 km) and high spectral resolution (10 nm) of the Ziyuan-1 AHSI image, the entire study area is covered. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram of a preferred embodiment of the present invention. DETAILED DESCRIPTION

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

[0049] It should be noted that the following detailed descriptions are illustrative and intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the art to which the present application belongs.

[0050] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application; as used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form, and it should be understood that when the terms "comprise" and / or "include" are used in this specification, they indicate the presence of features, steps, operations, devices, components and / or their combinations.

[0051] A eutrophication index inversion method based on Resource-1 hyperspectral data, referring to Figure 1 , including the following steps:

[0052] Step 1: Remote sensing data acquisition and preprocessing

[0053] Because satellite signals are susceptible to atmospheric influences during transmission, atmospheric correction is required to mitigate these effects and convert the pixel values ​​of satellite images into true water reflectance. This method converts the raw pixel values ​​of the Resource-1 AHSI imagery into radiance values. First, radiometric calibration is used to convert the raw reflectance recorded by the sensor into atmospheric surface reflectance. The FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) algorithm is then used to eliminate interference from aerosols, water vapor, and other factors, converting the atmospheric surface reflectance into true surface reflectance.

[0054] For the radiometric calibration of the Resource-1 AHSI image in this invention, you can directly click the Apply FLAASH Settings button in the Radiometric Calibration parameter setting panel in the ENVI software. The remaining parameters will be automatically adjusted to the data input format required by FLAASH atmospheric correction. Click OK to perform the radiometric calibration.

[0055] For the FLAASH atmospheric correction of the Resource-1 AHSI imagery described in this paper, set the SensorAltitude to 778 km, the Regional Ground Elevation to 0.023, the Pixel Size to 30 m, the Atmospheric Model to Tropical, the Water Retrieval to Yes, the Water Vapor Absorption Band to 940 nm, the Aerosol Model to Martime, the Spectral Polishing to Yes, and the Width to 3. After setting the parameters, click the Apply button to perform the FLAASH atmospheric correction.

[0056] The surface reflectance image after FLAASH atmospheric correction needs to remove the influence of sky light. For the Resource-1 AHSI image, this paper adopts an approximate estimation method based on image remote sensing reflectance, and the calculation formula is shown in formula (1):

[0057] (1)

[0058] Where, Represents the remote sensing reflectance estimated approximately on the image, represents the surface reflectivity, Indicates the smaller value of near infrared and short wave infrared surface reflectivity. Set to the average value of each band in 730-760nm, The NDWI is set as the average reflectance of each band in the 1530-1630nm range. The atmospheric correction image is used to crop the water area of ​​the study area using NDWI. This image is used to invert the concentration of water quality parameters and is input into the constructed CatBoost model. The NDWI calculation formula is as follows:

[0059]

[0060] Where b20 is the 20th band of the Ziyuan-1 AHSI, and b41 is the 41st band of the Ziyuan-1 AHSI. The waters of the study area are obtained by extracting the areas with pixel values ​​greater than 0.5 from the processed image.

[0061] Step 2: Resampling of measured data

[0062] Inverting water quality parameter concentrations requires not only satellite imagery but also real water reflectance. Using real water reflectance as a medium, resampling real water reflectance to the remote sensing reflectance of the Resource-1 AHSI (AHSI) allows a nonlinear relationship model to be established between satellite imagery and measured water quality parameter concentrations. This method uses a 1nm resolution spectral curve measured by an ASD spectrometer and resamples the measured spectral reflectance to the 10nm spectral interval of the Resource-1 AHSI (covering the 400-2500nm range) using linear interpolation.

[0063] That is, for each sample j, the original spectrum is defined as: λ=[λ1,λ2,...,λ n ], corresponding to the spectral value S j (λ)=[S j1 ,S j2 ,...,S jn ], define the target wavelength (resource number one hyperspectral AHSI) sequence as λ * =[λ1 * ,λ2 * ,...,λ n * ], then after interpolation, at each The spectral value at is:

[0064]

[0065] Where, is the target interpolation wavelength, is the jth sample at the target wavelength after interpolation The estimated spectral value on 、 For the closest Two adjacent wavelengths, 、 is the jth sample at wavelength 、 Finally, due to the measured water quality parameters COD, DIN, PO4 3- The true value range is small, and the three elements COD, DIN, and PO4 need to be 3- Take the logarithm with base 10 for processing, i.e. lg(COD), lg(DIN), lg(PO4 3- ) to ensure the accuracy of the modeling results.

[0066] Step 3: Construction of sensitive band combination

[0067] The accuracy of the model inversion can be improved by constructing four band combinations (ratio, normalization, three-band sum ratio, and three-band reciprocal difference product) from different satellite bands through conventional mathematical operations. The present invention calculates the correlation coefficients of all possible band combinations with lg(COD), lg(DIN), lg(PO4) based on the Pearson correlation coefficient (Formula 2). 3- ). Combinations with an absolute value of correlation coefficient > 0.6 were selected, and combinations with a band spacing < 20 nm were excluded to avoid spectral overlap leading to multicollinearity.

[0068] (2)

[0069] Where, is the covariance, is the product of the standard deviations of the two variables. For the present invention, the detailed formula for calculation is:

[0070]

[0071] Here, X is the combined result of the reflectivity of each band. The present invention ultimately screened out four band combinations: 567nm / 533nm, (585nm-473nm) / (585nm+473nm), (438nm+516nm) / 473nm, and (1 / 670nm-1 / 695nm)×816nm as characteristic variables for inversion.

[0072] Step 4: CatBoost model training analysis

[0073] The resampled AHSI reflectance and its corresponding water quality parameters were divided into a training set (39 samples) and a test set (10 samples) in a ratio of 8:2. The CatBoost library was imported into Python and the optimal parameters were determined using a network search method (e.g. As shown in the figure), the CatBoost model with optimal parameters is used to invert lg(COD), lg(DIN), lg(PO4 3- )concentration.

[0074] surface CatBoost optimal parameter settings

[0075] parameter value iterations 600 learning_rate 0.01 l2_leaf_reg 2 max_depth 3 loss_function RMSE eval_metric R2 verbose 0

[0076] Grid Search is an exhaustive search algorithm that tries all possible combinations in a predefined parameter space to find the parameter configuration that performs best on the validation set. For example, suppose there are two hyperparameters, Parme1 and Parme2, and the corresponding configurable values ​​are a1, a2 and b1, b2, and b3, respectively. That is:

[0077]

[0078] The network search will traverse all combinations The network search method will conduct model training analysis on all parameter combinations and finally select the parameter combination with the highest accuracy in the training set and test set as the optimal parameter of the model.

[0079] Step 5: Calculation of eutrophication index

[0080] According to the CatBoost inversion model obtained in the above steps, the three water quality parameters lg(COD), lg(DIN), lg(PO4 3- ) value, perform mathematical operations on it, convert it into a logarithmic value, and use the eutrophication index calculation formula as shown in formula (3):

[0081] (3)

[0082] Where, E is the eutrophication index, dimensionless; COD is the chemical oxygen demand concentration obtained by inversion (mg / L); PO4 3- is the active phosphate concentration obtained by inversion (mg / L); DIN is the inorganic nitrogen concentration obtained by inversion (mg / L).

[0083] The key points and protection points of this application proposal

[0084] 1. Sensitive band screening method for hyperspectral data

[0085] Based on the 330 hyperspectral bands of Resource No. 1 AHSI, a feature set was constructed through the following four combinations, and the band combination with the highest correlation coefficient and no redundancy was selected ( ):

[0086] Table 2 Characteristic band reflectivity combinations and corresponding formulas

[0087]

[0088] Existing techniques often use single-band ratios (such as NDCI) or raw band reflectivity, without optimizing the combination types required for multi-parameter inversion. This method addresses the feature redundancy caused by blind band selection in traditional methods by combining bands (such as the reciprocal difference product) and using a Pearson correlation coefficient threshold (|r|>0.6).

[0089] By utilizing the high spectral resolution (330 bands) of the Resource-1 AHSI image and combining it with the inversion process, the E value was calculated pixel by pixel and a spatial distribution map of the eutrophication index was generated.

[0090] Existing technologies (such as drone hyperspectral) have a small coverage range (<10 km²) and cannot meet the large-scale monitoring needs of coastal aquaculture areas. This invention solves the problem of Sentinel-2 and other wide-band satellites covering short-wave infrared (PO4 3- The problem of insufficient resolution of sensitive areas.

Claims

1. A eutrophication index inversion method based on Resource No. 1 hyperspectral data, characterized in that: The following steps are involved: Step 1: Acquire remote sensing data and preprocess it, including radiometric calibration, FLAASH atmospheric correction, and image cropping. Step 2: Resample the measured spectral reflectance; Step 3: Construct sensitive band combination; Step 4: Train the CatBoost model. Step 5: Calculate the eutrophication index; In step 2, the inversion of water quality parameter concentrations also requires real water body reflectance; using the real water body reflectance as a medium, a nonlinear relationship model between satellite images and measured water quality parameter concentrations is established by resampling the real water body reflectance to the Resource-1 AHSI remote sensing reflectance; using the 1 nm resolution spectral curve measured by the ASD spectrometer, the measured spectral reflectance is resampled to the 10 nm spectral interval of the Resource-1 AHSI using a linear interpolation method; That is, for each sample j, the original spectrum is defined as: λ=[λ1,λ2,...,λ n ], corresponding to the spectral value S j (λ)=[S j1 ,S j2 ,...,S jn ], define the target wavelength sequence as λ * =[λ1 * ,λ2 * ,...,λ n * ], then after interpolation, at each The spectral value at is: Where, is the target interpolation wavelength, The jth sample after interpolation at the target interpolation wavelength The estimated spectral value on 、 The closest Two adjacent wavelengths, 、 The jth sample at wavelength 、 Finally, the chemical oxygen demand concentration COD, inorganic nitrogen concentration DIN, active phosphate concentration PO4 3- Take the logarithm with base 10 for processing, i.e. lg(COD), lg(DIN), lg(PO4 3- ); In step 3, the ratio, normalization, three-band and ratio, three-band inverse difference satellite band combination is constructed, and the Pearson correlation coefficient is used to Calculate the band combination and lg(COD), lg(DIN), lg(PO4 3- ) correlation coefficient; select combinations with an absolute value of correlation coefficient > 0.6 and exclude combinations with a band interval < 20 nm: （2） Where, is the covariance, is the product of the standard deviations of the two variables; the detailed calculation formula is: Among them, X is the combined result of the reflectivity of each band, R COD 、R DIN 、R PO4 3- Represents COD, DIN, PO4 respectively 3- Correlation with the combined results of the band reflectances.

2. The eutrophication index inversion method based on Resource No. 1 hyperspectral data according to claim 1 is characterized in that: The radiation calibration in step 1 specifically includes the following processing steps: For the radiometric calibration of the Resource-1 AHSI image, directly click the Apply FLAASH Settings button in the Radiometric Calibration parameter setting panel of the remote sensing image processing platform ENVI. The remaining parameters will be automatically adjusted to the data input format required by FLAASH atmospheric correction. Click OK to perform the radiometric calibration.

3. The eutrophication index inversion method based on Resource No. 1 hyperspectral data according to claim 1 is characterized in that: The FLAASH atmospheric correction in step 1 specifically includes the following processing steps: In the remote sensing image processing platform ENVI, select the FLAASH atmospheric correction function. In the pop-up interface, set the satellite altitude Sensor Altitude to 778km, the regional average ground elevation to 0.023, the pixel size to 30m, the atmospheric model Atmospheric Model to Tropical, the water vapor inversion WaterRetrieval to Yes, and select the water vapor absorption band to 940nm, the Aerosol Model to Martime, the spectral smoothing to Yes, and the width to 3. After setting the parameters, click the Apply button to perform the FLAASH atmospheric correction process.

4. The eutrophication index inversion method based on Resource No. 1 hyperspectral data according to claim 1 is characterized in that: The image cropping in step 1 specifically includes the following processing steps: The image remote sensing reflectance approximate estimation method is used to remove the influence of sky light on the surface reflectance image after FLAASH atmospheric correction. The calculation formula is shown in formula (1): （1） Where λ is the band of the remote sensing image, that is, the original spectral value, Represents the remote sensing reflectance estimated approximately on the image, represents the surface reflectivity, represents the smaller value of the near-infrared and short-wave infrared surface reflectance, and π is the circumference of the circle. The atmospheric correction image is clipped using NDWI to identify the water area of ​​the study area. This image is used to invert the concentration of water quality parameters and is input into the constructed CatBoost model. The NDWI calculation formula is as follows: Where b20 is the 20th band of the Ziyuan-1 AHSI, and b41 is the 41st band of the Ziyuan-1 AHSI. The areas with pixel values ​​greater than 0.5 in the processed image are extracted to obtain the water area of ​​the study area.

5. The eutrophication index inversion method based on Resource No. 1 hyperspectral data according to claim 1 is characterized in that: In step 4, the resampled AHSI reflectance and its corresponding water quality parameters were divided into training set and test set in a ratio of 8:2; the CatBoost library was imported into Python, and the selected parameters were used to invert lg(COD), lg(DIN), lg(PO4 3- )concentration; Suppose there are two hyperparameters Parme1 and Parme2, and the parameter values ​​corresponding to the two hyperparameters are a1, a2 and b1, b2, b3, namely: Traverse all combinations , all parameter combinations are used to train and analyze the model, and finally the model parameters are selected.

6. The eutrophication index inversion method based on Resource No. 1 hyperspectral data according to claim 1 is characterized in that: In step 5, the CatBoost inversion model obtained in step 4 is used to input the atmospherically corrected AHSI remote sensing image, and the lg(COD), lg(DIN), lg(PO4 3- ) value, perform mathematical operations on it, convert it into a logarithmic value, and use the eutrophication index calculation formula as shown in formula (3): （3） Where, E is the eutrophication index, dimensionless; COD is the chemical oxygen demand concentration obtained by inversion; PO4 3- is the active phosphate concentration obtained by inversion; DIN is the inorganic nitrogen concentration obtained by inversion.

7. A eutrophication index inversion system based on Resource No. 1 hyperspectral data, characterized in that: Run the eutrophication index inversion method based on Resource No. 1 hyperspectral data as described in any one of claims 1 to 6.

Images