A method and system for atmospheric correction of remote sensing images for nearshore marine scenes

Abstract

This application provides a method and system for atmospheric correction of remote sensing images for nearshore marine scenes, relating to the field of atmospheric correction technology for marine remote sensing images. After acquiring remote sensing data, the method extracts a target mask and performs solar flare correction using a secondary simulation model of satellite signals in the solar spectrum. Then, it calculates multi-band luminosity indices and retrieves the optimal aerosol optical thickness (AOT) value using a radiative transfer model. Finally, it calculates atmospheric correction coefficients for each band based on the optimal AOT value to generate the correction result data. The method can perform robust flare correction based on a SWIR proxy, utilizing the physical characteristics of water reflection and path radiation under deep-water conditions to construct a scene-adaptive flare proxy, achieving flare subtraction. It also uses multi-band joint luminosity index inversion, leveraging the spectral characteristics of deep-water body reflection in the visible spectrum to smooth single-band random noise and improve the accuracy of AOT estimation results.

Description

Technical Field

[0001] This application relates to the field of atmospheric correction technology for marine remote sensing images, and in particular to a method and system for atmospheric correction of remote sensing images for nearshore marine scenes. Background Technology

[0002] In the field of marine remote sensing, apparent reflectance (TOA) refers to the reflectance value obtained by normalizing the radiation signal received by a satellite sensor at the top of the atmosphere to solar irradiance. Since the radiation signal received by satellite sensors is the result of multiple effects, including reflection from the ocean surface or water body, atmospheric scattering and absorption, and specular reflection from the water surface, its data quality directly determines the reliability of marine environmental monitoring, marine resource surveys, and global marine carbon cycle research.

[0003] To obtain accurate ocean surface reflectance, atmospheric correction is required for ocean remote sensing data. Atmospheric correction refers to the process of removing atmospheric scattering, absorption, and solar flare effects from the apparent reflectance received from satellite sensors to deduce the true reflectance of the ocean surface or water body. Accurate estimation of aerosol optical thickness (AOT) is a prerequisite for atmospheric correction. Furthermore, specular reflection from the water surface significantly increases apparent reflectance in the visible band. Without effective flare correction, this not only leads to inflated reflectance and color cast but also disrupts the stability of dark pixels in deep water regions, causing systematic biases in AOT estimation.

[0004] The dark pixel method can be used to estimate aerosol optical thickness (AOT). This involves finding the pixel with the lowest reflectance in a single band within an image and directly equating its observed value to atmospheric path radiation (AOT), thus inferring AOT based on the assumption of a single darkest point. However, the dark pixel method relies on the darkest point in a single blue band to infer AOT. Taking the B2 band of the OLI sensor as an example, its estimation error reaches 30% to 50% in near-shore scenes, and this error widens to over 60% under the influence of thin clouds. The dark pixel method is susceptible to interference from noise and sub-pixel bubbles. The darkest point is often an atypical outlier, and the selection of fixed low quantiles is arbitrary, leading to inaccurate and highly volatile AOT estimations and a lack of statistical robustness.

[0005] Regarding flare correction, forward models such as Cox-Munk rely on external parameters like wind speed and wave state, requiring wind speed data with an accuracy better than 2 m / s. When nearshore data is scarce or inaccurate, the correction residual can reach 0.02 to 0.05 reflectance units, and a single image takes over 30 minutes to process, making it difficult to adapt to the needs of large-scale, high-frequency observations. Furthermore, the band-by-band empirical regression method is sensitive to outliers and prone to cross-band color shifts, with color shift indices typically exceeding 0.1, failing to meet the color fidelity requirements for water quality inversion. Therefore, atmospheric correction of nearshore scene remote sensing images faces problems such as inaccurate AOT estimation, poor robustness of flare correction, and insufficient scene adaptability, making it difficult to support high-precision application requirements. Summary of the Invention

[0006] In view of this, embodiments of this application provide a method and system for atmospheric correction of remote sensing images for nearshore marine scenes, in order to solve the problem of inaccurate AOT estimation results.

[0007] According to a first aspect of this application, an atmospheric correction method for remote sensing images of nearshore marine scenes is provided, the method comprising:

[0008] Acquire remote sensing data, which includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light band, near-infrared band, and short-wave infrared band. The auxiliary parameters include solar geometric parameters and sea surface state parameters.

[0009] Target masks are extracted from the remote sensing data, including water masks and deep-water masks;

[0010] Solar flare correction is performed on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum. The secondary simulation model is a model relating the visible light band to the flare intensity proxy based on the flare coefficient. The flare intensity proxy uses a first shortwave infrared band as a specular reflection proxy, constructed in conjunction with reference anchor points. The reference anchor points are obtained by statistically analyzing a preset percentile of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained using robust regression fitting with the median slope method.

[0011] The multi-band luminosity index of the deep-water mask after solar flare correction is calculated. The multi-band luminosity index is calculated based on the pixel value of each pixel in the blue visible light, green visible light and red visible light bands of the deep-water mask.

[0012] The optimal aerosol optical thickness value is retrieved by inverting the multi-band luminosity index using a radiative transfer model. The radiative transfer model is a physical parameter calculation model based on the auxiliary parameters. The physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance.

[0013] The atmospheric correction coefficients for each band are calculated based on the optimal aerosol optical thickness value, and the correction result data is generated based on the atmospheric correction coefficients. The correction result data includes the corrected image, the difference map, and the color composite image.

[0014] By employing the above technical solutions, this application provides an atmospheric correction method for remote sensing images in nearshore marine scenarios. In the flare correction stage, a robust design for SWIR proxy can be achieved by utilizing the unique physical characteristics of the first shortwave infrared band in nearshore deep water regions, such as the B6 band of the Landsat 8 OLI sensor, with a spectral wavelength range of 1.560-1.660µm. The water reflectivity in this band is less than 0.01, and the atmospheric path radiation contribution is stable; therefore, its reflectivity signal is mainly dominated by specular reflection (flare), which can be used to avoid coupling interference between water reflection and flare signals. The first shortwave infrared band is used as the core carrier of the flare intensity proxy, instead of using the visible band or external parameters. Simultaneously, a reference anchor point is determined by a preset percentile within the deep-water mask, eliminating the influence of background radiation differences under different scenarios and solving the scene-dependent problem of flare subtraction.

[0015] In addition, the flare coefficients are robustly fitted using the median slope method with a fixed intercept of 0 and numerical constraints are applied, such as a preset numerical range of [0, 5]. This avoids the influence of outliers on the fitting results and allows for the selection of global scale parameters for fitting through cross-band spectral consistency constraints, thereby reducing the risk of color shift and ensuring the color fidelity of the corrected image.

[0016] In some embodiments, acquiring remote sensing data includes:

[0017] Acquire raw remote sensing data;

[0018] Multiple band data are read from the original remote sensing data using a mask array method, wherein the band data retains information on data-free values ​​and mask information;

[0019] Based on the key metadata of the source file of the original remote sensing data, the band data is cropped and resampled.

[0020] Value range processing is performed on the band data after data cropping and resampling to generate the remote sensing data.

[0021] In some embodiments, extracting a target mask from the remote sensing data includes:

[0022] Input bands are read from the remote sensing data, including green visible light band, red visible light band, near-infrared band, first shortwave infrared band and second shortwave infrared band;

[0023] The spectral index is calculated based on the reflectance value of the input band. The spectral index includes the improved normalized differential water index calculated based on the reflectance of the green visible light band and the first shortwave infrared band, and the normalized differential vegetation index calculated based on the reflectance of the near-infrared band and the red visible light band.

[0024] Water body screening conditions are set based on the spectral index;

[0025] Water body masks are extracted from the remote sensing data according to the water body screening conditions. The water body masks include loose water body masks and strict water body masks. The loose water body mask is used as a decoupling application domain. The strict water body mask is used as a sampling domain.

[0026] In some embodiments, the method further includes:

[0027] Effective values ​​are extracted from the strict water mask, the effective values ​​including a first effective value extracted in the first shortwave infrared band and a second effective value extracted in the second shortwave infrared band;

[0028] A preset quantile threshold is set based on the effective value, and the preset quantile threshold includes a first threshold and a second threshold; the preset quantile threshold is dynamically set according to the number of effective samples.

[0029] Based on the preset quantile threshold, a deep-water mask is extracted from the strict water mask. The deep-water mask is a set of pixels whose reflectance in the first shortwave infrared band is less than the first threshold and whose reflectance in the second shortwave infrared band is less than the second threshold.

[0030] Output statistical information, which includes the number of pixels and the percentage of pixels in the loose water mask, the strict water mask, and the deep water mask.

[0031] In some embodiments, solar flare correction is performed on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum, including:

[0032] Specular reflection proxy is set based on the reflectivity of the first shortwave infrared band;

