Inland water turbidity satellite remote sensing method based on optical classification and spectrum simulation

Abstract

The invention discloses an inland water body turbidity satellite remote sensing method based on optical classification and spectral simulation, which comprises the following steps: S1, acquiring water body spectral data, water body turbidity, satellite remote sensing images and meteorological data of a to-be-detected area, and calculating water body remote sensing reflectivity; s2, preprocessing the satellite remote sensing image to obtain a water body area; performing equivalent calculation on the remote sensing reflectivity data of the water body actual measurement spectrum to obtain the remote sensing reflectivity of each channel of the satellite sensor; performing water body optical classification by combining the actually measured water body spectral shape and turbidity distribution characteristics of the research area; s3, performing correlation analysis on different optical types of the actually measured spectrum according to the actually measured turbidity and the remote sensing reflectivity of multiple channels of the satellite, selecting channel remote sensing reflectivity data with high turbidity correlation as a sample, and constructing a turbidity remote sensing model; the turbidity remote sensing model is trained and tested, and the optimal turbidity remote sensing model is screened out. The method can improve the precision of remote sensing inversion of the turbidity of the inland water body.

Description

Technical Field

[0001] The present invention relates to the technical field of spectral analysis, and in particular to a satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation. Background Art

[0002] Turbidity is a key indicator of water quality. It is not only related to the concentration of suspended particulate matter in the water, but also influenced by the composition, size, shape, and surface reflectivity of these particles. The migration, sedimentation, and resuspension of suspended matter in water all affect turbidity. Turbidity is an optical effect that reflects the degree to which suspended matter in water obstructs light from entering the water due to absorption and scattering. Turbidity determines the depth of light penetration into the water, affects solar absorption and heat distribution, alters the vertical temperature structure of the water, and thus affects the growth of aquatic plants. Therefore, turbidity is of great significance in the study of aquatic ecosystems, hydrodynamic conditions, changes in the water environment, and the transport and circulation of substances.

[0003] Conventional water turbidity monitoring uses turbidity meters from underway ships or continuous observations from fixed buoy platforms. While this method offers high measurement accuracy, it suffers from limited observation sites and insufficient spatial representation of the data. This is particularly true for inland waters, where water quality fluctuates frequently. Conventional methods are unable to rapidly, efficiently, and comprehensively monitor, manage, and predict water quality parameters like turbidity on a large scale.

[0004] Satellite remote sensing technology provides high-temporal and spatial resolution data for water turbidity monitoring, which is unattainable by conventional methods. It is not limited to a single point in the water body's space and offers advantages such as speed, periodicity, and wide coverage. Currently, water turbidity remote sensing primarily relies on visible-to-near-infrared satellite sensor data, atmospherically corrected water surface reflectance, and empirical statistical models based on field turbidity data. This approach rarely considers the optical properties of water and lacks analysis of the relationship between the spectral characteristics of each satellite sensor channel and turbidity. Consequently, the resulting turbidity remote sensing models have limited generalizability. Furthermore, due to the low spectral reflectance of water and the fact that over 90% of satellite remote sensing imagery is noise, water-leaving radiance, a weak signal, is a key parameter in water remote sensing. Existing radiation transfer models are unable to accurately extract water-leaving radiance when atmospherically correcting remote sensing images. Water remote sensing primarily utilizes optical satellite imagery in the visible-to-near-infrared bands, but current approaches rarely consider the optical properties of water, resulting in poor applicability and generalizability of the resulting remote sensing models. Summary of the Invention

[0005] Purpose of the invention: The purpose of the present invention is to provide a satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation, to construct the best model for water turbidity remote sensing for different optical types, and to improve the accuracy of inversion of inland water turbidity remote sensing.

[0006] Technical solution: A satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation, including the following steps: S1, obtain water spectral data, water turbidity, satellite remote sensing images and meteorological data of the area to be measured; calculate the water remote sensing reflectance based on the acquired water spectral data; S2, pre-process the satellite remote sensing image to obtain the water body area; the remote sensing reflectance data R of the measured spectrum of the water body rs (λ) performs equivalent calculations to obtain the remote sensing reflectivity R of each channel of the satellite sensor rs (ch); cluster the measured spectra of water bodies, and perform optical classification of water bodies based on the measured spectral shapes and turbidity distribution characteristics of the water bodies in the study area; S3, for optical types with different measured spectra, classify them according to the measured turbidity and the equivalent calculated satellite multi-channel remote sensing reflectance, perform correlation analysis on all turbidities of the same type and the remote sensing reflectance of a single channel or a combination of multiple channels, select the channel remote sensing reflectance data with high correlation with turbidity as samples, use turbidity as the output parameter, and construct a turbidity remote sensing model; divide the samples into training and test sets, train and test the turbidity remote sensing model, and screen the best turbidity remote sensing model.

