Method, system, terminal and medium for obtaining water body attenuation coefficient

Abstract

The application provides a water body diffuse attenuation coefficient acquisition method, system, terminal and medium, wherein the method comprises: according to the water depth distribution and the bottom type distribution of a working area, screening a plurality of sampling points with the same bottom type and different depths as target sampling points; based on L1-level satellite image data to be inverted, extracting the TOA value corresponding to the target sampling points, combining the water depth to construct a sample data pair, using the sample data pair corresponding to the target sampling points to solve a water body diffuse attenuation coefficient model constructed in advance, and acquiring a water body diffuse attenuation coefficient value according to a solving result. The application can quickly and conveniently acquire the water body diffuse attenuation coefficient of an optical shallow water area without a complex atmospheric correction process, thereby effectively supporting long-time-series water quality monitoring of the optical shallow water area.

Description

Technical Field

[0001] This invention relates to the field of water color remote sensing analysis technology, and in particular to a method, system, terminal, and computer storage medium for obtaining the diffuse attenuation coefficient of water bodies. Background Technology

[0002] The diffuse attenuation coefficient of water is an important optical parameter in water color remote sensing analysis. It is mainly used to quantify the attenuation rate of light signals in water, thereby enabling the inversion of water depth or water quality information.

[0003] Currently, research on obtaining the diffuse attenuation coefficient of water bodies in deep water areas is relatively mature. However, due to interference from bottom reflection in shallow water areas, the acquisition methods applicable to deep water areas are not suitable for shallow water areas, resulting in relatively few studies on acquisition methods specifically for shallow water areas. Specifically, for shallow water areas, existing methods for obtaining the diffuse attenuation coefficient typically employ semi-analytical methods based on radiative transfer processes (such as the HOPE method). However, since semi-analytical methods require rigorously atmospherically corrected Rrs parameters as input, they rely on strictly atmospherically corrected water color remote sensing data. The accuracy of Rrs parameters is easily affected by atmospheric correction in nearshore high-turbidity and complex aerosol areas. If poorly corrected water color remote sensing images are used, the calculation error will be large, failing to meet the requirements of subsequent inversion and analysis. Furthermore, semi-analytical methods also have high requirements for the number of bands; for example, when using the HOPE method to invert water depth, at least five bands are required to achieve good inversion results.

[0004] Therefore, while low-cost water color remote sensing data, such as China's Gaofen (GF) and HJ series satellites, can support the needs of low-cost, long-term water monitoring, these data often lack specific bands required for rigorous atmospheric correction. For example, GF series imagery typically only has four bands: blue, green, red, and near-infrared, lacking the shortwave infrared or ultraviolet bands required for rigorous atmospheric correction algorithms, thus making rigorous atmospheric correction impossible. Even if some types (such as GF-5 hyperspectral) have shortwave infrared bands, their data acquisition costs are high and their spatiotemporal resolution is low, making them unsuitable for low-cost, long-term monitoring. Furthermore, the limited band data in these water color remote sensing data makes it difficult to support the band count requirements of semi-analytical methods, rendering them unsuitable for low-cost, high-time-series shallow water monitoring.

[0005] Therefore, how to obtain the diffuse attenuation coefficient of water bodies in shallow water areas based on remote sensing image data without strict atmospheric correction has become a technical problem that needs to be solved in this field. Summary of the Invention

[0006] In view of the shortcomings of the prior art described above, the purpose of this invention is to provide a method, system, terminal and computer storage medium for obtaining the diffuse attenuation coefficient of water bodies, in order to solve the problems that the existing acquisition methods rely on images with strict atmospheric correction, have poor compatibility with domestic high-resolution / environmental satellites that have not been strictly atmospheric calibrated, are difficult to support long-term water quality monitoring, and have low accuracy in obtaining the diffuse attenuation coefficient of optical shallow water areas.

[0007] To achieve the above and other related objectives, the present invention first provides a method for obtaining the diffuse attenuation coefficient of a water body, comprising:

[0008] Based on the water depth and sediment distribution in the working area, several sampling points with the same sediment type but different water depths were selected as target sampling points. Based on the L1-level satellite image data to be inverted, the TOA values ​​corresponding to the target sampling points were extracted, and combined with the water depth, sample data pairs corresponding to each target sampling point were constructed. Using the sample data pairs corresponding to the target sampling points, the pre-constructed water diffuse attenuation coefficient model was solved, and the water diffuse attenuation coefficient value corresponding to the working area was obtained based on the solution results.

[0009] In one optional embodiment, the method for obtaining the target sampling point includes:

[0010] The working area is divided into grids to obtain each grid. Based on the water depth distribution, the water depth distribution interval corresponding to each grid is extracted, and grids with water depth differences greater than a depth threshold are selected based on the water depth distribution interval. Based on the bottom sediment distribution, target grids with high reflectivity bottom sediment area greater than an area threshold are extracted from the selected grids. Several target sampling points are extracted from the target grids. The target sampling points are sampling points with different depths and high reflectivity bottom sediment types within the same grid.

[0011] In an optional embodiment, for a single target grid, the method for extracting a plurality of target sampling points from the target grid includes:

[0012] Select initial sampling points in the target grid and obtain the initial water depth corresponding to the initial sampling points; according to the water depth distribution, extract candidate sampling points in areas with depths different from the initial water depth to form a candidate sampling point set; based on the linear relationship between the water depth index and the water depth, filter the candidate sampling point set to obtain the target sampling point; wherein, the water depth index is the natural logarithm ratio between the reflectance values ​​of the green light band and the blue light band.

[0013] In an optional embodiment, the step of filtering the candidate sampling point set based on the linear relationship between the depth index and the water depth includes:

[0014] A sample pair is formed by combining the depth index and water depth corresponding to each candidate sampling point. Based on the sample pair, a linear fitting relationship between the depth index and water depth is constructed. According to the linear fitting relationship, the fitting residual corresponding to each candidate sampling point is calculated, and candidate sampling points with fitting residual values ​​greater than a preset residual threshold are removed from the candidate sampling point set. The residual threshold is an empirical value; for example, the residual threshold is 4.685.