[0033] Within the deep-water mask, a preset percentile of the first shortwave infrared band is counted to obtain a reference anchor point;

[0034] Based on the specular reflection proxy and the reference anchor point, a flare intensity proxy is constructed, wherein the flare intensity proxy is the non-negative part of the difference between the specular reflection proxy and the reference anchor point;

[0035] Extract the reflectance in the visible light band from the remote sensing data;

[0036] Based on the reflectivity of the visible light band and the flare intensity surrogate, the median slope method is used to robustly regress and estimate the relationship slope, and the relationship slope is subject to numerical constraints within a preset numerical range.

[0037] Based on the slope of the relationship, a first relationship model is constructed between the reflectivity of the visible light band and the flare intensity proxy.

[0038] In some embodiments, performing solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum further includes:

[0039] Obtain the relative weights, which are either default weight values ​​or set using an offline simulation table;

[0040] The fitted independent variables are calculated based on the relative weights and the flare intensity proxy.

[0041] Using the reflectivity of the visible light band as the dependent variable, a stacked regression was performed on the deep-water mask using the median slope method to estimate the global scale parameters.

[0042] Based on the global scale parameters, a second relationship model is constructed between the reflectivity of the visible light band and the flare intensity surrogate quantity.

[0043] In some embodiments, calculating the multi-band luminosity index of the deepwater mask after solar flare correction includes:

[0044] The dark image band is determined, which includes the blue visible light band, the green visible light band, and the red visible light band;

[0045] Extract the pixel values ​​of the dark image band from the deep-water mask;

[0046] For each pixel, the average value of the corresponding pixel value of the dark image band is calculated to obtain the multi-band darkness index;

[0047] Based on the multi-band luminance index, the index mean and index standard deviation of the pixel set are calculated;

[0048] Obtain a preset dark image coefficient, and calculate the theoretical dark pixel reflectance based on the dark image coefficient, the mean of the index, and the standard deviation of the index.

[0049] By employing the above technical solutions, this application provides an atmospheric correction method for remote sensing images of nearshore marine scenes. This method overcomes the limitations of dark pixel methods, which rely on single-band, single-point data. It proposes a multi-band darkness index and effectively smooths the influence of single-band random noise by jointly averaging the pixel values ​​of three visible bands. The technical approach of finding the darkest single pixel is optimized into identifying a set of dark pixels with consistent spectral characteristics, thereby improving the statistical representativeness and reliability of dark pixel signals from the perspective of technical principles.

[0050] Simultaneously, by calculating the mean and standard deviation of the multi-band darkness index of the dark pixel set, and combining it with the preset dark image coefficient k, the theoretical dark pixel reflectance is derived, and ρ is applied. toa,obs The physical constraint ≥0 is applied. The operation mode of setting thresholds based on human experience in the dark pixel method is optimized into a quantitative inference method based on statistical laws. Through clear statistical parameters and calculation logic, the subjective human intervention in the threshold selection process is reduced, effectively solving the technical defects of lack of unified standards for threshold setting and sensitivity to extreme anomalies. This provides more stable and reliable observational benchmark data for subsequent aerosol optical thickness (AOT) inversion.

[0051] In some embodiments, the optimal aerosol optical thickness value is retrieved based on the multi-band luminosity index using a radiative transfer model, including:

[0052] The geometric constraints of the radiative transfer model are set according to the auxiliary parameters;

[0053] The objective function of the radiative transfer model is constructed based on the theoretical dark pixel reflectance. The objective function is used to represent the difference between the actual observed dark pixel reflectance and the theoretical dark pixel reflectance.

[0054] Set the search interval and step size for the radiative transfer model;

[0055] According to the preset output criteria, an inversion is performed based on the radiative transfer model to determine the optimal aerosol optical thickness value; the preset output criteria are used to take the aerosol optical thickness value that minimizes the objective function as the optimal aerosol optical thickness value.

[0056] In some embodiments, the optimal aerosol optical thickness value is retrieved by inverting the radiative transfer model based on the multi-band luminosity index, and the method further includes:

[0057] The physical parameters are calculated band by band.

[0058] The apparent reflectance is extracted from the remote sensing data, and the surface reflectance of each pixel is calculated based on the apparent reflectance and the physical parameters.

[0059] The effective pixels are determined according to the surface reflectance.

[0060] The number of effective pixels and the percentage of effective pixels are statistically analyzed.

[0061] According to a second aspect of this application, an atmospheric correction system for remote sensing images of nearshore marine scenes is provided, the system comprising:

[0062] The data acquisition module is used to acquire remote sensing data, which includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light band, near-infrared band, and short-wave infrared band. The auxiliary parameters include solar geometric parameters and sea surface state parameters.

[0063] A mask extraction module is used to extract target masks from the remote sensing data, the target masks including water masks and deep-water masks;

[0064] The flare correction module is used to perform solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum. The secondary simulation model is a relationship model between the visible light band and the flare intensity proxy based on the flare coefficient. The flare intensity proxy uses a first shortwave infrared band as a specular reflection proxy, constructed in conjunction with a reference anchor point. The reference anchor point is obtained by statistically analyzing a preset percentile of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained using robust regression fitting with the median slope method.

[0065] The index calculation module is used to calculate the multi-band luminosity index of the deep-water mask after solar flare correction. The multi-band luminosity index is calculated based on the pixel value of each pixel in the blue visible light, green visible light and red visible light bands of the deep-water mask.

[0066] The inversion module is used to invert the optimal aerosol optical thickness value based on the multi-band luminosity index using a radiative transfer model; the radiative transfer model is a physical parameter calculation model based on the auxiliary parameter settings; the physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance.

[0067] The result output module is used to calculate the atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value, and to generate correction result data based on the atmospheric correction coefficients. The correction result data includes corrected images, difference maps, and color composite images.

[0068] According to a third aspect of this application, a computer device is provided, including a storage medium, a processor, and a computer program stored on the storage medium and executable on the processor, wherein the processor executes the program to implement the above-described atmospheric correction method for remote sensing images of nearshore marine scenes.

[0069] According to a fourth aspect of this application, a storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements the above-described atmospheric correction method for remote sensing images of nearshore marine scenes.

[0070] By employing the above technical solutions, this application provides a method and system for atmospheric correction of remote sensing images for nearshore marine scenes. After acquiring remote sensing data, the method extracts target masks and performs solar flare correction using a secondary simulation model of satellite signals in the solar spectrum. Then, it calculates multi-band luminosity indices and retrieves the optimal aerosol optical thickness value using a radiative transfer model. Finally, it calculates atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value to generate correction result data. This method can perform robust flare correction based on a SWIR proxy, utilizing the water reflection and path radiation physical characteristics of SWIR1 (B6) under deep-water conditions to construct a scene-adaptive flare proxy, achieving flare subtraction. Furthermore, it uses multi-band joint luminosity indices for inversion, leveraging the spectral characteristics of clean deep-water bodies in the visible spectrum reflectance to effectively smooth single-band random noise and improve the accuracy of AOT estimation results.

[0071] The above description is only an overview of the technical solution of this application. In order to better understand the technical means of this application and to implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of this application more obvious and understandable, the following are specific embodiments of this application. Attached Figure Description

[0072] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0073] Figure 1 A schematic diagram of the atmospheric correction method for remote sensing images of nearshore marine scenes provided in this application embodiment;

[0074] Figure 2 This is a schematic diagram of the overall process of atmospheric correction of remote sensing images provided in the embodiments of this application;

[0075] Figure 3 This is a schematic diagram of the solar flare correction process provided in an embodiment of this application;

[0076] Figure 4 This is a schematic diagram of an automated aerosol optical thickness estimation process provided in an embodiment of this application;

[0077] Figure 5 This is a schematic diagram of the structure of an atmospheric correction system for remote sensing images of nearshore marine scenes, provided in an embodiment of this application. Detailed Implementation

[0078] The present application will be described in detail below with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in the embodiments of the present application can be combined with each other.

[0079] In this application embodiment, remote sensing imagery refers to image data acquired by spectral sensors carried by spacecraft such as satellites. For example, in the field of marine remote sensing, satellite sensors can receive radiation signals. These radiation signals include apparent reflectance (Top-of-Atmosphere Reflectance, TOA). Apparent reflectance refers to the reflectance value obtained by normalizing the radiation signal received by the satellite sensor at the top of the atmosphere to solar irradiance.

[0080] Since the radiation signals received by satellite sensors are the superposition of multiple effects such as reflection from the ocean surface or water body, atmospheric scattering and absorption, and specular reflection from the water surface, the data quality of the radiation signals directly determines the reliability of subsequent marine environmental monitoring such as red tide identification and water quality parameter inversion, marine resource surveys such as fishery resource assessment and coastal zone mapping, and global marine carbon cycle research.