[0007] Furthermore, in step S2, the preprocessing of the satellite remote sensing image includes radiometric correction, geometric correction, atmospheric correction and water area extraction of the image; The specific steps for atmospheric correction of remote sensing images are as follows: S31, calculate the atmospheric Rayleigh scattering optical depth and aerosol scattering optical depth using measured atmospheric pressure and visibility data; calculate the ozone absorption optical depth based on ozone concentration; S32, combining the observation geometry data of remote sensing images, calculates the atmospheric direct transmittance and diffuse transmittance, atmospheric Rayleigh scattering radiance and aerosol scattering radiance of each channel of satellite images; S33, calculate the water surface foam reflection radiance and solar direct flare radiance based on the water surface wind speed data and observation geometry data, and calculate the water-leaving radiance based on the total radiance of the image; S34, combined with the solar irradiance of the top atmosphere of each channel, calculates the water surface remote sensing reflectivity of each channel and completes the atmospheric correction of the satellite images of each channel.

[0008] Furthermore, the normalized difference water index (NDWI) or the improved normalized difference water index (MNDWI) is constructed using the remote sensing reflectance of the green band and the near-infrared band or the short-wave infrared band after atmospheric correction. The water area and the non-water area are distinguished from the atmospheric corrected image by the threshold method. , , Where R rs-G 、R rs-NIR and R rs-SWIR Represent the remote sensing reflectance of the green light band, near infrared band and shortwave infrared band respectively.

[0009] Furthermore, in step S2, the spectral response function f(λ) of each channel of the satellite sensor is combined to calculate the water surface remote sensing reflectance data R of the measured spectrum of the water body. rs (λ) is calculated equivalently to obtain the water surface remote sensing reflectivity R of each channel of the satellite sensor rs (ch), the calculation formula is:

[0010] Where ch represents the central wavelength of each channel of the satellite sensor, λ max ,λ min They represent the maximum wavelength and minimum wavelength of the spectral response of each channel of the sensor.

[0011] Furthermore, in step S3, all samples of the same type are divided into training sets and test sets in a ratio of 2:1, and regression analysis, neural network, random forest and XGBoost algorithms are used to establish turbidity remote sensing models respectively; three evaluation indicators, namely correlation coefficient R, root mean square error RMSE and mean absolute error MAE, are used to screen out the best turbidity remote sensing model according to the test results.

[0012] Compared with the prior art, the present invention has the following significant effects: 1. The present invention uses meteorological data from satellite transits to accurately calculate atmospheric radiation, water surface foam radiation, and direct solar radiation on the water surface. This method solves the problem of inaccurate calculation of atmospheric radiation by several standard atmospheric models and aerosol models in current radiation transfer models, and compensates for the defect that current radiation transfer models cannot calculate water surface foam radiation and direct solar radiation on the water surface. 2. This invention constructs the optimal model for water turbidity remote sensing based on different optical types by optically classifying water spectra and satellite images, which has higher inversion accuracy and universality. 3. The present invention uses the measured water body spectrum and combines it with the spectral response function of each channel of the satellite sensor to equivalently calculate the water surface remote sensing reflectance of each sensor channel, thereby establishing a turbidity remote sensing model, overcoming the model error problem caused by inaccurate atmospheric correction of satellite remote sensing images. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 is a flow chart of the present invention; Figure 2This is a comparison chart of the atmospheric correction result of the OLCI image and the equivalent calculated remote sensing reflectance of the measured spectrum in an embodiment of the present invention; Figure 3 : The results of water body spectrum classification in the embodiment of the present invention are shown in Figure 1, (a) is type I, (b) is type II, and (c) is type III; (d) is the typical spectra of the three optical types; Figure 4 The OLCI image classification results of typical lakes in the embodiment of the present invention are shown in Figure 1, where (a) is Luoma Lake, (b) is Hongze Lake, and (c) is Gaoyou Lake. Figure 5 Schematic diagram of the optimal regression model and verification scatter plot of turbidity remote sensing for different water body types in an embodiment of the present invention, wherein (a) is the model of type I, (b) is the verification of the model of type I, (c) is the model of type II, (d) is the verification of the model of type II, (e) is the model of type III, and (f) is the verification of the model of type III; Figure 6 The scatter plots of the training results and test results of the random forest turbidity model for different water body types in the embodiment of the present invention are shown, where (a) is type I, (b) is type II, and (c) is type III. Figure 7 The scatter plots of the training and testing results of the XGBoost turbidity model for different water body types are shown below, where (a) is type I, (b) is type II, and (c) is type III. Figure 8 These are the spatial distribution maps of turbidity in three typical lakes in the summer of 2020, among which (a) is Luoma Lake, (b) is Hongze Lake, and (c) is Gaoyou Lake. DETAILED DESCRIPTION

[0014] The present invention will be described in further detail below with reference to the accompanying drawings and specific implementations.