[0015] In one optional embodiment, the step of solving the pre-constructed water attenuation coefficient model using sample data pairs corresponding to the target sampling points includes:

[0016] The water depth H corresponding to each target sampling point is extracted and set as the model input, and the TOA value corresponding to each target sampling point is set as the model output; the model parameters of the water body diffuse attenuation coefficient model are calculated using the least squares method; wherein, the water body diffuse attenuation coefficient model is:

[0017]

[0018] In the formula, , , These are the model parameters to be solved;

[0019] It is the sum of reflectance produced by atmospheric molecules, aerosols, and the interaction between atmospheric molecules and aerosols;

[0020] For the reflectivity of the substrate, The atmospheric equivalent transmittance coefficient; It is the sum of the absorption and backscattering coefficients;

[0021] This refers to the TOA value corresponding to the optical deep-water zone;

[0022] These are the subsurface solar zenith angle and the subsurface satellite observation angle, respectively.

[0023] A0 and A1 are constant terms.

[0024] In an optional embodiment, obtaining the water attenuation coefficient value corresponding to the working area based on the calculation result includes:

[0025] Based on the model parameters, the subsurface solar zenith angle, and the subsurface satellite observation angle, the water body attenuation coefficient is solved using the correlation between the water body attenuation coefficient and the model parameters, the subsurface solar zenith angle, and the subsurface satellite observation angle. The solution is as follows:

[0026]

[0027] In the formula, The diffuse attenuation coefficient of the water body;

[0028] These are the subsurface solar zenith angle and the subsurface satellite observation angle, respectively, satisfying:

[0029]

[0030] In the formula, , These are the surface solar zenith angle and the surface satellite observation angle, respectively.

[0031] In one optional embodiment, the method for obtaining the TOA value corresponding to the optical deep-water zone includes:

[0032] Identify the adjacent deep-water area of ​​the working area; select several effective deep-water pixels in the adjacent deep-water area based on the water depth distribution and image quality; and obtain the TOA value corresponding to the optical deep-water area by combining the TOA values ​​of each effective deep-water pixel.

[0033] A second aspect of this application provides a system for obtaining the diffuse attenuation coefficient of a water body, comprising:

[0034] The sampling point selection module is used to select several sampling points with the same bottom sediment type but different depths as target sampling points based on the water depth distribution and bottom sediment distribution in the working area; the sample pair construction module is used to extract the TOA value corresponding to the target sampling points based on the L1-level satellite image data to be inverted, and construct the sample data pairs corresponding to each target sampling point in combination with the water depth; the coefficient acquisition module is used to solve the pre-constructed water diffuse attenuation coefficient model using the sample data pairs corresponding to the target sampling points, and obtain the water diffuse attenuation coefficient value corresponding to the working area based on the solution results.

[0035] In a third aspect, this application provides a terminal comprising: a processor and a memory, wherein the memory and the processor are communicatively connected; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to enable the terminal to perform any of the above-described methods for obtaining the water body attenuation coefficient.

[0036] Finally, the present invention also provides a computer storage medium storing a computer program that, when executed by a processor, implements the water body attenuation coefficient acquisition method as described above.

[0037] As described above, the method, system, terminal, and computer storage medium for obtaining the water body diffuse attenuation coefficient provided in this application construct target sampling points with the same substrate type but different depths, and use the target sampling points to perform parameter calculation on the water body diffuse attenuation coefficient model. The calculation results are then used to obtain the water body diffuse attenuation coefficient value corresponding to the working area. This method eliminates the dependence of traditional methods on strictly atmospheric corrected image data. It can not only adapt to Level-1 domestic high-resolution / environmental satellite data, but also achieve high-precision inversion of the water body diffuse attenuation coefficient in shallow water areas, thereby better supporting long-term water quality monitoring. Attached Figure Description

[0038] Figure 1 The diagram shown is a flowchart of an embodiment of the method for obtaining the diffuse attenuation coefficient of water body provided in this application;

[0039] Figure 2 The diagram shows a flowchart of one embodiment of the method for obtaining the target sampling points described in this application.

[0040] Figure 3 The diagram shows a flowchart of one embodiment of the method for obtaining the TOA value corresponding to the optical deep-water region.

[0041] Figure 4 The plot is displayed as a scatter plot based on AC-S measured data.

[0042] Figure 5 The chart shows a scatter plot comparing the values ​​of water diffuse attenuation coefficients retrieved from HJ-2 / CCD and GF-6 / WFV data.

[0043] Figure 6 The diagram shown is a structural schematic of the water body diffuse attenuation coefficient acquisition system provided in this application in one embodiment.

[0044] Figure 7 The diagram shown is a structural schematic of the terminal provided by the present invention in one embodiment. Detailed Implementation

[0045] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0046] To address the limitations of existing methods for obtaining diffuse attenuation coefficients in shallow water areas, which rely on strictly atmospheric-corrected remote sensing imagery and thus have limited compatibility with water color remote sensing data, hindering low-cost, long-term water monitoring, and conversely, existing methods cannot achieve high accuracy in obtaining diffuse attenuation coefficients for low-cost imagery suitable for long-term monitoring, this invention provides a method, system, terminal, and computer storage medium for obtaining diffuse attenuation coefficients. This method uses uncorrected top-of-atmosphere (TOA) reflectance data to obtain diffuse attenuation coefficients for shallow water areas.

[0047] The shallow water area, also known as the optical shallow water area, is the region where underwater reflected signals can penetrate the water and be received by the sensor, such as nearshore shallow seas, island and reef waters, shallow lakes, and relatively clear waters inside reservoirs.

[0048] In some optional embodiments, the shallow water area is a region with a water depth of 0-30 meters that includes a highly reflective substrate; the highly reflective substrate is a substrate region with high bottom reflectivity, that is, in the optical shallow water area, the substrate region in which the bottom reflection contribution accounts for no less than 30% of the total water-leaving radiance received by the sensor; for example, the highly reflective substrate includes calcareous sand, coral sand, light-colored rocks, etc.