[0081] To obtain accurate ocean surface reflectance, atmospheric correction of remote sensing imagery is necessary. Atmospheric correction, a core component of ocean remote sensing data processing, removes influencing factors such as atmospheric scattering, absorption, and solar flares through correction algorithms, revealing the true reflectance of the ocean surface or water body. Accurate estimation of aerosol optical thickness (AOT) is both a prerequisite and a challenge for atmospheric correction. Meanwhile, specular reflection from the water surface significantly increases apparent reflectance in the visible band. Without effective flare correction, this not only leads to inflated reflectance and color cast but also disrupts the stability of dark pixels in deep water regions, causing systematic biases in aerosol optical thickness estimation.

[0082] In some embodiments, sea surface flare correction can be performed based on geometry- or wind-driven forward models such as Cox-Munk. This involves using solar and sensor geometry, wind speed, and wave slope statistics to predict specular reflection and subtract the flare effect. However, forward models rely on external priors and parameters such as wind speed, wave state, bidirectional reflection distribution function (BRDF), and white waves. Parameter or prior biases can be amplified by the observed geometry, easily leading to under-correction or over-correction problems. Furthermore, the implementation of geometry-driven forward model correction methods is complex, heavily reliant on computational resources and external data, making them unsuitable for large-scale, online, or cross-sensor stable deployment.

[0083] In some embodiments, sea surface flare correction can be performed using band-by-band empirical regression. This involves fitting a linear relationship to a flare surrogate quantity separately for each visible band to correct the flare. However, empirical regression is sensitive to outliers and thin clouds, its slope is easily manipulated, and independent fitting of each band leads to cross-band inconsistencies and color shifts. If the intercept is allowed, the fitting may also swallow scene-specific additive biases, making the results scene-dependent or unstable.

[0084] In some embodiments, the dark pixel method can also be used for aerosol optical thickness estimation. When estimating aerosol optical thickness, the pixel with the lowest single-band reflectance, such as the blue light band B2, in the image can be found, and its observed value can be directly equated to atmospheric path radiation, thus achieving atmospheric path radiation estimation based on the assumption of a darkest point.

[0085] However, the dark pixel method produces unstable outputs in marine scenes. A single pixel may exhibit extremely low values ​​due to residual noise, sub-pixel level foam, thin cloud pollution, or transient sensor anomalies. This darkest point is likely an atypical outlier and does not represent the true atmospheric signal of the entire clean deep-water area. Furthermore, because the algorithm uses fixed low quantiles such as the 1% or 5% quantile as observations, the selection is highly empirical and arbitrary, and it is very sensitive to a small number of extreme high-value outliers in the data. This can easily lead to large fluctuations in the AOT estimation results, lacking statistical robustness and resulting in inaccurate AOT estimation results.

[0086] Compared to open ocean, nearshore marine scenes present unique challenges such as water turbidity, uneven flare distribution, and frequent sub-pixel foam and thin cloud pollution, significantly reducing the accuracy and robustness of atmospheric correction methods. To address the issue of inaccurate AOT estimation results, this application provides an atmospheric correction method for remote sensing images in nearshore marine scenes in some embodiments. This method can perform atmospheric correction on remote sensing images for marine or nearshore water environments. Targeting the core challenges of complex nearshore scenes, the method achieves dual optimization of AOT estimation and flare correction through band physical characteristic mining, robust algorithm design, and hierarchical mask collaboration, providing high-precision data support for nearshore marine environmental monitoring and water quality parameter inversion.

[0087] The method may include key steps such as water body and deep-water mask extraction, flare correction using a first short-wave infrared (SWIR1) proxy, dark pixel-driven AOT estimation, calculation and inversion of surface reflectance using the second simulation of the satellite signal in the solar spectra (6S), difference analysis, and RGB synthesis. Based on deep-water dark pixel constraints, visible light multi-band flare removal, and physical radiative transfer simulation, the method achieves water color observation reflectance correction for multispectral sensors such as Landsat 8 OLI, and outputs correction results, difference statistics, and visualization products.

[0088] The method can be applied to electronic devices with data processing capabilities. These electronic devices include, but are not limited to, computers, servers, mobile terminals, smart wearable devices, and industrial control machines. For ease of description, this application embodiment uses an electronic device as the execution subject of the method. It should be understood that the method can also be applied to other types of execution subjects, which are not illustrated in this application embodiment. Figure 1 As shown, the method includes:

[0089] S101. Acquire remote sensing data.

[0090] When performing atmospheric correction on remote sensing imagery, remote sensing data can be acquired first. This data includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light, near-infrared, and short-wave infrared bands.

[0091] For example, using the Landsat 8 OLI satellite sensor, remote sensing images covering reflectance in multiple bands, including visible light, near-infrared, and SWIR, can be acquired. These multiple bands can be represented using B1–B7. Among them, the B1 band refers to the coastal aerosol band, with a spectral wavelength range of 0.433–0.453µm; the B2 band refers to the blue visible light band, with a spectral wavelength range of 0.450–0.515µm; the B3 band refers to the green visible light band, with a spectral wavelength range of 0.525–0.600µm; the B4 band refers to the red visible light band, with a spectral wavelength range of 0.630–0.680µm; the B5 band refers to the near-infrared band, with a spectral wavelength range of 0.845–0.885µm; the B6 band refers to the first shortwave infrared (SWIR1) band, with a spectral wavelength range of 1.560–1.660µm; and the B7 band refers to the second shortwave infrared (SWIR2) band, with a spectral wavelength range of 2.100–2.300µm.

[0092] Remote sensing imagery can be created by receiving reflectance data across multiple wavelengths using satellite sensors. Electronic devices can then acquire this remote sensing imagery by sending data acquisition requests to the satellite. To facilitate data transmission, remote sensing imagery data can be formatted into specific file formats. For example, a remote sensing image might include TOA reflectance ρ... toa The image data is converted into an image file in the format LC08_L1TP_…_TOA_B1B7.tif.

[0093] While acquiring remote sensing imagery, auxiliary parameters can also be obtained. These auxiliary parameters can include solar geometric parameters and sea surface state parameters. Solar geometric parameters include runtime interactive input parameters such as solar zenith angle and solar azimuth angle. Sea surface state parameters can include runtime interactive input parameters such as wind speed, wind direction, salinity, and chlorophyll concentration. Auxiliary parameters can be obtained by calculating the satellite's current position and the geographic information system of the measured area.

[0094] In some embodiments, data preprocessing can be performed during the data acquisition process to acquire remote sensing data, i.e., ... Figure 2 As shown, the process involves acquiring raw remote sensing data and then reading multiple band data from it using a mask array. The band data retains information about missing data values ​​and masking information. Based on the key metadata of the original remote sensing data source file, the band data is then cropped and resampled. Finally, value range processing is performed on the cropped and resampled band data to generate the final remote sensing data.

[0095] For example, in the preprocessing of remote sensing data, image reading can be performed first, that is, reading 7 bands using a mask array method, retaining no-data information and masking information. Then, data cropping and resampling are performed. During data cropping and resampling, the data can be written back in the source file format (profile) without changing the spatial reference and resolution. Finally, value range processing is used to adjust the input TOA reflectance ρ. toa Reasonable truncation is performed to limit the range of reflectivity to the [0, 1] interval, thereby enhancing robustness.

[0096] S102. Extract the target mask from the remote sensing data.

[0097] After acquiring remote sensing data, water and deepmask computation can be performed to extract target masks from the remote sensing data. These target masks include water masks and deep masks.

[0098] In some embodiments, when extracting a target mask from remote sensing data, an input band can be read from the remote sensing data first, wherein the input band includes a green visible light band, a red visible light band, a near-infrared band, a first shortwave infrared band, and a second shortwave infrared band.

[0099] For example, when performing water body and deep-water mask calculations, data input and preprocessing can be performed first. The electronic device can then safely read the input bands from the remote sensing data by executing the "get band data" command. Input bands can include the B3 (green) band, B4 (red) band, B5 (NIR) band, B6 (SWIR1) band, and B7 (SWIR2) band, and automatically handle masked arrays and invalid values ​​(Not a Number, NaN).

[0100] After reading the input band, spectral indices can be calculated based on the reflectance values ​​of the input band. These spectral indices include the improved normalized difference water index, calculated based on the reflectance of the green visible light band and the first shortwave infrared band, and the normalized difference vegetation index, calculated based on the reflectance of the near-infrared band and the red visible light band.

[0101] For example, to maintain the stability of values ​​in remote sensing data, when calculating spectral indices, you can first use the `np.errstate(divide='ignore', invalid='ignore')` command, and then use `np.nan_to_num` to set NaN and ±Inf to 0 to avoid overflow affecting the mask.

[0102] Next, the indices are defined as follows: the Modified Normalized Difference Water Index (MNDWI) and the Normalized Difference Vegetation Index (NDVI). The MNDWI and NDVI can be calculated using the following formulas:

[0103] MNDWI=(B3-B6) / (B3+B6);