[0015] To address the challenges of existing technologies, this paper proposes a method for atmospheric correction of remote sensing images using meteorological data from simultaneous satellite observations. This method utilizes measured hyperspectral reflectance data from inland water bodies, combined with the spectral response functions of each satellite sensor channel, to simulate the water surface reflectance of each channel. The K-means method is then used to perform optical classification of the hyperspectral reflectance data and satellite imagery. Statistical regression and machine learning methods are then used to construct optimal models for remote sensing turbidity of water bodies of different optical types. The turbidity of inland water bodies is then inverted from satellite remote sensing imagery.

[0016] The present invention provides a satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation. To verify the proposed method, three lakes in the Huaihe River Basin (Luoma Lake, Hongze Lake, and Gaoyou Lake) were used as the study area. The measured spectra of the lakes were optically classified and the water surface remote sensing reflectance of the Sentinel-3 satellite Ocean and Land Colorimeter (OLCI) was equivalently calculated. Turbidity remote sensing inversion models for different water optical types were constructed, and the spatial distribution of turbidity in the three lakes was inverted. The technical solution implementation process is as follows: Figure 1 As shown, the following steps are included: Step 1: Acquisition of sample data; The sample data for the Huaihe River Basin mainly includes water spectral data, water turbidity, satellite remote sensing images and meteorological data. The implementation steps are as follows: Step 11, measuring water body spectrum data; Water surface remote sensing reflectance data is the main variable of turbidity remote sensing. rs The water surface spectrum measurement method can be used to measure the water surface spectrum radiance L sw , sky spectral radiance L sky and the spectral radiance L of the gray reference plate p , calculate the water-leaving radiance L w and the total solar irradiance E on the water surface d (0 + ), and then calculated by the following formula: (1) Usually the total spectral radiance of the water surface is L sw The water-leaving radiance L contains the water component information w , the radiance L of sky light reflected on the water surface sky , radiance of water foam L wc , and the radiance L of the water surface reflecting direct sunlight g , the calculation formula is as follows: (2) Where: r is the reflectivity of the air-water interface, generally ranging from 2.1 to 5.0%, and is related to the sun's position. , observation geometry , wind speed and direction or water surface roughness and other factors.

[0017] In order to avoid the influence of direct sunlight and ship shadows, the angle between the spectrometer observation plane and the solar incident plane (azimuth) should be kept constant during the water surface radiance measurement. ) is [90°, 135°], the angle between the spectrometer and the normal direction of the water surface (zenith angle ) into [30°, 45°]; when measuring sky radiance, the spectrometer observation azimuth angle is maintained 135°, zenith angle When avoiding water foam and direct sunlight, the water surface spectral radiance L sw When measuring, L wc and L g are all 0, the water-leaving radiance L w It can be expressed as: (3) For calm water surface, the reflectivity r of the air-water interface can be taken as 2.2%; when the wind speed on the water surface is about 5 m / s, r can be taken as 2.5%.

[0018] Total solar irradiance on the water surface E d (0 + ) can use a spectrometer to vertically observe the specific reflectivity ρ under sunlight p The radiance L of the diffuse gray reference plate p Calculation is performed, the reflectivity ρ of the reference plate p It is usually around 0.3 and is strictly calibrated before leaving the factory. When measuring the spectral radiance of the reference plate, the height of the spectrometer probe is consistent with that of the water surface. Therefore, the total solar irradiance E on the water surface is d (0 + ) can be calculated as: (4) In this embodiment, field sampling was carried out in the main lakes and reservoirs in the Huaihe River Basin, and the water surface spectrum radiance L was measured using a ground object spectrometer and a water surface spectrum measurement method. sw , sky spectral radiance L sky and the spectral radiance L of the gray reference plate p , and further calculate the water surface remote sensing reflectivity. The number of samples collected in this embodiment is 181.

[0019] Step 12: Water turbidity measurement Water turbidity can be measured with a turbidity meter. Place the meter on the surface of the water and save the data after the instrument reading stabilizes. Water turbidity measurement is performed simultaneously with water spectrum measurement.

[0020] In this embodiment, while collecting the spectrum of the water body, the turbidity of the water body is measured on-site using a US HACH turbidimeter 1900C.

[0021] Step 13: Acquisition of satellite remote sensing images Satellite images used for remote sensing of inland water turbidity generally select clear-sky land resource satellite data with medium to high spatial resolution, such as the Landsat series, Sentinel-2 / 3 series, and Gaofen GF series. The sensors cover data from multiple spectral channels from visible to near-infrared.