[0049] It should be noted that the method described in this application is not suitable for all bodies of water, especially for areas where the water is too turbid, the bottom signal is completely invisible, or the water depth hardly changes.

[0050] The top-of-atmosphere (TOA) reflectance data that has not undergone strict atmospheric correction refers to Level-1 (hereinafter referred to as L1) TOA data, that is, top-of-atmosphere reflectance data that has not yet undergone atmospheric scattering / absorption correction; for example, including but not limited to Gaofen-4 (GF-4) data, Gaofen-6 (GF-6) data, Gaofen-1 (GF-1) data, or Jilin-1 data, etc.

[0051] Please see Figure 1 The diagram shows a flowchart of an embodiment of the method for obtaining the water body attenuation coefficient provided in this application.

[0052] S100: Acquire water depth distribution data, bottom sediment distribution data, and L1-level satellite image data to be retrieved in the working area, and perform spatial registration on each data.

[0053] The water depth distribution data is data with a spatial resolution of no more than 30 m, used to characterize the water depth at different locations within the working area, and it has a unified vertical reference.

[0054] The substrate distribution data is used to characterize the substrate type or bottom reflectance characteristics corresponding to different locations within the working area; in this application, the substrate distribution data includes substrate classification data and / or bottom reflectance characteristic data obtained from multispectral remote sensing image recognition.

[0055] The L1-level satellite imagery data to be inverted is an L1-level satellite remote sensing imagery covering the working area, and has undergone radiometric calibration, that is, the original image values ​​in the image have been converted into top atmospheric reflectance (hereinafter referred to as TOA value) to obtain the TOA value distribution corresponding to the working area.

[0056] Specifically, after acquiring water depth distribution data, bottom sediment distribution data, and L1-level satellite imagery data within the work area, spatial registration is performed on each data to ensure that each data has a consistent spatial coordinate reference.

[0057] In one optional embodiment, the water depth data is obtained from airborne lidar data; the bottom water data is obtained from Sentinel-2 / MSI Level-2 data, specifically image data that has undergone geometric correction, radiometric processing and atmospheric correction, and can clearly and accurately show the differences in bottom sediment in shallow water areas.

[0058] The L1-level satellite imagery data used is HJ-2 data, which has only undergone radiometric calibration and geometric correction processing, and has not yet undergone rigorous atmospheric correction.

[0059] To better understand the characteristics of the satellite data above, the information corresponding to each satellite data is summarized in Table 1.

[0060] Table 1 Satellite Data Processing Status

[0061] sensor MSI CCD Data types Surface reflectance DN Center wavelength (nanometers) 490 485 Wavelength range (nanometers) 458-523 430-520 Spatial resolution (meters) 10 16 Time resolution (days) 5 4

[0062] It should be noted that although sediment distribution data and water depth distribution data may be relatively stable within a certain period and can be reused for multiple satellite images to be processed, when the time interval between the acquisition of each data is too long, there may still be significant changes in sediment distribution and / or water optical conditions, which may affect the accuracy of the water diffuse attenuation coefficient acquisition results.

[0063] Based on this, in some preferred embodiments, during the data collection process, data matching based on the acquisition time is performed on each data point to ensure that the acquisition time of each data point is within a preset time interval. Specifically, water depth distribution data within the working area is acquired, and the acquisition time corresponding to this data is extracted. Based on the acquisition time, satellite imagery data with a time interval less than the acquisition time and a data processing level of L1 are collected. Similarly, sediment distribution data with a time interval less than the acquisition time is collected, thereby reducing the impact of water quality changes, sediment changes, and depth changes on the calculation results. For example, the time interval threshold is 5 years.

[0064] S200: Based on the water depth distribution and bottom sediment distribution, select several sampling points with the same bottom sediment type but different depths within the working area to construct a target sampling point set;

[0065] Specifically, the working area is divided according to a preset grid size to obtain multiple candidate grids. For example, the preset grid size is 500 m × 500 m to adapt to the spatial resolution of satellite data.

[0066] For each candidate grid, based on the water depth distribution data of the working area, the effective water depth values ​​located within that candidate grid are extracted, and the maximum water depth among the effective water depth values ​​is determined. and minimum water depth And obtain the water depth change corresponding to all candidate grids. ,for:

[0067]

[0068] Candidate grids with water depth variations greater than a preset depth difference threshold are selected and identified as the first candidate grids; for example, the preset depth difference threshold is 2 m, so as to better perform subsequent fitting processing.

[0069] Based on the substrate distribution data of the working area, target grids containing substrate areas with high reflectivity greater than an area threshold are extracted from the selected candidate grids; for example, the area threshold is 100 m. 2 .

[0070] For a single target grid, sampling points with significant differences in water depth and high reflectivity substrate are extracted from the same grid based on the water depth distribution corresponding to that grid. These points are then used as target sampling points. The high reflectivity substrate areas include bright sandy substrate areas or other areas where bottom reflection signals can be effectively identified by remote sensing imagery.

[0071] In one specific implementation, 30 to 50 sampling points with the same substrate type (all high reflectivity substrate) and different water depths are selected within a single target grid and designated as target sampling points; the 30 to 50 target sampling points are preferably distributed across multiple water depth ranges from shallow to deep to avoid the sampling points being concentrated in a narrow water depth range.

[0072] The methods for obtaining sampling points with the same / similar substrate types include:

[0073] Visual interpretation of Sentinel-2 / MSI multispectral image data (Level 2) and / or image classification or spectral feature recognition methods are used to determine the substrate type or bottom reflectance characteristics corresponding to each pixel; based on the determination results, target sampling points with the same substrate type or similar bottom reflectance characteristics are selected.

[0074] Furthermore, to more accurately and efficiently screen target sampling points with similar or identical substrate types, and to avoid interference from differences in reflectivity caused by significant differences in substrate type on the inversion of diffuse attenuation coefficients, in some optional embodiments, the method for obtaining the target sampling points is as follows: Figure 2 As shown, it includes the following steps:

[0075] S201, Select the initial sampling point in the target grid and obtain the initial water depth corresponding to the initial sampling point;