[0104] NDVI=(B5-B4) / (B5+B4);

[0105] Among them, MNDWI represents the Improved Normalized Difference Water Index; B3 represents the reflectance in the green visible light band; B6 represents the reflectance in the first shortwave infrared band; NDVI represents the Normalized Difference Vegetation Index; B5 represents the reflectance in the near-infrared band; and B4 represents the reflectance in the red visible light band.

[0106] After calculating the spectral index, water body screening conditions can be set based on the spectral index, thereby extracting water body masks from remote sensing data according to the water body screening conditions. These water body masks include loose water body masks and strict water body masks; loose water body masks are used as decoupling application domains, while strict water body masks are used as sampling domains.

[0107] For example, after calculating MNDWI and NDVI, water body screening conditions can be set based on MNDWI and NDVI to perform initial water body screening. Since typical open water bodies are relatively bright in the green band and relatively dark in the SWIR band, while water bodies are dark in the NIR band and vegetation is positive, the screening conditions can include: Condition 1: MNDWI > 0; Condition 2: NDVI < 0. Combining these conditions, the water body mask screening conditions can be obtained as: water_mask0 = (MNDWI > 0) or (NDVI < 0).

[0108] To accommodate the subsequent decoupling process, two levels of masks can be set up, serving as the decoupling application domain and the sampling domain, respectively. The two levels of masks include a loose water mask (application domain), denoted as `water_mask_app`. The selection criteria for the loose water mask are: `water_mask_app = (MNDWI > -0.2)` or `(NDVI < 0)`. The loose water mask can be used without superimposing the SWIR upper limit to ensure that strong flare surfaces are not excluded.

[0109] In some embodiments, terrestrial or cloud shadow protection can also be implemented, i.e., by overlaying external cloud, shadow, and snow quality assessments (QA) for exclusion. The Normalized Difference Snow Index (NDSI) and NDVI upper limit (e.g., NDVI ≥ 0.3) can also be used to label terrestrial protected areas. At the morphological level, opening operations (one erosion followed by one expansion) can be used to reduce noise and maintain narrow waterway connectivity. Area filtering is then applied to remove small patches smaller than 100 pixels.

[0110] For a strict water mask (sampling area), it can be represented as `water_mask_strict`. The screening criteria for a strict water mask are defined as: `water_mask_strict = water_mask0`; (B6 < 0.03) and (B7 < 0.02). The thresholds used are empirical values ​​for OLI TOA, which can be relaxed to 0.035 or 0.025 depending on the nearshore turbidity, and can also be adaptively set. At the morphological level, opening operations and area filtering are also employed to avoid systematic shrinkage caused solely by corrosion.

[0111] After screening out the water body mask according to the screening conditions, a deep water mask can be further extracted based on the water body mask. In some embodiments, the electronic device first extracts valid values from the strict water body mask, where the valid values include a first valid value extracted in a first short-wave infrared band and a second valid value extracted in a second short-wave infrared band.

[0112] Then, a preset quantile threshold is set according to the valid values, and a deep water mask is extracted from the strict water body mask based on the preset quantile threshold. The preset quantile threshold includes a first threshold and a second threshold; the deep water mask is a set of pixels whose reflectance in the first short-wave infrared band is less than the first threshold and whose reflectance in the second short-wave infrared band is less than the second threshold.

[0113] For example, since the reflectance of deep water in the SWIR band is extremely low and stable, the deep water mask can be extracted according to the SWIR band. By extracting the valid values of the B6 and B7 bands within water_mask_strict and then calculating the 20th quantile threshold respectively, the first threshold B6_20th and the second threshold B7_20th can be obtained. Then, by defining the screening conditions for the deep water mask, that is: deep_mask = water_mask_strict; (B6 < B6_20th) and (B7 < B7_20th), the extraction of the deep water mask is achieved.

[0114] To improve the robustness of the extracted parameters, the preset quantile threshold can be dynamically set according to the number of valid samples, that is, an adaptive SWIR threshold is performed to replace the fixed thresholds of 0.03 or 0.02. Then, the brightness upper limit is determined by quantiles within water_mask0, such as B6 < P80 and B7 < P80 as the SWIR constraints of water_mask_strict, so as to better fit the specific application scenario.

[0115] In some embodiments, if the number of samples of the remote sensing data is insufficient, a fallback operation that supports dynamic threshold adjustment can be performed. That is, when the number of valid samples in water_mask_strict is less than 1000, or the number of valid samples in deep_mask is less than 200, it can be determined that the number of valid samples is insufficient. At this time, the quantile can be reduced, that is, from P20 to P15 or P10, and the SWIR threshold can also be relaxed. When the number of valid samples is extremely insufficient, it can fallback to the fixed threshold loading and record an alarm.

[0116] After extracting the target mask from the remote sensing data, quality control and statistical information output can also be performed. That is, the electronic device can output statistical information, where the statistical information includes the number of pixels and the pixel ratio of the loose water body mask, the strict water body mask, and the deep water mask.

[0117] For example, when generating statistical output, the electronic device can print the number of pixels and the percentage of pixels for water_mask_app, water_mask_strict, and deep_mask, respectively, to form a sample statistical report. This sample statistical report is used to determine the final sample size for AOT or flare regression.

[0118] In the downstream explicit decoupling process, the strict deep-water sample (deep_mask) can be used for AOT estimation and flare coefficient regression to achieve robust fitting of visible light and flare intensity surrogate X. The relaxed application domain (water_mask_app) is used to apply flare correction to a wider water body area containing strong flares, but not to land areas or cloud shadows.

[0119] In some embodiments, for turbid nearshore phenomena, missed detections occur due to MNDWI > 0 but high B6 / B7 ratios. Since the application domain depends on water_mask_app to cover strong flares, the sampling domain can be relaxed for SWIR or adaptive quantiles can be used. For shallow sandy bottom phenomena, the visible light band is bright, but the SWIR band is still relatively dark, easily leading to contamination issues. Therefore, low quantiles B2 or B3 or shoreline constraints can be superimposed on the deep_mask. For cloud shadows or shading, setting NDVI < 0 may result in misjudgment of water bodies; therefore, a QA cloud shadow mask or shading mask can be superimposed.

[0120] Based on the content described in the above embodiments, by designing a hierarchical mask system including a loose water body mask, a strict water body mask, and a deep water mask, a clear separation between the calibration application domain and the parameter sampling domain can be achieved. Specifically, for the loose water body mask, i.e., the calibration application domain, a screening condition of improved normalized difference water index (MNDWI) > -0.2 or normalized difference vegetation index (NDVI) < 0 can be used, without additional constraint on the upper limit of shortwave infrared (SWIR) band reflectivity. This effectively avoids the erroneous rejection of areas covered by strong flares on the water surface and at the water body edges. Actual measurements have verified that its coverage of the target water body can reach over 95% of the total measured water area.

[0121] For strict water body masking, i.e. parameter sampling domain, the basic screening condition of "MNDWI > 0 or NDVI < 0" can be combined with the SWIR band reflectance threshold for secondary purification, i.e., the first shortwave infrared band B6 < 0.03 and the second shortwave infrared band B7 < 0.02, to ensure that the purity of the screened water body sample meets the requirements for core parameter estimation.

[0122] For the deep-water mask, i.e., the core sample domain, an adaptive quantile threshold can be used for extraction based on SWIR band reflectivity, with the quantile dynamically adapting to the effective sample size. When the effective sample size is ≥500, the 20th percentile (P20) is used; when 200 ≤ effective sample size < 500, the 15th percentile (P15) is used; and when 100 ≤ effective sample size < 200, the 10th percentile (P10) is used. The resulting deep-water sample purity is ≥90%.

[0123] Comparative tests were conducted using remote sensing images of turbid, complex nearshore waters. A single-body masking method was used as the baseline. Regarding parameter estimation accuracy, the hierarchical masking system increased the coefficient of determination (R²) for flare coefficient regression from 0.65 to 0.88, reduced the aerosol optical thickness (AOT) estimation error by over 40%, and significantly improved the stability and accuracy of core parameter estimation. In terms of noise suppression, morphological opening operations (one erosion followed by one expansion) and area filtering (removing isolated small patches with an area <100 pixels) effectively suppressed the impact of isolated noise points, land debris, and other interference factors on sample quality, providing highly reliable sample support for subsequent solar flare correction and AOT inversion.

[0124] Comparative experimental results show that the hierarchical mask system not only ensures the coverage of solar flare correction in complex nearshore scenarios, but also provides high-purity samples for core parameter estimation. It resolves the contradiction between sample mixing leading to parameter estimation bias and incomplete coverage leading to omission of correction areas in traditional single-body masks, and significantly improves the adaptability and robustness of atmospheric correction technology in complex nearshore scenarios.

[0125] S103. Perform solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum.