[0022] This example uses the medium-resolution multispectral Sentinel-3 Ocean and Land Colorimeter (OLCI) imagery. Sentinel-3A and Sentinel-3B form a network, increasing the temporal resolution to two days. The OLCI sensor is equipped with 21 bands, with parameters for each band shown in Table 1. Bands 2-12 are specifically designed for water color remote sensing, offering a high signal-to-noise ratio, enabling more accurate capture of the spectral characteristics of water bodies. The remaining bands are for oxygen and water vapor absorption, as well as atmospheric correction. The spatial resolution is 300 meters. The satellite imagery covers the entire Huaihe River basin and can be used to infer turbidity in major lakes. Image data can be downloaded directly from the European Space Agency's website.

[0023] Table 1 OLCI band parameters

[0024] Step 14, meteorological data acquisition; Meteorological data are mainly used for atmospheric correction of satellite remote sensing images, mainly including atmospheric pressure, visibility and wind speed. The data can be obtained from meteorological stations in various places.

[0025] Although OLCI images include meteorological data from the European Centre for Medium-Range Weather Forecasts (ECMWF) model reanalysis, their accuracy is inferior to that of measured meteorological data. There are multiple meteorological observation stations distributed across the Huaihe River Basin, including Suqian Station (58131) near Luoma Lake, Hongze Station (58139), Xuyi Station (58138), Siyang Station (58132), and Sihong Station (58135) around Hongze Lake, and Gaoyou Station (58241) near Gaoyou Lake. This example uses atmospheric pressure, visibility, and wind speed data from these stations to calibrate the Sentinel-3 OLCI remote sensing imagery.

[0026] Step 2, processing of sample data; Sample data processing mainly includes preprocessing of satellite remote sensing images, equivalent calculation of OLCI water surface remote sensing reflectance based on measured spectra, and optical classification of water bodies.

[0027] Step 21, preprocessing of OLCI remote sensing images; The preprocessing of OLCI remote sensing images mainly includes radiation correction, geometric correction, atmospheric correction and water area extraction, and is completed with the help of remote sensing image processing software ENVI / IDL programming.

[0028] Step 211, radiometric correction of satellite images; The radiation correction of OLCI remote sensing images converts the DN value of each channel image into radiance L according to a linear relationship by combining its attribute scale factor (scale_factor) and offset (add_offset). The calculation formula is: L = a + b·DN (5) Where a and b are calibration coefficients; DN represents the original digital signal value of the image.

[0029] Step 212, geometric correction of satellite images; Since the original satellite image does not have geometric coordinates, a map or image of the same area with a high-precision coordinate system can be used as a reference for geometric correction of the image. The coordinates of several ground objects with the same name are selected from the image to be corrected and the reference image, and the row and column numbers of the pixel of the image to be corrected and the geographic coordinate control point information are established. By establishing an image geometric correction model and performing pixel resampling, the geometric correction of the original satellite image can be achieved.

[0030] When geometrically correcting OLCI remote sensing images, since the row and column numbers of the image are consistent with those of the built-in geographic coordinate data (geo_coordinates.nc), a set of latitude and longitude data is read every 50 rows and 50 columns, starting from row 1 and column 1, to form the control points (GCPs) for image geometric correction. This embodiment uses a polynomial model to establish the relationship between the longitude and latitude of image pixels and their row and column numbers, and uses the nearest neighbor interpolation method to calculate the pixel values to complete the geometric correction of OLCI remote sensing images.

[0031] Step 213, atmospheric correction of satellite images; In inland water remote sensing, the water's low spectral reflectivity means that the water-off-water radiance received by satellite sensors is a weak signal. Over 90% of the radiation in remotely sensed images comes from atmospheric pathlength radiation. Therefore, this paper proposes using meteorological data synchronously observed during satellite transits to perform atmospheric correction on remotely sensed water images.

[0032] According to the atmospheric radiation transfer theory, the total radiance L received by the satellite sensor above the water body is t (λ) can be divided into the following 5 parts: atmospheric Rayleigh scattering radiance L r (λ), aerosol scattering radiance L a (λ), radiance of direct solar flare on water surface L g (λ), radiance of water foam reflection L wc (λ) and the water-leaving radiance L containing water component information w (λ). It can be expressed as: (6) Where: λ is the center wavelength of each channel of the sensor; T (λ) and t(λ) are the direct transmittance and diffuse transmittance of the atmosphere respectively, which can be calculated based on the atmospheric Rayleigh scattering optical thickness. , aerosol scattering optical depth and ozone absorption optical depth , and combined with the sensor zenith angle θ vcalculate: (7) (8) Rayleigh optical depth of the atmosphere It can be calculated directly from the surface atmospheric pressure: (9) Where P is the local atmospheric pressure and P0 is the sea level pressure, 1013.25 pha.

[0033] Aerosol Optical Depth Meteorological visibility can be calculated using: (10) Where v is the visibility at sea level, z is the altitude, and H1=0.886+0.0222v(km).

[0034] Ozone optical depth It can be calculated according to the following formula: or (11) Where: is the ozone absorption coefficient, in cm -1 DU is the ozone concentration, the unit is Dobson, H oz is the characteristic thickness of ozone.