[0076] Within the target grid cell, based on the water depth distribution data, a cell with a suitable water depth is pre-selected as the initial sampling point, and its corresponding water depth value is obtained.

[0077] S202, Based on the water depth distribution, candidate sampling points are extracted in areas with depths different from the initial water depth to form a candidate sampling point set;

[0078] Within the same grid cell, based on the water depth distribution characteristics, in the region where the water depth difference from the initial sampling point is greater than the depth difference threshold, a new sampling point is extracted as a candidate pixel, and the corresponding water depth value is obtained.

[0079] Repeat this step to obtain multiple candidate sampling points to form a candidate point set.

[0080] S203, based on the linear relationship between the depth index and the water depth, the candidate sampling point set is screened to obtain the target sampling point set.

[0081] Specifically, the green band reflectance and blue band reflectance of each candidate sampling point in the candidate sampling point set are obtained in the Sentinel-2 / MSI level multispectral image (L2 level);

[0082] Based on the reflectance of the green light band and the reflectance of the blue light band, the water depth index (WDI) corresponding to each candidate sampling point is calculated respectively; the water depth index (WDI) is used to characterize the ratio between the reflectance of the candidate sampling point in the green light band and the reflectance in the blue light band.

[0083] In one specific embodiment, the water depth index (WDI) is the ratio of the natural logarithm of the reflectance values ​​of the green light band to the blue light band, i.e.:

[0084]

[0085] In the formula, Represents the water depth index of the i-th candidate sampling point; This represents the reflectance value of the green light band corresponding to the i-th candidate sampling point; This represents the reflectance value of the blue light band corresponding to the i-th candidate sampling point.

[0086] For shallow water areas with similar substrate types and uniform optical properties, the WDI value exhibits a stable linear relationship with water depth:

[0087] H=f(WDI)

[0088] Based on this, in this embodiment, the water depth index corresponding to the candidate sampling point is... With water depth Form a sample pair to obtain the sample pair corresponding to each candidate sampling point; based on the sample pair, construct a linear fitting relationship between the depth index and the depth, as follows:

[0089]

[0090] In the formula, The slope of the linear fitting relationship is denoted as . The intercept represents the linear fit.

[0091] Based on the linear fitting relationship, the fitting residuals corresponding to each candidate sampling point are calculated, so as to identify abnormal sampling points based on the fitting residuals.

[0092] More specifically, after obtaining the fitting residuals corresponding to each candidate sampling point, candidate sampling points whose fitting residual values ​​are greater than a preset residual threshold are removed from the candidate sampling point set as abnormal sampling points in order to construct a new candidate sampling point set.

[0093] Based on the new candidate sampling point set after removing outlier sampling points, the above linear fitting process is re-executed to establish a new linear fitting relationship; and the outlier sampling point removal process is re-executed until the fitting residual of each candidate sampling point satisfies the preset residual condition; then, the multiple candidate sampling points that satisfy the linear fitting relationship are determined as the final target sampling points, and the target sampling point set is obtained.

[0094] Among them, the target sampling points in the target sampling point set have the same or similar bottom reflection characteristics, similar water optical properties, and different water depths.

[0095] It should be noted that after removing abnormal sampling points, it is checked whether the number of remaining candidate sampling points is less than a preset threshold. If so, the process returns to step S201, new candidate sampling points are extracted again within the target grid, and steps S201 to S203 are repeated until the requirement is met. The threshold value is adapted to the size of a single grid; for example, when the area of ​​a single grid is 100 m²... 2 At that time, the quantity threshold is 30.

[0096] By utilizing the spectral characteristics of candidate sampling points, it is possible to quickly verify whether the candidate sampling points have the same substrate type, ensuring that the sampling point set constructed within the same grid cell has the same substrate reflectance characteristics. This effectively eliminates the interference of substrate differences on the water diffuse attenuation coefficient inversion model and effectively improves the inversion accuracy.

[0097] S300: Based on the L1-level satellite imagery data to be inverted, extract the TOA values ​​corresponding to the target sampling points; combine the TOA values ​​corresponding to the target sampling points with the water depth to construct sample data pairs corresponding to each target sampling point;

[0098] Specifically, based on the coordinate information of the target sampling point, the TOA value of the pixel corresponding to the coordinate is extracted from the L1 level satellite image data to be inverted and set as the TOA value corresponding to the target sampling point.

[0099] For each target sampling point, the water depth value obtained by S200 is combined with its corresponding TOA value to construct a sample data pair corresponding to that target sampling point.

[0100] S400, using the sample data pairs corresponding to the target sampling points, the pre-constructed water body diffuse attenuation coefficient model is solved, and the water body diffuse attenuation coefficient value corresponding to the working area is obtained based on the solution results.

[0101] The water body attenuation coefficient model is a function characterizing the relationship between TOA value and water depth, including:

[0102]

[0103]

[0104] In the formula, , , These are the three parameters that need to be fitted, and they change with variations in water depth and TOA value. , , Remain unchanged; among which, It has an initial value and a corresponding floating range is set during the fitting process to constrain the fitting results.

[0105] It is the sum of reflectance produced by atmospheric molecules, atmospheric aerosols, and the interactions between atmospheric molecules and aerosols, that is:

[0106]

[0107] R r R a and R ra These represent the reflectance produced by atmospheric molecules, atmospheric aerosols, and the interaction between atmospheric molecules and aerosols, respectively.

[0108] For the reflectivity of the substrate, The atmospheric equivalent transmittance coefficient;

[0109] It is the sum of the absorption and backscattering coefficients;

[0110] This refers to the TOA value corresponding to the optical deep-water zone; in this application, this... Set as model parameters The initial values ​​are determined based on the sample data corresponding to each target sampling point, and the model parameters are adjusted using the nonlinear least squares method. , , Perform a joint solution to achieve the desired result. Constraints corresponding to the solution results.

[0111] These are the subsurface solar zenith angle and the subsurface satellite observation angle, i.e., the solar zenith angle below the water surface and the satellite observation angle below the water surface;