[0126] After extracting the target mask from the remote sensing data, solar flare removal can be performed based on SWIR1. Solar flare correction can be performed on the remote sensing image corresponding to the target mask using the secondary simulation (6S) model of satellite signals in the solar spectrum.

[0127] The secondary simulation model is a model relating visible light bands and flare intensity surrogate quantities based on the flare coefficient. The flare intensity surrogate quantity is obtained by using the first shortwave infrared band as a specular reflection surrogate and combining it with reference anchor points. The reference anchor points are obtained by statistically analyzing the preset percentiles of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained by robust regression fitting using the median slope method.

[0128] Since SWIR1 (B6) contributes weakly to water reflection and has a significant response to sunlint, it can be used as a proxy variable for flare. Therefore, in the deep water region, robust regression (Theil–Sen) can be used to fit the linear relationship between the visible light band and (B6-B6 deep water quantile) to estimate the flare coefficient k of each visible band, and then remove it in the water region according to glint=k·max(B6-proxy_min,0).

[0129] To perform solar flare correction, in some embodiments, when performing solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum, a specular reflection proxy can first be set based on the reflectivity of the first shortwave infrared band, and a preset percentile of the first shortwave infrared band can be statistically analyzed within the deep-water mask to obtain a reference anchor point. Then, based on the specular reflection proxy and the reference anchor point, a flare intensity proxy is constructed. The flare intensity proxy is the non-negative portion of the difference between the specular reflection proxy and the reference anchor point.

[0130] Next, the reflectance in the visible light band is extracted from the remote sensing data. Then, based on the reflectance in the visible light band and the surrogate magnitude of the flare intensity, a robust regression estimation of the relationship slope is performed using the median slope method, wherein the relationship slope is numerically constrained within a preset range. Thus, based on the relationship slope, a first relationship model between the reflectance in the visible light band and the surrogate magnitude of the flare intensity is constructed.

[0131] For example, such as Figure 3 As shown, when performing solar flare correction in the visible light band, a proxy band and reference anchor point can be set first. For the proxy band, since the outgoing radiation from clean deep water is approximately zero in SWIR, SWIR mainly carries specular reflection and path radiation. Therefore, SWIR1 (B6) can be used as a specular reflection proxy. As for the reference anchor point (proxy_min), the 5th percentile (P5) of B6 within the deep water mask can be used as the reference anchor point proxy_min, which approximately represents a background without flares or with weak flares, thus eliminating common bias.

[0132] After setting the proxy band and reference anchor point, the flare intensity proxy can be constructed. That is, by defining the flare intensity proxy X=max(B6-proxy_min,0), it represents the non-negative part of the difference between the specular reflection proxy and the reference anchor point, so as to avoid negative values ​​and only model the specular component that exceeds the background.

[0133] By selecting and cleaning samples, the sample domain is determined, namely, the deep_mask from the strict water mask water_mask_strict, to ensure that the water outgoing radiation is weakest in SWIR and the regression relationship is more stable.

[0134] When determining the sample domain, anomaly removal can be performed, i.e., removing extreme values ​​in deep water, such as X > 95th percentile, or any band of B1–B4 > 99th percentile, to reduce the influence of anomalies such as foam, ship tracks, and cloud shadows. Furthermore, when the sample size is too large, downsampling can be used, randomly downsampling to no more than 20,000 sample points to improve computational efficiency and reduce the risk of overfitting.

[0135] When performing robust regression estimation, the target relationship can be determined first. That is, for each visible band Bi (i∈{1,2,3,4}), the first relationship model is determined by fitting the linear relationship between the visible band Bi and the flare intensity surrogate X. The linear relationship of the first relationship model is expressed as:

[0136] Bi≈k i X+b i ;

[0137] in, k i Indicates the slope of the relationship; b i This represents the relation intercept. To reduce the bias caused by the instability of the constant term, and because the common bias has already been approximated by using the reference anchor point `proxy_min`, the relation intercept can be fixed. b i ≈0, in practical applications only the slope of the relationship is used. k i .

[0138] Therefore, the slope of the relationship can be estimated using regression methods. k i Because simple linear regression requires strong constraints and outlier removal, its robustness is insufficient. Therefore, the median slope method (Theil–Sen) can be used for robust regression estimation. k i To ensure robustness against outliers.

[0139] You can also limit the amplitude using parameters for each... k i Numerical constraints are applied in the range of [0, 5] to avoid overcorrection due to residual anomalies or sample bias. Tighter upper limits can also be set according to the band, such as smaller for blue and green and larger for red and near-red. The specific limit value can be determined according to the sensor and sea conditions.

[0140] In some embodiments, solar flare correction can also be performed through robust regression estimation with spectral consistency constraints. This involves first obtaining relative weights, where the relative weights are default values ​​or set using an offline simulation table. Fitting independent variables are then calculated based on the relative weights and flare intensity surrogates. Next, with visible light reflectance as the dependent variable, stacked regression is performed on a deep-water mask using the median slope method to estimate global scale parameters. Finally, a second relationship model between visible light reflectance and flare intensity surrogates is constructed based on these global scale parameters.

[0141] For example, in the robust regression estimation process with spectral consistency constraints, cross-band spectral consistency constraints can be introduced, estimating only a global scale parameter kg. Then, by determining the target relation, for Bi (i∈{1,2,3,4}), the second relation model can be expressed as:

[0142] Bi≈k g ( c i X ) + 0;

[0143] Where, k g Represents the global scale parameter; c i This represents the relative weight. Based on the near-gray assumption, the relative weight defaults to c. i =1, or c can be set using an offline simulation table. i .

[0144] Then, regression analysis can be used to stack regressions y=concat(B1, B2, B3, B4) and x=concat(c1X, c2X, c3X, c4X) on deep-water samples, and the global scale parameter k can be estimated using the median slope method (Theil–Sen). g And the forced intercept is 0.

[0145] Similarly, for global scale parameters, parameter limiting can also be used to limit the global scale parameter k. g The constraint is set within the interval [0, 5] to avoid overcorrection caused by residual anomalies. Furthermore, the relative weight c... i It can be a fixed value or a weak constraint, such as [1, 1, 1, 1] etc.

[0146] like Figure 4 As shown, in some embodiments, a constraint strategy can also be set to adjust the relationship slope based on sample sufficiency. That is, if the number of deep-water samples is less than 50, or the proportion of valid samples (X>0) is too low (e.g., <20 samples), then k can be set. i=0 to skip flare correction for that band. The flare correction process can also be set based on mask integrity. If the water or deep-water mask is empty or has a very small area, the entire flare correction process can be skipped, and a quality marker recorded.

[0147] Following the secondary simulation model shown in the above embodiment, after determining the relational model, solar flare correction can be performed using the relational model and applied to the water body area. The application domain mask_app can then include: water_mask_app, effective pixels, and non-cloud or shadow masks. Correction can then be performed according to the following correction formula based on the relational model:

[0148] ρ corri =ρ i -k i max(B6-proxy_min, 0);

[0149] Where, ρ corri ρ represents the corrected reflectance of the i-th band, i∈{1,2,3,4}; i k represents the raw reflectance, i.e., the apparent reflectance in remote sensing data. i B6 represents the slope of the relationship; B6 represents the reflectivity of the first shortwave infrared band; and proxy_min represents the reference anchor point.

[0150] After solar flare correction, solar flare phenomena in remote sensing images can be removed, and SWIR B6 is retained as a proxy without flare correction, retaining the original value. Since the near-infrared (B5) band is not a visible light band, its effect on remote sensing images is relatively weak. Therefore, the near-infrared (B5) band does not need to be corrected or should be evaluated separately before deciding whether to correct it.

[0151] The solar flare correction method described in the above embodiments utilizes the physical characteristics of SWIR1 (B6) under deep-water conditions—near-zero water reflection and primarily bearing specular reflection and path radiation—to construct a scene-adaptive flare proxy and perform unified and robust subtraction. By taking the low quantile P5 of B6 in deep water as the scene baseline, a non-negative flare intensity X = max(B6 - proxy_min, 0) is defined. A common slope across the visible band is estimated using a robust regression method passing through the origin on the deep-water sample, and then subtracted from each visible band according to the second relationship model.

[0152] The method guarantees physical consistency and non-negativity constraints, ensuring that zero flare levels are not deducted and that brightness is not mistakenly increased. By sharing slopes across bands and supplementing with a slight spectral weight ci, color shift and banding differences are significantly reduced, maintaining color fidelity. Furthermore, robustness is improved by suppressing outlier effects through quantile baselines and robust regression. Moreover, it achieves complete scene adaptation without external priors; when there are no significant flares, X is near zero, and no correction action is required. Therefore, the computation is simple, efficient, and easily parallelized and processed in blocks, making it suitable for large-scale and on-orbit or online applications. In addition, it can be coupled with band-by-band slope consistency and residual quality control and backoff strategies to improve diagnostics and engineering usability.