[0035] Atmospheric Rayleigh scattering radiation brightness L r (λ) and aerosol scattering brightness L a (λ) is collectively referred to as the atmospheric path radiance, which can be calculated according to the following formula: (12) Where: i=r represents atmospheric Rayleigh scattering, i=a represents aerosol scattering, is the single scattering albedo, p i is the scattering factor, is the instantaneous solar irradiance at the top of the atmosphere after two ozone decays , The calculation formula is: (13) Scattering factor The scattering phase function P i And the water surface diffuse reflectivity ρ is calculated: (14) Where: 、 are the water surface diffuse reflectance in the satellite observation direction and the sun incidence direction, 、 are the forward and backward scattering angles of the atmosphere, respectively, and are the forward and backward scattering phase functions of the atmosphere, respectively.

[0036] The diffuse reflectance of the water surface in the satellite observation direction and the solar incident direction is calculated according to the Fresnel reflection formula (15): (15) In the formula: and are respectively expressed as the incident angle and the refraction angle, , n w is the refractive index of water.

[0037] The atmospheric scattering angle can be calculated according to the solar geometric position and the sensor observation geometric position as follows: (16) In the formula: and are the zenith angles of the sun and the satellite, respectively, and are the azimuth angles of the sun and the satellite, respectively.

[0038] The relationship between the atmospheric Rayleigh scattering phase function P r and the scattering angle can be expressed as: (17) The aerosol scattering phase function P a is determined by using the two-term Henyey-Greenstein function: (18) where [[ID=6:2]] , g1 and g2 are the asymmetry factors, which take different values for different aerosol particles and can be calculated by the following formula:

[0039] For the water surface foam radiance L wc (λ), it can be calculated according to the water surface wind speed W: (19) In the formula, ρ wc is the water surface foam reflectance, which is related to the wind speed W; represents the diffuse transmittance in the solar incident direction; When W ≤ 4 m / s, ρ wc = ; When 4 < W ≤ 7 m / s, ρwc =2.2×10 -5 ρ a C D W 2 -4.0×10 -4 , C D =(0.62+1.56W-1) ×10 -3 ; When W>7 m / s, ρ wc =(4.5×10 -5 ρ a C D -4.0×10 -5 )W 2 , C D =(0.49+0.065W) ×10 -3 ; Among them, ρ a is the atmospheric density (1.2×10 3 g / m 3 ), C D is the tension coefficient; E0 is the solar irradiance at the top of the atmosphere, and θ0 is the solar zenith angle.

[0040] Surface solar flare radiance L g (λ) can be determined based on the geometric position of the sun , sensor observation geometry And the water surface wind speed W is calculated: (20) Where, is the reflectivity of the solar flare, , where ω is the angle at which direct sunlight hits a small plane on the water surface, is the Fresnel reflectivity at the non-polarized incident angle ω, and β is the angle relative to the lowest point of the water surface slope that can produce reflection to the sensor. ω and β can be solved by equations (21) and (22): (twenty one) (twenty two) p w The sun is at a certain position θ0, , the satellite is in a certain observation direction θ v 、 The probability of seeing a solar flare is related to the surface wind speed W, which satisfies the Gram Charlier distribution: (twenty three) Among them, σ c , σ uare the root mean square of the slope of the water surface wave in the horizontal and vertical directions, c 21 and c 03 is the skewness coeffients of the GramCharlier distribution, c 40 、c 22 and c 04 is the peakedness coefficient, c 40 =0.40, c 22 =0.12, c 04 =0.23, ξ and η represent the slope index of the wave; the specific values of each parameter are as follows: , , , , .

[0041] When performing atmospheric correction of OLCI remote sensing images, the atmospheric Rayleigh scattering optical depth is first calculated using the measured atmospheric pressure and visibility data. and aerosol scattering optical depth ; Take the ozone concentration as 300 Dobson and calculate the optical depth of ozone absorption Combined with the observation geometry data (tie_geometries.nc) of the OLCI image, the atmospheric direct transmittance of each channel of the OLCI image is calculated. T (λ) and diffuse transmittance t(λ), atmospheric Rayleigh scattering radiance L r (λ) and aerosol scattering radiance L a (λ); Calculate the water surface foam reflection radiance L from the water surface wind speed data and observation geometry data wc (λ) and the radiance of the solar flare L g (λ), calculate the water-leaving radiance L from the total image radiance L w (λ); Combined with the solar irradiance E0(λ) of the top atmosphere of each channel, the water surface remote sensing reflectivity R of each channel is calculated rs (λ), complete the atmospheric correction of each channel image of OLCI. The atmospheric correction result of OLCI image in this embodiment is compared with the remote sensing reflectance calculated by equivalent calculation of measured spectrum. Figure 2 shown.