[0112] A0 and A1 are two constant terms, which are 1 and 1 / π respectively when the bottom of the water body is a Lambertian surface.

[0113] It should be noted that, in this application, since the area of ​​a single grid is small, the atmospheric conditions, substrate reflectivity and water optical properties corresponding to different sampling points within the same grid can be considered the same.

[0114] Specifically, the water depth H corresponding to the target sampling point is used as the model input, and the TOA value corresponding to the same target sampling point is used as the model output, as follows:

[0115]

[0116] Then, based on the extracted sample data, the above equation is evaluated using the least squares method to obtain... , , The corresponding value;

[0117] Based on the correlation between the water body attenuation coefficient and c2, the water body attenuation coefficient is calculated to obtain the final water body attenuation coefficient, which is:

[0118]

[0119] In the formula, Let be the subsurface solar zenith angle and the subsurface satellite observation angle, respectively. These angles, along with their corresponding surface solar zenith angle and surface satellite observation angle, satisfy the following condition:

[0120]

[0121] in, , The surface solar zenith angle and surface satellite observation angle are metadata included with the TOA data.

[0122] It should be noted that, due to the exponential nonlinear fitting of the water depth inversion model, the fitting parameters include... , , The initial value and upper and lower limits directly affect the fitting convergence and inversion accuracy.

[0123] In one specific implementation, the TOA value corresponding to the optical deep-water area adjacent to the working area in the same image is set as... The initial values; based on this, the fitting parameters and The methods for obtaining the initial value include: First, Set its corresponding initial value to obtain ln( )and The linear relationship between them; Then, by performing a linear fit between the water depth at the target sampling point and the TOA value corresponding to the same target sampling point, the solution is obtained. and The initial value.