[0153] To validate the flare correction scheme, three typical nearshore marine scenarios were selected: Landsat 8 OLI images of open waters with strong flares, turbid nearshore waters, and shallow sandy-bottomed waters, for comparative experiments. Mainstream techniques such as the Cox-Munk forward model and band-by-band empirical regression were used as controls. Experimental comparison results were obtained through testing indicators such as flare correction residuals, cross-band color shift index, and computation time.

[0154] Regarding the flare correction residuals, the residuals of visible band (B2-B4) flares corrected by the method are all ≤0.01 (units of reflectivity); while the Cox-Munk model has residuals of 0.02 to 0.05 when no high-precision wind speed data is available; and the residuals of the band-by-band empirical regression method are 0.03 to 0.06.

[0155] Regarding the cross-band color cast index, the corrected color cast index of the method is ≤0.03, ensuring the fidelity of image colors; while the color cast index of the band-by-band empirical regression method is 0.08 to 0.12, and the Cox-Munk model is 0.07 to 0.15 when wind speed data is insufficient.

[0156] Regarding computation time, for a single 1024×1024 pixel image, the computation time of the method is ≤5 minutes; while the Cox–Munk model, due to its reliance on external parameter calculation and complex geometric simulation, takes more than 30 minutes; and the band-by-band empirical regression method takes about 15 minutes.

[0157] Experimental results show that the proposed method does not rely on external prior parameters such as wind speed and wave state. It can achieve accurate flare subtraction using only the image's own band data and robust regression algorithm. It is significantly better than other schemes in terms of correction accuracy, color fidelity and computational efficiency. It is especially suitable for practical applications where nearshore scene data is scarce and large-scale rapid processing is required.

[0158] S104. Calculate the multi-band luminosity index of the deep-water mask after solar flare correction.

[0159] After performing solar flare correction, AOT estimation with deep-water dark pixel constraints can be performed based on the deep-water mask after solar flare correction. To do this, multi-band occultation indices of the deep-water mask after solar flare correction can be calculated. The multi-band occultation indices are calculated based on the pixel values ​​of each pixel in the blue visible light, green visible light, and red visible light bands of the deep-water mask.

[0160] Assuming that the surface reflectivity in the visible band of offshore deep water is extremely low, approaching that of dark pixels, its TOA observations can be used as an approximation of atmospheric contributions. By matching the minimum difference with forward simulations under different AOTs in 6S, the scene AOT can be calculated inversely. Therefore, the ocean boundary conditions corresponding to a uniform sea surface model and the atmospheric scenario of marine aerosols can be used to make the AOT estimation more closely resemble the ocean scene.

[0161] In some embodiments, when calculating the multi-band octane index of a deep-water mask after solar flare correction, the dark image bands can be determined first. These dark image bands include the blue visible light band, the green visible light band, and the red visible light band. Pixel values ​​of the dark image bands are then extracted from the deep-water mask. For each pixel, the average value of the corresponding pixel value in the dark image band is calculated to obtain the multi-band octane index. Based on the multi-band octane index, the index mean and index standard deviation of the pixel set are then calculated. A preset dark image coefficient is then obtained, and the theoretical dark pixel reflectance is calculated based on the dark image coefficient, the index mean, and the index standard deviation.

[0162] For example, when using dark pixels and 6S forward matching for AOT estimation, the observation statistics of dark pixels can be performed first. Since the overall reflectance of deep water in the corresponding spectral ranges of B2 (blue), B3 (green), and B4 (red) bands is extremely low, dark pixel statistics can be performed by selecting bands and using the three visible light bands, namely B2 (blue), B3 (green), and B4 (red). By combining multiple bands, single-band anomalies caused by noise, foam, and thin clouds can be suppressed.

[0163] For the pixel values ​​of bands B2, B3, and B4 within the deep_mask of the sample source, calculate the multi-band darkness index for each pixel, i.e.:

[0164] D = (B2 + B3 + B4) / 3;

[0165] Where D represents the multi-band luminance index; B2 represents the reflectance of the blue visible light band; B3 represents the reflectance of the green visible light band; and B4 represents the reflectance of the red visible light band.

[0166] Then, through robust statistics, using a statistical model based on probability distribution, the mean μD and standard deviation σD of the pixel set D are calculated, and then... toa，obs The theoretical dark pixel reflectance is obtained by calculating μD - k·σD. It can be seen that by calculating multi-band luminosity indices, deviations caused by high-value anomalies such as minor noise, thin clouds, and residual flares can be effectively suppressed, making the calculation results more stable and reliable. k defaults to 2 and can be adjusted within the range of 1.5–2.5. Furthermore, physical constraints can be applied during the calculation of the theoretical dark pixel reflectance, i.e., ρ toa，obs =max(ρ toa ,0).

[0167] Similarly, the calculation process can be determined based on sample sufficiency. That is, if the number of pixels in deep_mask is less than 10, the default AOT=0.08 is rolled back and a warning is issued. Alternatively, extreme value processing (winsorizing) or removing outliers in the top 1% can be performed on the samples before calculating quantiles.

[0168] S105. The optimal aerosol optical thickness value is obtained by inverting the radiative transfer model based on the multi-band octane index.

[0169] After calculating the multi-band luminosity index, 6S coefficient-driven physical atmospheric correction can be performed, that is, the optimal aerosol optical thickness value can be retrieved based on the multi-band luminosity index using a radiative transfer model. The radiative transfer model is a physical parameter calculation model based on auxiliary parameter settings; the physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance.

[0170] For example, if the radiative transfer model is a Py6S model, then the key physical parameters for each band, including the path reflectivity ρ, are calculated using Py6S. path Spherical albedo S, Total downward transmittance T d Total Upward Transmittance T u Then, the surface reflectance is inverted using the following formula, which takes into account both atmospheric scattering and multiple reflection terms during the inversion process:

[0171]

[0172] in, This represents the surface reflectance of each pixel; Indicates apparent reflectance; Indicates atmospheric path reflectivity; Indicates the total downlink transmittance; Indicates the total upward transmittance; This represents the albedo of the balloon surface.

[0173] Therefore, in some embodiments, when retrieving the optimal aerosol optical thickness value based on a multi-band luminosity index using a radiative transfer model, the geometric constraints of the radiative transfer model can be set first according to auxiliary parameters, and then the objective function of the radiative transfer model can be constructed based on the theoretical dark pixel reflectance. The objective function is used to represent the difference between the actual observed dark pixel reflectance and the theoretical dark pixel reflectance.

[0174] Then, the search interval and step size of the radiative transfer model are set, and an inversion is performed based on the radiative transfer model according to a preset output criterion to determine the optimal aerosol optical thickness value. The preset output criterion is used to select the aerosol optical thickness value that minimizes the objective function as the optimal aerosol optical thickness value.

[0175] For example, in determining the optimal aerosol optical thickness value, a 6S forward modeling configuration can be performed first, i.e., configuring geometric constraints, including solar geometry, viewing angle, altitude, atmosphere and aerosols, and sea surface or seawater. Solar geometry can be configured based on user-inputted information such as `solar_zenith` and `solar_azimuth`. The viewing angle can be approximated by `nadir(view_z=0, view_a=0)`, using OLI to approximate the vertical viewing angle sufficiently. Altitude can be set to equal the sensor altitude and the target altitude to equal the sea level, i.e., `Altitudes.set_sensor_satellite_level / set_target_sea_level`. Atmospheric and aerosol constraints are represented as `atmos_profile=Midlatitude Summer`, which can be replaced by season or latitude. `aero_profile=Maritime`. Sea surface or seawater constraints can be represented as `GroundReflectance.Lambertian`, which assumes the remaining surface approximates a Lambertian surface after flare removal.

[0176] Then, set the objective function and search strategy, where the objective function can be expressed as:

[0177] J(AOT) = |ρ toa，sim (AOT)-ρ toa，obs |;

[0178] Where J(AOT) represents the output value of the objective function; ρ toa，sim (AOT) represents the actual observed dark pixel reflectance; ρ toa，obs This represents the theoretical dark pixel reflectance.

[0179] By setting the search range and step size, i.e., AOT_RANGE=[0.01, 0.8]; AOT_STEP=0.01, and selecting the criterion, the AOT with the smallest J is taken as the optimal aerosol optical thickness value best_aot; if multiple solutions occur, the smallest AOT is selected, and the minimum difference is printed.

[0180] Similarly, when inverting the optimal aerosol optical thickness value through the radiative transfer model, a failure handling mechanism can also be set. That is, if some AOTs fail and throw errors due to 6S, the point can be skipped and the search can continue; if all fail or there are no valid points, the default AOT can be rolled back.