[0042] Step 214, extracting water areas from satellite images; The normalized difference water index (NDWI) or improved normalized difference water index (MNDWI) is constructed using the remote sensing reflectance of the green band and near-infrared band or short-wave infrared band after atmospheric correction. The water area and non-water area are distinguished from the atmospheric corrected image through the threshold method.

[0043] Step 214, extracting water areas from satellite images; The normalized difference water index (NDWI) or improved normalized difference water index (MNDWI) is constructed using the remote sensing reflectance of the green band and near-infrared band or short-wave infrared band after atmospheric correction. The water area and non-water area are distinguished from the atmospheric corrected image through the threshold method.

[0044] (twenty four) (25) Where R rs-G 、R rs-NIR and R rs-SWIR Represent the remote sensing reflectance of the green light band, near infrared band and shortwave infrared band respectively.

[0045] Since the OLCI image lacks shortwave infrared data, this example uses the remote sensing reflectance of the OLCI green light Oa6 and near-infrared Oa17 channels to construct the Normalized Difference Water Index (NDWI). A certain threshold (generally, the threshold is greater than 0) is set to extract the water bodies in the Huaihe River basin.

[0046] Step 22, equivalent calculation of OLCI water surface remote sensing reflectance based on the measured spectrum; Combined with the spectral response function f(λ) of each channel of the satellite sensor, the water surface remote sensing reflectance data Rrs(λ) calculated by the measured spectrum of inland water bodies is equivalent to obtain the water surface remote sensing reflectance R of each channel of the satellite sensor rs (ch), the calculation formula is: (26) Where ch represents the central wavelength of each channel of the satellite sensor, λ max ,λ min Indicates the maximum and minimum wavelengths of the sensor's spectral response for each channel.

[0047] The satellite remote sensing reflectance calculated equivalently from the measured spectrum has high accuracy and can be used as the true value to verify the accuracy of the water surface remote sensing reflectance obtained after atmospheric correction of satellite images.

[0048] The 181 remote sensing reflectance spectra sampled on-site from the main lakes and reservoirs in the Huaihe River Basin were combined with the spectral response function of the OLCI 21 channels to equivalently calculate the remote sensing reflectance of each spectrum corresponding to each OLCI channel, laying the foundation for the optical classification of OLCI images and the construction of turbidity remote sensing models.

[0049] Step 23, optical classification of water bodies; Water optical classification uses the K-means method to automatically cluster N water spectral or image data into K clusters, minimizing variance within the same cluster and maximizing variance between different clusters. During the classification process, K samples are randomly selected from the N data sets as initial cluster centers. The distances from the remaining NK data sets to these K cluster centers are then calculated, and the data sets are assigned to the cluster closest to the cluster center, completing a clustering process. The average value of all data within each cluster is then used as the new cluster center for a new clustering process. This cycle continues until the intra-class mean value remains unchanged or the maximum number of iterations is reached, at which point the clustering result is output.

[0050] For the measured water spectral data, the distance function uses the spectral angle distance (SAD) to measure the difference in spectral shape: (27) Where, X (x1, x2,…,x n ) and Y(y1, y2,…,y n ) are two spectral data, x n and y n are the remote sensing reflectances of the nth wavelength in the two spectra respectively.

[0051] For satellite remote sensing images, a method combining Euclidean distance and spectral angle is used to describe the brightness difference of multi-channel remote sensing images. The calculation expression is as follows: (28) Where, G is the distance metric value, x i is the pixel value of the remote sensing image in channel i, y i is the remote sensing reflectance value of the ith channel of the remote sensing image equivalent to the measured spectrum, and n is the number of channels (i.e., the number of bands) of the remote sensing image; is the Euclidean distance between the satellite remote sensing image spectrum and the equivalent simulated value of the measured spectrum; It represents the degree of spectral angle matching between the satellite remote sensing image spectrum and the equivalent simulated value of the measured spectrum.

[0052] For the measured spectral data of water bodies, 181 spectra with wavelengths ranging from 400 to 900 nm were classified using the spectral angle distance as the classification indicator. The K-means clustering algorithm (k-means++) was used, with the initial number of categories set to 10, the maximum number of iterations to 300, and the iteration termination accuracy to 10. -4 After clustering the measured spectra of water bodies, the optical categories of water bodies were finally determined to be three categories (k=3) by examining the measured spectral shapes and turbidity size distribution characteristics in the Huaihe River Basin. The water body spectrum classification results are shown in the figure. Figure 3 shown.

[0053] For the pre-processed OLCI images of the Huaihe River Basin, the image data of channels 2-12 were selected for optical classification. G, which is a combination of Euclidean distance and spectral angle, was used as the classification index. The K-means clustering algorithm (k-means++) was also used to divide the image water body categories into three categories. The classification results of the OLCI images of the three major lakes in the Huaihe River Basin are shown below. Figure 4 shown.