[0124] And, the fitting parameters and The constraint range is determined as follows: The lower limit corresponds to the pure water diffuse decay coefficient. The minimum value, The upper limit corresponds to the substrate reflectivity. Atmospheric equivalent transmittance The value that simultaneously reaches the theoretical maximum value; Furthermore, to ensure the physical rationality of the fitting process, the TOA value corresponding to the optical deep-water region (i.e., The value is set to be less than the minimum TOA value of the model output to ensure that the fitting process does not exceed the actual physical range of remote sensing observation.

[0125] Because in the water body diffuse attenuation coefficient model, The initial value is the TOA value corresponding to the optical deep-water region adjacent to the working area, i.e. To improve the accuracy of model solutions, in some optional embodiments, the method for obtaining the TOA value corresponding to the optical deep-water region is as follows: Figure 3 As shown, it includes:

[0126] S401, Determine the adjacent deep water area of ​​the work area;

[0127] The adjacent deep-water area refers to an optical deep-water area located within a preset space around the work area, and having the same or similar water optical properties and atmospheric conditions as the work area.

[0128] It should be noted that the adjacent deep-water area and the working area are located in the same satellite image to be inverted, and have atmospheric top reflectivity data in the same band.

[0129] S402, based on the water depth distribution and image quality, select several effective deep-water pixels in the adjacent deep-water area;

[0130] The effective deep-water pixel is defined as a water depth greater than a preset deep-water threshold, so that the influence of underwater reflection on the atmospheric top reflectivity of the candidate deep-water pixel can be ignored; for example, the deep-water threshold is 30 meters.

[0131] Specifically, in the adjacent deep water area, several pixel points with a water depth greater than a preset deep water threshold are extracted as candidate deep water pixels; the candidate deep water pixels are subjected to quality screening to remove pixels affected by land, islands, reefs, cloud shadows, waves, suspended sediment, etc., in order to obtain an effective deep water pixel set.

[0132] S403, by combining the TOA values ​​of each effective deep-water pixel, obtains the corresponding TOA value for the optical deep-water region.

[0133] Specifically, the TOA values ​​of each effective deep-water pixel in the effective deep-water pixel set are extracted in the target band, the extracted TOA values ​​are statistically processed, and the statistical results are determined as the TOA values ​​of the adjacent deep-water area.

[0134] For example, the median or average TOA value corresponding to each effective deepwater pixel is calculated to reduce the impact of random noise and local anomalous water bodies on the corresponding TOA value in the optical deepwater region.

[0135] For satellite imagery such as HJ-2 / CCD and GF-6 / WFV, due to their low radiometric resolution and limited response to low-emissivity signals, the measured values ​​of atmospheric top reflectivity may be too high or fluctuate significantly in optical deep-water regions. This measurement error reduces the accuracy of TOA values ​​in deep-water regions, thus affecting the accuracy of the water diffuse attenuation coefficient inversion. To improve the measurement accuracy of the satellite imagery to be inverted in optical deep-water regions, and further improve the accuracy of water diffuse attenuation coefficient inversion, in some optional embodiments, the TOA data acquisition method includes:

[0136] Obtain radiometric calibration data that matches the satellite image data to be inverted, and use this radiometric calibration data to perform relative radiometric calibration on the satellite image data to be inverted, so as to perform subsequent steps using the relatively radiometrically calibrated satellite image data.

[0137] Among them, the radiometric calibration data is the image data that has completed the radiometric calibration at the sensor end, such as the L1B data of Sentinel-3 / OLCI. The OLCI L1B (Full Resolution / Reduced Resolution) is the TOA data obtained from the instrument digital values ​​after processing such as geometric positioning, spectral calibration, nonlinear correction and absolute radiometric calibration.

[0138] In one specific embodiment, HJ-2 and GF-6 data located within the working area are acquired, and radiometric calibration of the images is performed on them, that is, the DN values ​​in the original satellite images are converted into TOA (Top of Atmosphere Reflectivity).

[0139] Acquire Sentinel-3 / OLCI L1B data that matches the satellite image to be inverted in terms of spatial range, imaging time, and spectral band; wherein, the Sentinel-3 / OLCI L1B data is reference radiation data that is comparable to the reflectivity of the top of the atmosphere to be corrected;

[0140] Spatial registration, spatial resolution matching, and band matching are performed on the satellite image to be inverted and the Sentinel-3 / OLCI L1B data, and the relative radiometric conversion relationship between the two is established based on the matching results.

[0141] Based on the relative radiative conversion relationship, the atmospheric top reflectivity to be corrected is subjected to relative radiative correction to obtain radiatively corrected atmospheric top reflectivity data; wherein, the corrected atmospheric top reflectivity data is used to extract the TOA value of the deep water area and the TOA value of each target sampling point, so as to solve the water body diffuse attenuation coefficient model.

[0142] The relative radiation calibration process described herein employs existing relative radiation calibration methods, which will not be elaborated upon here.

[0143] Since Sentinel-3 / OLCI L1B image data has high radiometric quality, using it as a reference for relative radiometric correction of HJ-2 / CCD or GF-6 / WFV images can reduce the high low reflectivity and abnormal fluctuations in the optical deep water region of the satellite image to be inverted. This improves the accuracy of C3 in the fitting parameter, making the fitting curve closer to the actual water decay process, thereby improving the accuracy of the water diffuse decay coefficient obtained by fitting.

[0144] In this application, the water diffuse attenuation coefficient model is constructed based on the optical shallow water radiative transfer mechanism. To further demonstrate the objectivity of the model construction and the objectivity of the water diffuse attenuation coefficient derived from it, the construction process is explained in detail below:

[0145] The actual reflectance obtained by the satellite is the TOA (Total Optical Angle) at the top of the atmosphere, which is:

[0146] (3-1)

[0147] make The reflectivity TOA at the top of the atmosphere is:

[0148] (3-2)

[0149] Among them, R r R a and R ra These represent the reflectance produced by atmospheric molecules, atmospheric aerosols, and the interaction between atmospheric molecules and aerosols, respectively. Atmospheric path reflectance is the sum of reflectance produced by atmospheric molecules, atmospheric aerosols, and the interaction between atmospheric molecules and aerosols.

[0150] The remote sensing reflectance above the water surface at a water depth of H;

[0151] is the atmospheric equivalent transmission coefficient, used to characterize the atmospheric transmission effect when a signal reflected from the seabed propagates to the top of the atmosphere.

[0152] For shallow water areas, since the water properties are vertically homogeneous and the bottom substrate is of high reflectivity, it conforms to the Lambertian reflectance characteristics. Therefore, the subsurface remote sensing reflectance model for shallow water areas is as follows:

[0153] (3-3)

[0154] in, The remote sensing reflectance of the subsurface layer (below the water surface) at a water depth of H;

[0155] Represents the remote sensing reflectance of the subsurface layer (below the water surface) in the optical deep water region;

[0156] A0 and A1 are two constant terms, which are 1 and 1 / π respectively when the bottom of the water body is a Lambertian surface;

[0157] ρ is the water reflectance;

[0158] H represents water depth;

[0159] This is the equivalent water body attenuation parameter;

[0160] The path factor is formed by the downward propagation of light and the upward propagation of light reflected from the bottom of the water.

[0161] After simplifying the above model, it becomes:

[0162] (3-4)

[0163] For path factors Although a complete shallow water radiative transfer model should include an optical path extension factor in the path factor, to reduce dependence on prior water optical parameters such as absorption coefficient and backscattering coefficient, and under the condition that the optical properties of the water within the target grid are approximately uniform, the optical path extension factor can be simplified by reducing the shallow water radiative transfer relationship to a single exponential model, as follows:

[0164] (3-5)

[0165] in, The subsurface solar zenith angle; This is the observation angle for subsurface satellites.

[0166] It should be noted that the subsurface angles (solar zenith angle and satellite observation angle) are not directly equivalent to the surface angles (solar zenith angle and satellite observation angle) in the satellite image metadata. They satisfy the refraction relationship as follows:

[0167]

[0168]

[0169] in, These are the surface solar zenith angle and the surface satellite observation angle, respectively. The refractive index of water ranges from 1.33 to 1.4; for example, when... If we take 1.33, then:

[0170]

[0171]

[0172] Under low reflectivity conditions, r in shallow water areas rs Typically less than 0.1, based on r rs and R rs The transformation relationship between them for R rs After simplification, we get R. rs Approximately 0.5 times r rs ,Right now:

[0173] (3-7)

[0174] Substituting Equation 3-7 into the aforementioned shallow water model (Equation 3-3), we obtain:

[0175] (3-8)

[0176] Will Substituting into the atmospheric top reflectivity model (3-2), we get...

[0177] (3-9)

[0178] In deep water areas, as the water depth becomes sufficiently large, the contribution of bottom reflection approaches zero, i.e. As the value of approaches zero, the reflectivity of the top of the atmosphere approaches:

[0179] (3-10)

[0180] The atmospheric top reflectivity model can then be further simplified to:

[0181] (3-11)

[0182] To facilitate fitting using water depth and satellite TOA data, the physical parameters in the model are merged into corresponding model parameters, i.e.:

[0183] (3-12)

[0184] Therefore, the water attenuation coefficient calculation model used in the application embodiment is obtained as follows:

[0185] (3-13)

[0186] (3-14)

[0187] Due to the water body diffuse decay coefficient K d The total water body attenuation κ satisfies:

[0188] (3-15)

[0189] Combining with equation 3-14 above: We can obtain:

[0190] (3-16)

[0191] In the formula, ρ, ∑R, and κ are the same in the same wavelength band. TOA dp It is the atmospheric top reflectivity in the optical deep-water region; within a small area, atmospheric conditions (t and () can be considered the same.

[0192] To verify the effectiveness and accuracy of the water body attenuation coefficient acquisition method provided in this application, this embodiment uses in-situ measurement data to verify the water body attenuation coefficient inversion results obtained based on the above method.

[0193] The in-situ measurement equipment used is the AC-S water absorption attenuation meter, which is a water absorption attenuation meter developed by Sea-Bird Scientific (formerly WET Labs). After pure water calibration, the AC-S water absorption attenuation meter can be used to measure the water absorption coefficient and attenuation coefficient at different depths and multiple wavebands, and obtain the corresponding measured water diffuse attenuation coefficient based on this.

[0194] Specifically, six in-situ measurement stations were set up in the Qilianyu Islands study area of ​​the South China Sea, and AC-S in-situ measurement data were collected on November 30, 2022. The water quality of the Qilianyu Islands is relatively clear, with the measured water diffuse attenuation coefficients at each station ranging from 0.03 to 0.07 m⁻¹, and the water depths in the areas where each station is located ranging from 10 to 30 m. A GF-6 / WFV Level 1 image from November 29, 2022, which is close to the in-situ measurement time, was selected as the image to be inverted. Exponential fitting was performed using the water depth and atmospheric top reflectivity at each target sampling point to obtain the corresponding water diffuse attenuation coefficients. The inversion results were compared with the measured results from the six stations. Figure 4 As shown, the fitting results of the water body attenuation coefficient near each station are in good agreement with the measured results, as evidenced by the fitting determination coefficient R. 2 The value is 0.94, and the root mean square error (RMSE) is 0.0021 m. -1 The mean absolute percentage difference (MAPD) ​​was 2.89%. This result indicates that in optically shallow water areas with a water depth of no more than 30 m, where the bottom reflectivity can be effectively identified by remote sensing images and the water depth difference meets the fitting requirements, this method can obtain a water diffuse attenuation coefficient that matches the in-situ measurement results.

[0195] To further verify the applicability of this method to different satellite imagery, this embodiment also utilizes HJ-2 / CCD and GF-6 / WFV imagery for inter-satellite consistency verification. The verification area includes the Florida Keys and the Qilianyu Islands, and the verification period is from January 1, 2021 to December 31, 2024. The temporal matching window for the two types of satellite imagery was set to 5 days, resulting in 42 sets of temporally and spatially matched water attenuation coefficient inversion results. The 42 sets of matching results were compared, as follows... Figure 5 As shown, the coefficient of determination R of the inversion results of HJ-2 / CCD and GF-6 / WFV is... 2 The mean absolute percentage difference (MAPD) ​​was 0.89, and the overall distribution of the results across all groups was near the 1:1 reference line. These results indicate that, under the same region and similar imaging time conditions, the inversion results obtained from HJ-2 / CCD and GF-6 / WFV imagery show good consistency. This method is applicable to the inversion of water diffuse attenuation coefficients in different regions and from different satellite imagery.

[0196] It should be noted that the effective application conditions corresponding to the above verification results include: the bottom reflectivity of the target area is high, preferably high reflectivity substrates such as bright sand; the water quality and atmospheric conditions are the same in the same target grid, and the substrate reflectivity characteristics of the selected target sampling points are similar; there is sufficient water depth difference within the target grid, preferably a water depth variation greater than 2m; and the image does not have severe clouds, cloud shadows, solar flares, or other anomalies that affect the extraction of water reflectivity.

[0197] To address the technical problems existing in the prior art, this invention also provides a system for obtaining the water body diffuse attenuation coefficient to implement the above method, such as... Figure 6 As shown, the water body attenuation coefficient acquisition system 800 includes: a sampling point screening module 810, a sample pair construction module 820, and a coefficient acquisition module 830;

[0198] The sampling point screening module 810 is used to select several sampling points with the same bottom sediment type but different water depths as target sampling points based on the water depth distribution and bottom sediment distribution in the working area. The sample pair construction module 820 is used to extract the TOA value corresponding to the target sampling point based on the L1 level satellite image data to be inverted, and combine it with the water depth to construct sample data pairs corresponding to each target sampling point.

[0199] The coefficient acquisition module 830 is used to solve the pre-constructed water body diffuse attenuation coefficient model using the sample data pairs corresponding to the target sampling point, and obtain the water body diffuse attenuation coefficient value corresponding to the working area based on the solution result.

[0200] Specifically, the execution steps performed by each module to achieve its respective function are the same as those in the above embodiments, and will not be repeated here.

[0201] To address the aforementioned technical problems in the prior art, the present invention also provides a terminal for obtaining the diffuse attenuation coefficient of a water body.

[0202] Please see Figure 7 The diagram shows the structure of the terminal; as shown below. Figure 7 As shown, the terminal 40 includes a processor 41 and a memory 42 connected to each other; the memory is used to store computer programs, and the processor is used to execute the computer programs stored in the memory, so that the terminal can implement all the steps of the water body attenuation coefficient acquisition method as described above when it is executed.

[0203] The specific steps of the method for obtaining the water body attenuation coefficient and the beneficial effects obtained by applying the terminal provided in this embodiment are the same as those in the above embodiments, and will not be repeated here.

[0204] It should be noted that the memory may include random access memory (RAM) and may also include non-volatile memory, such as at least one disk storage device. Similarly, the processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0205] This invention also provides a computer-readable storage medium storing a computer program that, when called by a processor, implements the method for obtaining the water body attenuation coefficient as described in the above embodiments.

[0206] A computer-readable storage medium can be a tangible device capable of holding and storing instructions used by an instruction execution device. Computer-readable storage media can be, for example, (but not limited to) electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, and mechanical encoding devices.

[0207] The computer-readable program described herein can be downloaded from a computer-readable storage medium to various computing / processing devices, or downloaded via a network, such as the Internet, a local area network, a wide area network, and / or a wireless network, to an external computer or external storage device. A network adapter card or network interface in each computing / processing device receives computer-readable program instructions from the network and forwards these instructions to the computer-readable storage medium in the respective computing / processing device.

[0208] The above embodiments are merely illustrative of the principles and effects of the present invention and are not intended to limit the invention. Any person skilled in the art can modify or alter the above embodiments without departing from the spirit and scope of the present invention. Therefore, all equivalent modifications or alterations made by those skilled in the art without departing from the spirit and technical concept disclosed in the present invention should still be covered by the claims of the present invention.

Claims

1. A method for obtaining the diffuse attenuation coefficient of a water body, characterized in that, include: Based on the water depth and sediment distribution in the work area, several sampling points with the same sediment type but different water depths were selected as target sampling points. Based on the L1-level satellite imagery data to be inverted, the TOA values ​​corresponding to the target sampling points are extracted, and combined with the water depth, sample data pairs corresponding to each target sampling point are constructed. Using the sample data pairs corresponding to the target sampling points, the pre-constructed water body diffuse attenuation coefficient model is solved, and the water body diffuse attenuation coefficient value corresponding to the working area is obtained based on the solution results. The step of using the sample data pairs corresponding to the target sampling points to solve the pre-constructed water body diffuse attenuation coefficient model includes: extracting the water depth H corresponding to each target sampling point as the model input, setting the TOA value corresponding to each target sampling point as the model output; and using the least squares method to solve the model parameters of the water body diffuse attenuation coefficient model. The water body diffuse attenuation coefficient model is as follows: ； ； In the formula, A model for the diffuse attenuation coefficient of water bodies; , , These are the model parameters to be solved; It is the sum of reflectance produced by atmospheric molecules, aerosols, and the interaction between atmospheric molecules and aerosols; For the reflectivity of the substrate, The atmospheric equivalent transmittance coefficient; It is the sum of the absorption and backscattering coefficients; This refers to the TOA value corresponding to the optical deep-water zone; These are the subsurface solar zenith angle and the subsurface satellite observation angle, respectively. A0 and A1 are constant terms; H is the water depth.

2. The method for obtaining the water body diffuse attenuation coefficient according to claim 1, characterized in that, The methods for obtaining the target sampling points include: The work area is divided into grids to obtain the resulting grids; Based on the water depth distribution, the water depth distribution interval corresponding to each grid is extracted, and grids with water depth differences greater than a depth threshold are selected based on the water depth distribution interval; and based on the bottom sediment distribution, target grids with high reflectivity bottom sediment area greater than an area threshold are extracted from the selected grids. Several target sampling points are extracted from the target grid. The target sampling points are sampling points within the same grid with different depths and all of them being high-reflectivity substrates.

3. The method for obtaining the water body diffuse attenuation coefficient according to claim 2, characterized in that, The method for extracting several target sampling points from the target grid includes: Select initial sampling points in the target grid and obtain the initial water depth corresponding to the initial sampling points; Based on the water depth distribution, candidate sampling points are extracted in regions with depths different from the initial water depth to form a candidate sampling point set; Based on the linear relationship between the depth index and the water depth, the candidate sampling point set is filtered to obtain the target sampling point; The water depth index is calculated using the natural logarithm ratio between the reflectance values ​​of the green and blue light bands.

4. The method for obtaining the water body diffuse attenuation coefficient according to claim 3, characterized in that, The process of filtering the candidate sampling point set based on the linear relationship between the depth index and water depth includes: The water depth index corresponding to each candidate sampling point is paired with the water depth to obtain the sample pair corresponding to each candidate sampling point. Based on the sample pairs, a linear fitting relationship between the depth index and the water depth is constructed. Based on the linear fitting relationship, the fitting residual corresponding to each candidate sampling point is calculated, and candidate sampling points with fitting residual values ​​greater than a preset residual threshold are removed from the candidate sampling point set.

5. The method for obtaining the water body diffuse attenuation coefficient according to claim 1, characterized in that, The step of obtaining the water body attenuation coefficient value corresponding to the working area based on the solution results includes: based on model parameters The water attenuation coefficient is calculated using the subsurface solar zenith angle and the subsurface satellite observation angle, as follows: ； In the formula, The diffuse attenuation coefficient of the water body; These are the subsurface solar zenith angle and the subsurface satellite observation angle, respectively, satisfying: ； In the formula, , These are the surface solar zenith angle and the surface satellite observation angle, respectively.

6. The method for obtaining the water body diffuse attenuation coefficient according to claim 1, characterized in that, The method for obtaining the TOA value corresponding to the optical deep-water zone includes: Identify the adjacent deep-water areas within the work area; Based on the water depth distribution and image quality, several effective deep-water pixels were selected in the adjacent deep-water area. By combining the TOA values ​​corresponding to each effective deep-water pixel, the TOA value corresponding to the optical deep-water region is obtained.

7. A system for obtaining the diffuse attenuation coefficient of a water body, characterized in that, include: The sampling point screening module is used to select several sampling points with the same bottom sediment type but different water depths as target sampling points based on the water depth distribution and bottom sediment distribution in the working area. The sample pair construction module is used to extract the TOA value corresponding to the target sampling point based on the L1 level satellite image data to be inverted, and combine it with the water depth to construct the sample data pair corresponding to each target sampling point. The coefficient acquisition module is used to solve the pre-constructed water body diffuse attenuation coefficient model using the sample data pairs corresponding to the target sampling point, and obtain the water body diffuse attenuation coefficient value corresponding to the working area based on the solution result. The step of using the sample data pairs corresponding to the target sampling points to solve the pre-constructed water body diffuse attenuation coefficient model includes: extracting the water depth H corresponding to each target sampling point as the model input, setting the TOA value corresponding to each target sampling point as the model output; and using the least squares method to solve the model parameters of the water body diffuse attenuation coefficient model. The water body diffuse attenuation coefficient model is as follows: ； ； In the formula, A model for the diffuse attenuation coefficient of water bodies; , , These are the model parameters to be solved; It is the sum of reflectance produced by atmospheric molecules, aerosols, and the interaction between atmospheric molecules and aerosols; For the reflectivity of the substrate, The atmospheric equivalent transmittance coefficient; It is the sum of the absorption and backscattering coefficients; This refers to the TOA value corresponding to the optical deep-water zone; These are the subsurface solar zenith angle and the subsurface satellite observation angle, respectively. A0 and A1 are constant terms; H is the water depth.

8. A terminal, characterized in that, include: A processor and a memory, wherein the memory and the processor are communicatively connected; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to cause the terminal to perform the method for obtaining the water body attenuation coefficient as described in any one of claims 1 to 6.

9. A computer storage medium having a computer program stored thereon, characterized in that, When executed by the processor, the program implements the method for obtaining the water body attenuation coefficient as described in any one of claims 1 to 6.

Images