[0181] In some embodiments, when retrieving the optimal aerosol optical thickness value based on a multi-band octane index using a radiative transfer model, physical parameters can be calculated band by band, and apparent reflectance can be extracted from remote sensing data. Furthermore, the surface reflectance of each pixel can be calculated based on the apparent reflectance and physical parameters. Effective pixels are then determined according to the surface reflectance, and the number and percentage of effective pixels are statistically analyzed.

[0182] For example, during the 6S coefficient calculation and atmospheric correction process, the 6S settings can be configured first, including setting geometric parameters such as the user's zenith, azimuth sun angle, and viewing angle to nadir (0°). Then, altitude settings can be configured, setting the sensor to satellite altitude and the target to sea level. Next, atmospheric and aerosol settings can be configured, specifying mid-latitude summer and marine aerosols. The marine surface is represented by a Lambertian surface. Aerosol optical thickness is calculated using best_aot from the inversion process.

[0183] The coefficients are then output by calculating ρpath, S, Td, and Tu band by band. Next, the surface reflectance ρs is calculated for each pixel according to the inversion formula, with invalid pixels or those with excessively small denominators set to nodata. Finally, effectiveness statistics are performed, calculating the proportion of effective pixels band by band to facilitate quality assessment.

[0184] Based on the AOT estimation method shown in the above embodiments, a fundamental shift from finding the darkest point to identifying the darkest area, and from empirical thresholds to statistical inference, is achieved by introducing a multi-band luminosity index D and a robust statistical model based on probability distribution. By calculating the overall luminosity index D in the visible light range for each deep-water pixel, the spectral characteristics of clean deep water bodies with extremely low reflectance across the entire visible spectrum are utilized to effectively smooth single-band random noise and lock the target to dark areas with consistent spectral morphology. Furthermore, a robust statistical model is applied to the set of D values ​​for all pixels in this dark area to estimate representative samples.

[0185] Therefore, the method possesses strong anti-interference capabilities, significantly reducing the impact of single-point anomalies on the results through multi-band averaging and statistics based on the overall distribution. Furthermore, the estimated values ​​are derived based on probability distributions, and the parameter k has clear statistical significance, reducing the influence of human subjectivity and making the estimation process more rigorous and reliable. It also enhances physical consistency; multi-band joint constraints ensure that the selected signal better matches the spectral characteristics of real water bodies, thereby improving the accuracy and applicability of AOT inversion in complex scenarios.

[0186] Using measured AOT data from three AERONET sites in the nearshore area, with an observation time deviation of ≤1 hour from the image transit time, the accuracy of AOT inversion was verified. The dark pixel method was used, with the darkest point in a single blue light band as a control, to test the absolute error, stability, and anomaly suppression capabilities of AOT estimation.

[0187] Regarding the absolute error of estimation, the proposed method achieves an absolute error of ≤0.05 for AOT estimation, while the dark pixel method exhibits an error of 0.08 to 0.15 in nearshore scenes. Furthermore, in terms of stability, AOT inversion was performed on 10 consecutive nearshore images, and the standard deviation of the proposed method's estimation results was ≤0.03; while the dark pixel method, due to noise and thin cloud interference, showed a standard deviation of 0.06 to 0.12. Regarding outlier suppression capability, when thin cloud contamination is present in the images, the proposed method's AOT estimation error remains ≤0.06, without significant amplification; while the dark pixel method's error increases to 0.15 to 0.20, highlighting its sensitivity to outliers.

[0188] The verification results confirm that by using multi-band luminosity indices and robust statistical models, single-band random noise can be effectively smoothed, avoiding bias caused by single darkest point anomalies, significantly improving the accuracy and stability of AOT inversion, and providing reliable core parameter support for subsequent atmospheric correction.

[0189] S106. Calculate the atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value, and generate correction result data based on the atmospheric correction coefficients.

[0190] After estimating the optimal aerosol optical thickness value, mask robustness and visualization evaluation can be performed. This involves calculating atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value, and generating correction result data based on these coefficients. The correction result data includes corrected images, difference maps, and color composite images.

[0191] For example, by combining MNDWI and NDVI with a SWIR threshold, water masks and deep-water masks are constructed, and then morphological operations are used to suppress isolated noise. GeoTIFFs and heatmaps showing the differences before and after correction (B2–B5) are then saved for comparison and quality control. Simultaneously, true-color PNGs (B4–B3–B2) are output for visual demonstration of the effect.

[0192] The corrected image can be output as a 32-bit floating-point, LZW compressed image with nodata=-9999, or a 7-band GeoTIFF image using the original nodata. The difference image can include GeoTIFF values ​​of the difference between corrected and uncorrected values ​​for B2, B3, B4, and B5, outputting the mean, standard deviation, extreme values, and effective pixel count for each band. The difference heatmap is symmetrically stretched using a divergent color band (RdBu_r) at 2%–98% quantiles, with invalid values ​​marked as light gray; a PNG image is output for quick viewing. The true-color composite image is linearly stretched from B4–B3–B2 at 2%–98% quantiles with an optional gamma (default 1.1) before being output as a PNG, providing a visually improved display.

[0193] By applying the technical solutions of the above embodiments, the atmospheric correction method for remote sensing images of nearshore marine scenes described in the above embodiments can be based on robust flare correction using SWIR proxy. It utilizes the physical characteristics of water reflection and path radiation of SWIR1(B6) under deep-water conditions to construct a scene-adaptive flare proxy, achieving flare subtraction. Furthermore, it uses multi-band joint luminosity index inversion, leveraging the spectral characteristics of the visible spectrum reflectance of clean deep water bodies to effectively smooth single-band random noise and improve the accuracy of AOT estimation results.

[0194] In some embodiments, as a specific implementation of the atmospheric correction method for remote sensing images of nearshore marine scenes described in the above embodiments, some embodiments of this application also provide an atmospheric correction system for remote sensing images of nearshore marine scenes, such as... Figure 5 As shown, the system includes:

[0195] The data acquisition module is used to acquire remote sensing data, which includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light band, near-infrared band, and short-wave infrared band. The auxiliary parameters include solar geometric parameters and sea surface state parameters.

[0196] A mask extraction module is used to extract target masks from the remote sensing data, the target masks including water masks and deep-water masks;

[0197] The flare correction module is used to perform solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum. The secondary simulation model is a relationship model between the visible light band and the flare intensity proxy based on the flare coefficient. The flare intensity proxy uses a first shortwave infrared band as a specular reflection proxy, constructed in conjunction with a reference anchor point. The reference anchor point is obtained by statistically analyzing a preset percentile of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained using robust regression fitting with the median slope method.

[0198] The index calculation module is used to calculate the multi-band luminosity index of the deep-water mask after solar flare correction. The multi-band luminosity index is calculated based on the pixel value of each pixel in the blue visible light, green visible light and red visible light bands of the deep-water mask.

[0199] The inversion module is used to invert the optimal aerosol optical thickness value based on the multi-band luminosity index using a radiative transfer model; the radiative transfer model is a physical parameter calculation model based on the auxiliary parameter settings; the physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance.

[0200] The result output module is used to calculate the atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value, and to generate correction result data based on the atmospheric correction coefficients. The correction result data includes corrected images, difference maps, and color composite images.

[0201] By applying the technical solutions of the above embodiments, the atmospheric correction system for remote sensing images of nearshore marine scenes described in the above embodiments can, after the data acquisition module acquires remote sensing data, the mask extraction module extracts the target mask, and the flare correction module then uses a secondary simulation model of satellite signals in the solar spectrum to perform solar flare correction. Then, the index calculation module calculates multi-band occultation indices so that the inversion module can invert the optimal aerosol optical thickness value through a radiative transfer model. The result output module then calculates atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value to generate correction result data. The system can perform robust flare correction based on a SWIR proxy, utilizing the water reflection and path radiation physical characteristics of SWIR1 (B6) under deep-water conditions to construct a scene-adaptive flare proxy to achieve flare subtraction. Furthermore, it uses multi-band joint occultation indices for inversion, utilizing the spectral characteristics of the reflectance of clean deep water bodies in the visible spectrum to effectively smooth single-band random noise and improve the accuracy of AOT estimation results.

[0202] It should be noted that other corresponding descriptions of the functional units involved in the remote sensing image atmospheric correction system for nearshore marine scenes provided in this application embodiment can be found in the corresponding descriptions in the remote sensing image atmospheric correction method for nearshore marine scenes provided in the above embodiments, and will not be repeated here.

[0203] This application also provides a computer device, specifically a personal computer, server, network device, etc. The computer device includes a bus, processor, memory, and communication interface, and may also include input / output interfaces and a display device. The processor of the computer device provides computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database of the computer device stores location information. The network interface of the computer device is used for communication with external terminals via a network connection. When the computer program is executed by the processor, it implements the steps in the various method embodiments.

[0204] Those skilled in the art will understand that the structure of the computer device described above is only a partial structure related to the solution of this application, and does not constitute a limitation on the computer device to which the solution of this application is applied. A specific computer device may include more or fewer components, or combine certain components, or have different component arrangements.