[0054] Step 3: Construction of inland water turbidity remote sensing model; Based on the optical classification of water bodies, measured turbidity and equivalently calculated satellite multi-channel remote sensing reflectance were classified according to different optical types of measured spectra. Correlation analysis was performed between all turbidity values of the same type and the remote sensing reflectance of a single channel or a combination of channels. Channel remote sensing reflectance data with high correlations with turbidity were selected as samples. Turbidity remote sensing models were established using regression analysis and three machine learning methods (neural network, random forest, and XGBoost algorithms), with turbidity as the output parameter. All samples of the same type were divided into training and test sets in a 2:1 ratio. The optimal turbidity remote sensing model was selected based on the test results using three evaluation metrics: correlation coefficient (R), root mean square error (RMSE), and mean absolute error (MAE).

[0055] In this example, the 181 spectra and turbidity values measured in the Huaihe River basin were classified into different optical categories, with 21 samples for Type I and 80 samples for Types II and III. For each category, the remote sensing reflectance of each OLCI channel, calculated from the measured spectra, was equivalently calculated using water color channel data from channels 2-12. Different band combinations, such as single-band data, two-band ratios, and multi-band combinations, were constructed. Correlation analysis was then performed with the measured turbidity data to determine the sensitive bands for turbidity remote sensing. Table 2 lists the bands or band combinations that correlate well with turbidity under different optical categories.

[0056] Table 2 Correlation between different OLCI band combinations and turbidity

[0057] As shown in Table 2, in Type I water bodies, the three band combinations of Oa11-Oa12, Oa7+Oa11, and Oa7 / (Oa2 / Oa11) have strong correlations with turbidity, with correlation coefficients exceeding 0.8. In Type II water bodies, the three band combinations of Oa8 / Oa6, Oa8 / (Oa6+Oa12), and (Oa6+Oa10) / Oa8 have strong correlations with turbidity, with correlation coefficients exceeding 0.8. In Type III water bodies, the four band combinations of Oa8 / Oa2, Oa10+Oa11, (Oa2-Oa8) / Oa6, and Oa8 / (Oa2+Oa7) have good correlations with turbidity, with correlation coefficients exceeding 0.74.

[0058] For the three water body optical types with good correlation with turbidity band combination, with turbidity as the dependent variable and the remote sensing reflectance of each band combination as the independent variable, linear regression statistics and machine learning (random forest and XGBoost algorithms) were used to establish turbidity remote sensing models. When modeling, 14 samples of type I water bodies were randomly selected as training samples and the remaining 7 were used as test samples; 54 samples of type II and III water bodies were randomly selected as training samples and the remaining 26 were used as test samples. The coefficient of determination (R 2 ), root mean square error (RMSE) and mean absolute error (MAE) are used to calibrate the best turbidity remote sensing model for different water body types.

[0059] (i) Construction of turbidity remote sensing model based on statistical regression; For each band combination with good correlation with turbidity selected from Table 2, a linear regression model was built for the training samples, and the accuracy of the model was verified using the test samples. The results are shown in Table 3. The best regression model and verification scatter plot of turbidity remote sensing for different water body types are shown in Table 3. Figure 5 shown.

[0060] Table 3 Results of turbidity remote sensing regression models for different water body types

[0061] (ii) Construction of turbidity remote sensing model based on random forest; The remote sensing reflectance of the 3-4 band combinations with high correlation with turbidity selected from Table 2 was used as the input feature variable, and turbidity was used as the output variable. The turbidity remote sensing model was constructed using the random forest algorithm. Through model parameter tuning, the two main parameters of the turbidity random forest model: the number of decision trees was 100 and the minimum number of leaf nodes was 5. The scatter plots of the training results and testing of the turbidity model for the three water body types are shown in Figure 2. Figure 6 .

[0062] (iii) Construction of turbidity remote sensing model based on XGBoost algorithm; The remote sensing reflectance of the 3-4 band combinations with high correlation with turbidity selected from Table 2 was used as the input feature variable, and turbidity was used as the output variable. The XGBoost algorithm (eXtreme Gradient Boosting) was used to construct the turbidity remote sensing model. During model construction, the model parameters were set to 100 decision trees (n_estimators), 0.1 learning rate (learning_rate), 3 maximum tree depth (max_depth), 1 regularization parameter (reg_lambda), and 0 leaf node weight coefficient (gamma). The scatter plots of the training and testing results of the turbidity model for the three water body types are shown in Figure 2. Figure 7 .

[0063] Using the coefficient of determination (R 2 ), root mean square error (RMSE) and mean absolute error (MAE) three evaluation indicators, by comparing the turbidity remote sensing models established by statistical regression, random forest and XGBoost algorithms for each water body type, the results showed that the turbidity remote sensing model of type I water body was the best when the statistical regression model with Oa7 / (Oa2 / Oa11) as the variable, the random forest model with Oa8 / Oa6, Oa8 / (Oa6+Oa12) and (Oa6+Oa10) / Oa8 as the variables was the best when the type II water body was the best when the XGBoost model with Oa8 / Oa2, Oa10+Oa11, (Oa2-Oa8) / Oa6 and Oa8 / (Oa2+Oa7) as the variables was the best when the type II water body was the best.