[0205] In one embodiment, a computer-readable storage medium is also provided, which may be non-volatile or volatile, and a computer program is stored thereon, which, when executed by a processor, implements the steps in the above method embodiments.

[0206] In one embodiment, a computer program product is also provided, including a computer program that, when executed by a processor, implements the steps in the above method embodiments.

[0207] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties.

[0208] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above methods.

[0209] Any references to memory, database, or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc.

[0210] Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can take many forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0211] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, distributed databases based on blockchain. The processors involved in the embodiments provided in this application may be, but are not limited to, general-purpose processors, graphics processors, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc.

[0212] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0213] The embodiments described above are merely examples of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of this patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these modifications and improvements all fall within the protection scope of this application.

Claims

1. An atmospheric correction method for remote sensing images of nearshore marine scenes, characterized in that, The method includes: Acquire remote sensing data, which includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light band, near-infrared band, and short-wave infrared band. The auxiliary parameters include solar geometric parameters and sea surface state parameters. Extracting target masks from the remote sensing data, the target masks including water masks and deep-water masks; extracting target masks from the remote sensing data includes: reading input bands from the remote sensing data, the input bands including a green visible light band, a red visible light band, a near-infrared band, a first shortwave infrared band, and a second shortwave infrared band; calculating spectral indices based on the reflectance values ​​of the input bands, the spectral indices including an improved normalized differential water index calculated based on the reflectance of the green visible light band and the first shortwave infrared band, and a normalized differential vegetation index calculated based on the reflectance of the near-infrared band and the red visible light band; setting water body screening conditions based on the spectral indices; extracting water body masks from the remote sensing data according to the water body screening conditions, the water body masks including loose water body masks and strict water body masks; the loose A water mask is used as a decoupling application domain; a strict water mask is used as a sampling domain; effective values ​​are extracted from the strict water mask, the effective values ​​including a first effective value extracted in the first shortwave infrared band and a second effective value extracted in the second shortwave infrared band; a preset quantile threshold is set according to the effective values, the preset quantile threshold including a first threshold and a second threshold; the preset quantile threshold is dynamically set according to the number of effective samples; a deep water mask is extracted from the strict water mask based on the preset quantile threshold, the deep water mask being a set of pixels whose reflectance in the first shortwave infrared band is less than the first threshold and whose reflectance in the second shortwave infrared band is less than the second threshold; statistical information is output, the statistical information including the number of pixels and the percentage of pixels in the loose water mask, the strict water mask, and the deep water mask; Solar flare correction is performed on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum. The secondary simulation model is a relationship model between the visible light band and the flare intensity proxy based on the flare coefficient. The flare intensity proxy uses a first shortwave infrared band as a specular reflection proxy, constructed in conjunction with reference anchor points. The reference anchor points are obtained by statistically analyzing a preset percentile of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained using robust regression fitting with the median slope method. The solar flare correction on the remote sensing image corresponding to the target mask using the secondary simulation model of satellite signals in the solar spectrum includes: based on the first shortwave infrared band... A specular reflection proxy is set for reflectivity; a preset percentile of the first shortwave infrared band is statistically analyzed within the deep-water mask to obtain a reference anchor point; a flare intensity proxy is constructed based on the specular reflection proxy and the reference anchor point, where the flare intensity proxy is the non-negative part of the difference between the specular reflection proxy and the reference anchor point; the reflectivity of the visible light band is extracted from the remote sensing data; based on the reflectivity of the visible light band and the flare intensity proxy, a robust regression estimation of the relationship slope is performed using the median slope method, where the relationship slope is numerically constrained within a preset numerical range; a first relationship model between the reflectivity of the visible light band and the flare intensity proxy is constructed based on the relationship slope. The multi-band luminosity index of the deep-water mask after solar flare correction is calculated. The multi-band luminosity index is calculated based on the pixel value of each pixel in the blue visible light, green visible light and red visible light bands of the deep-water mask. The optimal aerosol optical thickness value is retrieved by inverting the multi-band luminosity index using a radiative transfer model. The radiative transfer model is a physical parameter calculation model based on the auxiliary parameters. The physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance. The atmospheric correction coefficients for each band are calculated based on the optimal aerosol optical thickness value, and the correction result data is generated based on the atmospheric correction coefficients. The correction result data includes the corrected image, the difference map, and the color composite image.

2. The method according to claim 1, characterized in that, Acquiring remote sensing data, including: Acquire raw remote sensing data; Multiple band data are read from the original remote sensing data using a mask array method, wherein the band data retains information on data-free values ​​and mask information; Based on the key metadata of the source file of the original remote sensing data, the band data is cropped and resampled. Value range processing is performed on the band data after data cropping and resampling to generate the remote sensing data.

3. The method according to claim 1, characterized in that, Performing solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum also includes: Obtain the relative weights, which are either default weight values ​​or set using an offline simulation table; The fitted independent variables are calculated based on the relative weights and the flare intensity proxy. Using the reflectivity of the visible light band as the dependent variable, a stacked regression was performed on the deep-water mask using the median slope method to estimate the global scale parameters. Based on the global scale parameters, a second relationship model is constructed between the reflectivity of the visible light band and the flare intensity surrogate quantity.

4. The method according to claim 1, characterized in that, Calculate the multi-band luminosity index of the deep-water mask after solar flare correction, including: The dark image band is determined, which includes the blue visible light band, the green visible light band, and the red visible light band; Extract the pixel values ​​of the dark image band from the deep-water mask; For each pixel, the average value of the corresponding pixel value of the dark image band is calculated to obtain the multi-band darkness index; Based on the multi-band luminance index, the index mean and index standard deviation of the pixel set are calculated; Obtain a preset dark image coefficient, and calculate the theoretical dark pixel reflectance based on the dark image coefficient, the mean of the index, and the standard deviation of the index.

5. The method according to claim 4, characterized in that, Based on the aforementioned multi-band luminosity index, the optimal aerosol optical thickness value is retrieved through a radiative transfer model, including: The geometric constraints of the radiative transfer model are set according to the auxiliary parameters; The objective function of the radiative transfer model is constructed based on the theoretical dark pixel reflectance. The objective function is used to represent the difference between the actual observed dark pixel reflectance and the theoretical dark pixel reflectance. Set the search interval and step size for the radiative transfer model; According to the preset output criteria, an inversion is performed based on the radiative transfer model to determine the optimal aerosol optical thickness value; the preset output criteria are used to take the aerosol optical thickness value that minimizes the objective function as the optimal aerosol optical thickness value.

6. The method according to claim 5, characterized in that, Based on the aforementioned multi-band luminosity index, the optimal aerosol optical thickness value is inverted using a radiative transfer model, and the method further includes: The physical parameters are calculated band by band. The apparent reflectance is extracted from the remote sensing data, and the surface reflectance of each pixel is calculated based on the apparent reflectance and the physical parameters. The effective pixels are determined according to the surface reflectance. The number of effective pixels and the percentage of effective pixels are statistically analyzed.

7. An atmospheric correction system for remote sensing images of nearshore marine scenes, characterized in that, The system is applied to the method according to any one of claims 1-6; the system comprises: The data acquisition module is used to acquire remote sensing data, which includes remote sensing images and auxiliary parameters. The remote sensing images are reflectance images covering the visible light band, near-infrared band, and short-wave infrared band. The auxiliary parameters include solar geometric parameters and sea surface state parameters. A mask extraction module is used to extract target masks from the remote sensing data, the target masks including water masks and deep-water masks; The flare correction module is used to perform solar flare correction on the remote sensing image corresponding to the target mask using a secondary simulation model of satellite signals in the solar spectrum. The secondary simulation model is a relationship model between the visible light band and the flare intensity proxy based on the flare coefficient. The flare intensity proxy uses a first shortwave infrared band as a specular reflection proxy, constructed in conjunction with a reference anchor point. The reference anchor point is obtained by statistically analyzing a preset percentile of the first shortwave infrared band within the deep-water mask. The flare coefficient is obtained using robust regression fitting with the median slope method. The index calculation module is used to calculate the multi-band luminosity index of the deep-water mask after solar flare correction. The multi-band luminosity index is calculated based on the pixel value of each pixel in the blue visible light, green visible light and red visible light bands of the deep-water mask. The inversion module is used to invert the optimal aerosol optical thickness value based on the multi-band luminosity index using a radiative transfer model; the radiative transfer model is a physical parameter calculation model based on the auxiliary parameter settings; the physical parameters include atmospheric path reflectance, spherical albedo, total downlink transmittance, and total uplink transmittance. The result output module is used to calculate the atmospheric correction coefficients for each band based on the optimal aerosol optical thickness value, and to generate correction result data based on the atmospheric correction coefficients. The correction result data includes corrected images, difference maps, and color composite images.

Images