[0064] Step 4: remote sensing inversion of turbidity in three typical lakes in the Huaihe River Basin; For the satellite remote sensing images of inland water bodies after preprocessing and optical classification, the satellite images are first divided into regions according to the optical classification results, and the water turbidity is inverted for the pixels in the region according to the corresponding turbidity remote sensing model. Finally, the turbidity results inverted in different regions are spliced and fused to obtain the turbidity remote sensing inversion results of the entire water area.

[0065] For the OLCI remote sensing images containing Hongze Lake, Gaoyou Lake and Luoma Lake after preprocessing and optical classification, the water bodies in the OLCI images are first divided into regions, and the corresponding optimal turbidity remote sensing model is selected for the pixels in the region to invert the water turbidity. Finally, the turbidity results inverted from different water images are spliced and fused to obtain the turbidity remote sensing inversion results of the entire water area. This example inverts the turbidity spatial distribution of three typical lakes in summer 2020 based on OLCI images. Figure 8 .

Claims

1. A satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation, characterized in that: The steps are as follows: S1, obtain water spectral data, water turbidity, satellite remote sensing images and meteorological data of the area to be measured; calculate the water remote sensing reflectance based on the acquired water spectral data; S2, pre-process the satellite remote sensing image to obtain the water area; remote sensing reflectance data of the measured spectrum of the water body Perform equivalent calculations to obtain the remote sensing reflectivity of each channel of the satellite sensor Clustering the measured spectra of water bodies, and combining the measured spectral shapes and turbidity size distribution characteristics of water bodies in the study area to perform optical classification of water bodies; S3, for optical types with different measured spectra, classify them according to the measured turbidity and the equivalent calculated satellite multi-channel remote sensing reflectance, perform correlation analysis on all turbidities of the same type and the remote sensing reflectance of a single channel or a combination of multiple channels, select the channel remote sensing reflectance data with high correlation with turbidity as samples, use turbidity as the output parameter, and construct a turbidity remote sensing model; divide the samples into training and test sets, train and test the turbidity remote sensing model, and screen the best turbidity remote sensing model.

2. The inland water turbidity satellite remote sensing method based on optical classification and spectral simulation according to claim 1 is characterized in that: In step S2, the preprocessing of the satellite remote sensing image includes radiometric correction, geometric correction, atmospheric correction and water area extraction of the image; The specific steps for atmospheric correction of remote sensing images are as follows: S31, calculate the atmospheric Rayleigh scattering optical depth and aerosol scattering optical depth using measured atmospheric pressure and visibility data; calculate the ozone absorption optical depth based on ozone concentration; S32, combining the observation geometry data of remote sensing images, calculates the atmospheric direct transmittance and diffuse transmittance, atmospheric Rayleigh scattering radiance and aerosol scattering radiance of each channel of satellite images; S33, calculate the water surface foam reflection radiance and solar direct flare radiance based on the water surface wind speed data and observation geometry data, and calculate the water-leaving radiance based on the total radiance of the image; S34, combined with the solar irradiance of the top atmosphere of each channel, calculates the water surface remote sensing reflectivity of each channel and completes the atmospheric correction of the satellite images of each channel.

3. The satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation according to claim 2 is characterized in that: The normalized difference water index (NDWI) or modified normalized difference water index (MNDWI) is constructed using the remote sensing reflectance of the green band and near-infrared band or short-wave infrared band after atmospheric correction. The water area and non-water area are distinguished from the atmospheric corrected image by the threshold method. ， ， Where, 、 and Represent the remote sensing reflectance of the green light band, near infrared band and shortwave infrared band respectively.

4. The satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation according to claim 1 is characterized in that: In step S2, the spectral response function of each channel of the satellite sensor is combined , water surface remote sensing reflectance data of the measured spectrum of water bodies Perform equivalent calculations to obtain the water surface remote sensing reflectivity of each channel of the satellite sensor , the calculation formula is: ， Where ch represents the central wavelength of each channel of the satellite sensor, 、 They represent the maximum wavelength and minimum wavelength of the spectral response of each channel of the sensor.

5. The satellite remote sensing method for inland water turbidity based on optical classification and spectral simulation according to claim 1 is characterized in that: In step S3, all samples of the same type are divided into training set and test set in a ratio of 2:1, and the turbidity remote sensing model is established using regression analysis, neural network, random forest and XGBoost algorithms respectively; three evaluation indicators, namely correlation coefficient R, root mean square error RMSE and mean absolute error MAE, are used to screen the best turbidity remote sensing model according to the test results.

Images