A universal on-orbit radiometric calibration method for ocean color satellites

Abstract

The present invention discloses a universal on-orbit radiation calibration method for water color satellites, belonging to the field of radiation calibration. The method comprises a data analysis and quality control module, a radiation calibration module, a calibration coefficient comprehensive analysis module and a working status recording module. The present invention optimizes the algorithm of the hyperspectral atmospheric absorption coefficient and realizes the compatibility of the on-orbit radiation calibration algorithm with multiple types of water color satellite payloads.

Description

Technical Field

[0001] The present invention relates to the field of radiation calibration, and in particular to a universal on-orbit radiation calibration method for a water color satellite. Background Art

[0002] The extensive research and application of ocean color remote sensing (OCRS) products (such as chlorophyll a concentration, suspended particulate matter, and colored dissolved organic matter) have brought significant scientific and economic benefits to our deeper understanding of the ocean. To obtain high-quality products of these optically active components, the uncertainty of ocean radiometric measurements (such as water-leaving radiance (Lw) or remotely sensed reflectance (Rrs)) must be controlled to within 5% through atmospheric correction (Hooker et al., 1992). However, because Lw in the visible (VIS) band accounts for only 10% or less of the total radiance (Lt) measured by ocean color sensors in the top-of-atmosphere (TOA) region, the uncertainty of Lt must be less than 0.5%. To meet this requirement, ocean color sensors undergo a series of radiometric calibration tasks to couple and reduce the uncertainties introduced by atmospheric correction and ocean color retrieval algorithms. While existing publicly available OCRS on-orbit radiometric calibration systems already include the aforementioned calibration methods, these systems are typically sensor-specific and difficult to directly apply to other similar payloads. Furthermore, they are unable to meet the requirements of future high-spectral and high-spatial-resolution sensors. Given the current operational status of numerous OCRS systems internationally and the upcoming launch of more advanced ocean color sensors, this research aims to develop a universal on-orbit radiometric calibration system, termed the Universal Ocean Color Radiometric Calibration Software. Summary of the Invention

[0003] The purpose of the present invention is to solve the defects in the prior art and to propose a universal water color satellite on-orbit radiation calibration method.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] A universal on-orbit radiation calibration method for water color satellites, comprising a data analysis and quality control module, a radiation calibration module, a calibration coefficient comprehensive analysis module, and a working status recording module;

[0006] The specific working steps of the universal ocean color satellite on-orbit radiometric calibration method based on the above module are as follows:

[0007] Step 1: Input the data into the data analysis and quality control module to analyze the input calibration load data and various auxiliary data;

[0008] Step 2: Select the corresponding calibration algorithm based on the input data;

[0009] Step 3: When it is cross calibration, the radiation calibration module is scheduled to dynamically select the calibration area;

[0010] Step 4: Perform quality control screening to extract pixels that are suitable for the calibration method and meet the matching criteria, and generate the corresponding matching data set;

[0011] Step 5: The data analysis and quality control module outputs the scheduling instructions and matching data sets, and the radiation calibration module calculates the calibration coefficients according to different calibration algorithms;

[0012] Step 6: Summarize the calibration coefficients and generate corresponding daily or periodic reports and calibration coefficient adoption recommendations through the calibration coefficient comprehensive analysis module. Then, the working status recording module records system working information, tracks the working status of other modules, and generates a work log.

[0013] As a further solution of the present invention, the quality control principles of the data analysis and quality control module specifically include:

[0014] For field observation data, the QA scoring system is used for scoring. If the score generated by the QA scoring system is greater than 0.8, the data is directly adopted; otherwise, it is converted into remote sensing reflectance;

[0015] For satellite image data, when L2 level data is provided, it is screened according to the provided internationally accepted 32-bit flag to eliminate the influence of land, solar flares, saturated radiance, and cloud factors; otherwise, the satellite impact data is screened based on the principles of wind speed less than 6m / s, 5×5 grid top of atmosphere reflectivity variation coefficient within 0.15, solar zenith angle less than 70 degrees, satellite observation angle less than 56 degrees, and the longest near-infrared band top of atmosphere remote sensing reflectivity less than 0.027.

[0016] As a further solution of the present invention, the specific steps of the quality control screening in the fourth step are as follows:

[0017] S1.1: If the calibration is for system replacement, exclude data with a time difference greater than 1 hour and pixels in satellite images whose distance from the on-site observation site exceeds the spatial resolution of the calibration payload;

[0018] S1.2: If cross-calibration is used, data with a time difference greater than 30 minutes will be excluded, as will data where the spatial distance between the nearest pixels of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel, and the overlap range of the two payload images is less than 15×15 pixels.

[0019] As a further solution of the present invention, the specific steps of calculating the calibration coefficients by the radiation calibration module according to different calibration algorithms in the fifth step are as follows:

[0020] S2.1: Hyperspectral Rayleigh LUTs are established based on OSOAA. By simulating different solar zenith angles, satellite zenith angles, and wind speeds, the three Stokes vector coefficients of Rayleigh reflectance are obtained. The relative azimuth angle is then Fourier expanded based on each set of Stokes vector coefficients.

[0021] S2.2: Based on the 80 aerosol models classified by Ahmad, the single scattering albedo and scattering phase function of the aerosol are calculated using the Mie scattering model. The quadratic coefficient fitting values ​​of the aerosol single and multiple scattering albedo are obtained by simulation using the OSOAA radiative transfer model at different aerosol optical depths.

[0022] S2.3: Using the spectral response correction method, the spectral correction factors generated by the spectral response functions of the two satellites based on the long-term hyperspectral MOBY water-leaving radiance dataset stored in the system are used to construct the spectral correction slope factors and intercept factors between the bands. After the initial generation, these factors are stored. The water-leaving radiance provided by the reference payload is then converted to the water-leaving radiance at the central wavelength of the calibration payload.

[0023] S2.4: Obtain gas absorption coefficients from the Hitran database. Calculate the transmittance of each gas by accumulating the beam spread function based on the absorption coefficients and combining them with the gas density. Then, use an aerosol model containing aerosol lookup tables (LUTs) to build an aerosol model at different zenith angles, covering a dynamic range of 0° to 80°, to obtain the diffuse transmittance of each gas. The total atmospheric absorption transmittance is then calculated as the product of the absorption coefficients of each gas.

[0024] S2.5: For the coupled ocean-atmosphere system, in the absence of solar flare and sea whitecaps, the total radiance is calculated by combining it with atmospheric reflected radiation using a radiative transfer model.

[0025] As a further solution of the present invention, the specific calculation formula for the total radiant brightness described in S2.5 is as follows:

[0026] Where, represents the total radiance; Represents the Rayleigh scattered radiation produced by scattering of air molecules; represents aerosol scattered radiation; represents the diffuse transmittance from the water surface to the sensor; Represents the total atmospheric absorption and transmittance.

[0027] As a further solution of the present invention, the specific calculation formula for the water-leaving radiance at the central wavelength of the calibration load described in S2.3 is as follows

[0028]

[0029] Where, Represents the water-leaving radiance of the reference source in the corresponding band; represents the reference source spectral response function; Represents the spectral response function of the load to be calibrated;

[0030] The specific calculation formula for the total atmospheric absorption and transmittance mentioned in S2.4 is as follows:

[0031]

[0032]

[0033]

[0034]

[0035]

[0036]

[0037]

[0038] Where, From the absorption coefficient Calculated in represents the gas density; The representative interval can be selected according to the bandwidth of the calibration payload band; and Respectively Gas permeability from the Sun to the ocean surface and from the ocean surface to the satellite; and represent the solar zenith angle and satellite observation angle respectively; and They represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite, respectively.

[0039] As a further solution of the present invention, the calibration coefficient comprehensive analysis module is divided into two modes: daily and periodic. In the daily mode, the calibration report displays the process quantity of each calibration algorithm on that day. In the periodic mode, different calibration coefficients obtained by different calibration methods at different times are evaluated to achieve a comprehensive analysis of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within a specified time, and provide a set of calibration coefficients finally adopted.

[0040] As a further solution of the present invention, the specific calculation steps of the calibration coefficient finally adopted are as follows:

[0041] S3.1: Compare the relative percentage differences between the calibration coefficients obtained by different calibration methods or the calibration coefficients obtained by the same method within a specified time period. If the relative percentage difference is less than 5%, arrange all calibration gain coefficients in ascending order and calculate the average of the semi-quartile range as the final calibration coefficient;

[0042] S3.2: If the relative percentage difference is greater than 5%, then based on the stability of the trend over time, remove the abnormal data points that deviate from the mean of the calibration gain coefficient by one standard deviation, and recalculate the relative percentage difference to obtain the final calibration coefficient.

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

[0044] The present invention replaces the single lookup table in various calibration algorithms with a 1nm interval hyperspectral Rayleigh and aerosol lookup table established based on the OSOAA radiation transfer model and the long-term observation data of multiple water body stations in AERONET-OC by Ahmad et al., optimizes the algorithm of hyperspectral atmospheric absorption coefficient, and constructs spectral correction factors for the spectral response differences of the water-leaving radiance of each sensor for the cross-calibration algorithm, thus realizing the compatibility of the on-orbit radiation calibration algorithm for multiple types of water color satellite payloads. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0046] Figure 1 This is a system block diagram of a universal on-orbit radiometric calibration method for water color satellites proposed in the present invention;

[0047] Figure 2 This is an operational flow chart of a universal on-orbit radiometric calibration method for water color satellites proposed in the present invention.

[0048] Figure 3 Spectral variation diagram of three parameters of Rayleigh reflectance (wind speed is 2m / s) under different observation geometries (triangles represent the convolution results of the hyperspectral Rayleigh lookup table; circles represent the MODIS lookup table);

[0049] Figure 4 Spectral changes of extinction coefficient and scattering phase function (scattering angle is 2°) at different relative humidity (triangles represent the convolution results of the hyperspectral Rayleigh lookup table; circles represent the MODIS lookup table);

[0050] Figure 5 This is a comparison of the water-offset radiance scatter plot of the Sentinel-3A OLCI spectral response after correction and the HY3A COCTS2 central band;

[0051] Figure 6 This is the timing diagram of the calibration coefficients on the COCTS-HY1C / D satellite;

[0052] Figure 7 This is the COCTS-HY1C / D cross-calibration coefficient spectrum. DETAILED DESCRIPTION

[0053] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0054] Reference Figure 1-2 A universal on-orbit radiation calibration method for water color satellites includes a data analysis and quality control module, a radiation calibration module, a calibration coefficient comprehensive analysis module and a working status recording module.

[0055] Among them, the data analysis and quality control module reads and organizes according to the dedicated analysis unit of each data source, and at the same time determines the calibration method based on the input data source, and then performs quality control on the data in a unified format. The quality control criteria adopt corresponding data matching criteria for different data sources.

[0056] For in-situ observation data (MOBY), data with the highest quality rating were directly adopted. Otherwise, they were converted to remote sensing reflectance and scored using the QA scoring system (Wei, J., Lee, Z.-P., and Shang, S. (2016). A system to measure the data quality of spectral remote sensing reflectance of aquatic environments. J. Geophys. R., 121, 8189–8207.). Data with a score greater than 0.8 were adopted. For satellite imagery data, when L2-level data were provided, they were screened based on the internationally recognized 32-bit flags provided to exclude factors such as land, solar flares, saturated radiance, and clouds. Otherwise, they were screened based on the following criteria: 1) wind speed less than 6 m / s; 2) coefficient of variation of top-of-atmosphere reflectance within a 5×5 grid within 0.15; 3) solar zenith angle less than 70°; 4) satellite observation angle less than 56°; and 5) top-of-atmosphere remote sensing reflectance less than 0.027 in the longest near-infrared band.

[0057] For pixels that meet the quality control conditions, the corresponding spatiotemporal matching criteria are selected based on the calibration algorithm.

[0058] System alternative calibration: 1) Eliminate data with a time difference greater than 1 hour; 2) Eliminate pixels in satellite images whose distance from the on-site observation site exceeds the spatial resolution of the calibration payload.

[0059] Cross-calibration: 1) Eliminate data with a time difference greater than 30 minutes; 2) Eliminate data when the overlap range of the two payload images is less than 15×15 pixels; 3) Eliminate pixels when the spatial distance between the nearest pixels of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel.

[0060] The radiometric calibration module includes four algorithms: on-board calibration, system replacement calibration, cross calibration, and land calibration. The principles of each algorithm are introduced below:

[0061] a) On-board calibration

[0062] The onboard calibration method of the system is solar calibration. The main calibration devices of the onboard calibration system include the solar diffuse reflector (SD), the solar diffuse reflector stability detector (SDSM), the solar attenuation screen, the rotating telescope (RTA) and the half-angle mirror (HAM).

[0063] The sunlight is attenuated by the solar attenuation screen and enters the solar diffuse reflection plate. It then passes through the rotating telescope and is reflected by the half-angle mirror into the calibration sensor. At this time, the solar radiance can be expressed as:

[0064]

[0065] in: Indicates the use of a half-angle mirror Solar radiance at the pupil when facing; is the extraterrestrial solar irradiance; Expressed as the solar-terrestrial correction factor; and are the solar zenith angle and azimuth angle in the SD coordinate system respectively;

[0066] and are the observation zenith angle and azimuth angle of the calibration payload observation SD, respectively; is the bidirectional reflectance distribution function of the solar diffuse reflector;

[0067] and are the solar zenith angle and azimuth angle of the sun irradiating the solar attenuation screen, respectively; is the transmittance of the solar attenuation screen; BRDF correction factor for the solar diffuse reflection.

[0068] The above formula (1) can calculate the radiance value entering the calibration load. On this basis, the solar calibration gain coefficient in each band can be calculated by the following formula: .

[0069]

[0070] in: To use a half-angle mirror DN value when observing the solar diffuse reflection plate, To use a half-angle mirror Observe the DN value of cold air at surface time.

[0071] b) System alternative calibration

[0072] The system replacement calibration and cross calibration gain coefficient is defined as the ratio of the theoretical top-of-atmosphere radiance to the actual observed top-of-atmosphere radiance, as shown in formula (3).

[0073]

[0074] The actual observed top-of-atmosphere radiance is obtained by combining the DN value observed by the calibration payload data and the absolute radiation calibration coefficient calculated in the laboratory before launch with formula (3), and then converted into the top-of-atmosphere reflectivity using formula (4). The theoretical observed top-of-atmosphere radiance is derived to the top of the atmosphere through the system replacement calibration algorithm combined with field-measured water body information and meteorological auxiliary data.

[0075]

[0076] in: is the solar irradiance; Corrected distance for the Sun and Earth; is the solar zenith angle.

[0077] The system substitution calibration algorithm is often described as the inverse of the atmospheric correction algorithm. Atmospheric correction involves removing the atmospheric path reflectance from the top-of-atmosphere reflectance received from water color satellite sensors to obtain the water signal. However, the atmospheric correction algorithm itself is affected by multiple factors, including the optical properties of the water, the atmosphere, and the Sun-satellite observation geometry, which can introduce errors and increase the uncertainty of the inverted water-offset radiance beyond the specified range. By incorporating the atmospheric correction algorithm into the system substitution calibration algorithm, the uncertainty of atmospheric correction is coupled into the calibration process, minimizing the deviation between the remote sensing reflectance obtained from the top-of-atmosphere reflectance inversion and the measured data. The cross-calibration algorithm follows the same process as the system substitution calibration algorithm.

[0078] For the ocean-atmosphere coupled system, in the absence of solar flares, the top-of-atmosphere reflectivity The equation can be approximately written as:

[0079]

[0080] in: and represent Rayleigh scattering and aerosol (including the interaction between molecules and aerosol) reflectivity respectively; represents the white-hat reflectivity; represents remote sensing reflectance; and represent the atmospheric diffuse transmittance from the sun to the water and from the water to the sensor, respectively; represents the atmospheric transmittance.

[0081] Formula (5) shows that the top-of-atmosphere reflectivity includes contributions from Rayleigh scattering, aerosol scattering, water remote sensing reflectivity, and sea surface white caps. The contribution from solar flares is ignored in the above calculation formula. To avoid errors, pixels greater than 0 are removed by normalizing the solar flare coefficient. Each parameter is described below.

[0082] First of all, the white cap on the sea surface is caused by the white bubbles formed by the breaking of waves under the action of wind. It affects the reflectivity of the top of the atmosphere by reflecting sunlight and sky light. It is mainly manifested in the relative relationship with the wind speed on the sea surface. The white cap reflectivity of water bodies with a speed less than 6.33 m / s is 0, and the specific calculation formula is shown in (6). Normalized white hat reflectance, is the sea surface wind speed.

[0083]

[0084] In order to simplify the algorithm calculation process and improve the calculation efficiency, the system alternative calibration and cross calibration are all performed by lookup table. The Rayleigh lookup table is divided into multiple files by band. It contains parameters such as Rayleigh optical thickness, solar zenith angle, satellite observation angle, depolarization factor, wind speed and Rayleigh reflectivity coefficient. It can be calculated based on the wind speed and air pressure provided by the auxiliary data and the observation geometry provided by the calibration satellite data. .

[0085] Atmospheric transmittance calculation requires atmospheric pressure, Rayleigh optical thickness and observation geometry, where atmospheric pressure and observation geometry are provided by auxiliary data and uncalibrated data respectively, and Rayleigh optical thickness

[0086] The calculation is shown in formula (7), the atmospheric single gas transmittance It is expressed as the product of atmospheric gas content and gas absorption coefficient, which is used to correct the influence of observation geometry on atmospheric transmittance. The specific calculation formula is shown in (8).

[0087]

[0088]

[0089] in: are the current atmospheric pressure and the standard atmospheric pressure respectively; is the satellite observation angle; Indicates the type of gas.

[0090] The system's alternative calibration process separates the calculations for the NIR-SWIR and UV-VIS bands, performing a two-step calculation. The first step is to implement NIR-SWIR calibration based on the "dark pixel" hypothesis and known aerosol types in the water. The second step, using an atmospheric correction algorithm, extrapolates the visible band calibration gain coefficient from the two NIR / SWIR bands. The following details the calculation of aerosol scattering contributions based on this process.

[0091] Based on the two prerequisites of a system replacement calibration algorithm, this paper selected the South Pacific Gyre and the South Indian Ocean Gyre, both of which possess relatively clean, uniform, and stable non-nearshore marine atmospheric conditions, for near-infrared and shortwave infrared (NIR) calibration. The aerosol types in these two gyres are assumed to be a 50 / 50 mixed aerosol model consisting of a 70% relative humidity marine aerosol (M70) and a 90% relative humidity marine aerosol model (M90). First, in these waters with known aerosol types, an initial system calibration gain of 1 was assumed for the longer shortwave infrared band (1640 nm). The aerosol reflectance in the remaining shortwave infrared band (1245 nm) was estimated using an aerosol lookup table. Subsequently, aerosol optical properties in the NIR band were estimated using a "dark pixel" atmospheric correction algorithm based on the NIR band. Finally, the NIR-SWR calibration gain was derived based on the calculations of the other components mentioned above. At this point, the 1640nm band system replacement calibration gain coefficient is 1. Since the system replacement calibration algorithm assumes that the longer near-infrared band (865nm) is 1, the 1640nm band calibration gain coefficient is adjusted based on the 865nm band system replacement calibration gain coefficient. This process is repeated multiple times until the 865nm band system replacement calibration gain coefficient is 1. The final obtained near-infrared-short-wave infrared system replacement calibration gain coefficient is the final value.

[0092] The near-infrared-band infrared system is then used to replace the calibration gain coefficient for the water area of ​​the on-site observation site. At this time, based on the two near-infrared / short-wave infrared band calibration gain coefficients, the aerosol type and aerosol multiple scattering reflectivity of the area are determined, and then the near-infrared / short-wave infrared aerosol single scattering reflectivity is obtained by combining the lookup table calculation. Finally, according to the conversion formula of aerosol single scattering reflectivity and aerosol optical thickness, the extinction coefficient of different bands And the UV-visible aerosol scattering reflectivity can be calculated by aerosol single and multiple conversions.

[0093]

[0094]

[0095] in: is the aerosol optical depth; is the scattering phase function;

[0096] is the aerosol single scattering albedo.

[0097] The remote sensing reflectivity data is provided by the on-site measurement site and the calculation formula is shown below.

[0098]

[0099] in: Provide remote sensing reflectance converted from off-water radiance for ocean optical buoys, is the calibration load spectral response function. is the remote sensing reflectivity corresponding to the central wavelength of the calibration payload.

[0100] c) Cross-calibration

[0101] The cross-calibration algorithm differs from the system replacement calibration algorithm in that it calculates the aerosol scattering reflectance using a reference payload to provide the aerosol optical depth and aerosol type for the reference band in the cross-region. The estimation of the aerosol reflectance in the UV-Vis band is based on the estimation of the aerosol scattering reflectance in the near-infrared band based on the "dark pixel" assumption. The process is as follows:

[0102] The aerosol single scattering reflectance is estimated using an aerosol lookup table and top-of-atmosphere reflectance in two near-infrared bands observed by a reference payload. Because aerosol single scattering reflectance is related to inherent aerosol properties such as the scattering phase function, single scattering albedo, and the aerosol attenuation coefficient, the aerosol single scattering reflectance is calculated for each aerosol model. This allows the most appropriate aerosol model to be selected, ensuring the best match between the calculated radiance and the measured value. Furthermore, based on the direct determination of aerosol type through statistical analysis, the two known aerosol models are directly extrapolated to estimate aerosol reflectance in the UV-visible band.

[0103] Remote sensing reflectance is directly provided by the reference payload and is obtained by performing an atmospheric correction procedure on the reference payload data. This involves differences in the central wavelength deviation and spectral response between the reference band and the payload to be calibrated. This can be achieved by employing a band shift algorithm based on a bio-optical model, selecting the nearest bands between the calibration payload and the reference payload band for conversion.

[0104] d) Land calibration

[0105] The 6S radiation transfer model is used for simulation in land calibration. It is assumed that under cloudless atmospheric conditions, the top of the atmosphere radiance is It can be expressed as:

[0106]

[0107]

[0108] in, is the relative azimuth angle between the sun and the satellite, is the satellite zenith angle, is the solar zenith angle, is the spectral reflectance of the ground target, is the albedo of the large balloon surface, The transmittance caused by the absorption of atmospheric molecules.

[0109] are the diffuse transmittances on the two paths from the sun and from the sensor to the ground target, and are the top-of-atmosphere radiance and top-of-atmosphere reflectivity, is the cosine of the solar zenith angle, is the solar irradiance at the top of the atmosphere,

[0110] is the distance between the Earth and the Sun (unit: AU).

[0111] The 6S radiation transfer model can be used to calculate the top-of-atmosphere radiance in the visible and ultraviolet bands. Land substitution calibration factor It involves the differences between different detectors. The calculation formula is:

[0112]

[0113] Based on the above calculation principle, the algorithm is improved. Due to the differences in the spectral response functions of various water color satellite sensors, in order to make the algorithm adaptive to multiple types of water color satellite sensors during the calibration process, the system uses the OSOAA (ocean successive orders with atmosphere-advanced) radiation transfer model to establish a 1nm interval hyperspectral Rayleigh lookup table and an aerosol lookup table. The calculation formulas for each parameter in the lookup table are as follows:

[0114]

[0115] in: and They are the lookup table parameters for the central band of the hyperspectral and calibration loads respectively; is the spectral response function.

[0116] Atmospheric transmittance is calculated line by line using the line-by-line (LBL) method, which uses the Hitran database to provide gas absorption lines. The K distribution method is introduced to improve LBL calculation efficiency. First, the gas absorption coefficient is obtained from the Hitran database. The transmittance of each gas is calculated using the cumulative beam spread function of the absorption coefficient and the gas density, using the following formula.

[0117]

[0118]

[0119] From the absorption coefficient Calculated in. is the gas density (provided by auxiliary data). The selection can be made based on the bandwidth of the calibration load band. The above formula can be used to calculate the absorption transmittance of different gases, and combined with the observation geometry, the total gas transmittance can be obtained:

[0120]

[0121]

[0122]

[0123]

[0124]

[0125] in: and Respectively The gas permeability from the sun to the ocean surface and from the ocean surface to the satellite; and Indicates the solar zenith angle and satellite observation angle respectively; and They represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite, respectively.

[0126] In the cross-calibration algorithm, due to the difference in spectral response between the calibration payload and the reference payload, and considering the uncertainty and computational efficiency caused by the uncorrected spectral function difference in the band shift algorithm based on the bio-optical model, this paper is based on the spectral response correction method proposed by Doelling et al. (Doelling, DR, et al. (2011). Spectral Reflectance Corrections for Satellite Intercalibrations UsingSCIAMACHY Data. IEEE Geoscience and Remote Sensing Letters, 9(1), 119-123.), and uses the long-time series hyperspectral data obtained by field measurements of the ocean optical buoy MOBY to construct the spectral correction slope factor between bands. and the intercept factor , thereby converting the water-free radiance provided by the reference payload into the water-free radiance at the central wavelength of the calibration payload. For the two satellite payloads that first establish the cross-calibration algorithm in the system, this algorithm needs to rely on the spectral response functions of the two satellites to generate spectral correction factors based on the long-term hyperspectral MOBY water-free radiance dataset stored in the system. After the first generation, it will be stored. Formula (23) shows the spectral correction slope factor and the intercept factor How to use it.

[0127]

[0128] in is the water-leaving radiance of the reference payload in the corresponding band, Consider the load spectrum response function, is the spectral response function of the load to be calibrated.

[0129] The specific working steps of a universal on-orbit radiometric calibration method for ocean color satellites based on the above modules are as follows:

[0130] Step 1: Input the data into the data analysis and quality control module to analyze the input calibration load data and various auxiliary data.

[0131] It should be further explained that the quality control principles of the data analysis and quality control module include:

[0132] For field observation data, the QA scoring system is used for scoring. If the score generated by the QA scoring system is greater than 0.8, the data is directly adopted; otherwise, it is converted into remote sensing reflectance;

[0133] For satellite image data, when L2 level data is provided, it is screened according to the provided internationally accepted 32-bit flag to eliminate the influence of land, solar flares, saturated radiance, and cloud factors; otherwise, the satellite impact data is screened based on the principles of wind speed less than 6m / s, 5×5 grid top of atmosphere reflectivity variation coefficient within 0.15, solar zenith angle less than 70 degrees, satellite observation angle less than 56 degrees, and the longest near-infrared band top of atmosphere remote sensing reflectivity less than 0.027.

[0134] Step 2: Select the corresponding calibration algorithm based on the input data.

[0135] Step 3: When it is cross calibration, the radiation calibration module is scheduled to dynamically select the calibration area.

[0136] Step 4: Perform quality control screening, extract pixels that are suitable for the calibration method and meet the matching criteria, and generate the corresponding matching data set.

[0137] Specifically, if it is a system replacement calibration, data with a time difference greater than 1 hour and pixels in the satellite image whose distance from the on-site observation site exceeds the spatial resolution of the calibration payload are excluded. If it is a cross calibration, data with a time difference greater than 30 minutes and the spatial distance between the nearest pixels of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel are excluded, and the overlap range of the two payload images is less than 15×15 pixels are excluded.

[0138] Step 5: The data analysis and quality control module outputs the scheduling instructions and matching data sets, and calculates the calibration coefficients according to different calibration algorithms through the radiation calibration module.

[0139] Specifically, a hyperspectral Rayleigh LUTs was established based on OSOAA. By simulating different solar zenith angles, satellite zenith angles and wind speeds, the three coefficients of the Stokes vector of the Rayleigh reflectance were obtained, and the relative azimuth angle was Fourier expanded according to the obtained Stokes vector coefficients. Based on the 80 groups of aerosol models divided by Ahmad, the single scattering albedo and scattering phase function of the aerosol were obtained by Mie scattering model calculation. The OSOAA radiation transfer model was used to simulate at different aerosol optical depths to obtain the quadratic coefficient fitting values ​​of the single and multiple scattering albedo of the aerosol. The spectral response correction method was used to construct the spectral correction factor generated by the spectral response function of the two satellites based on the hyperspectral MOBY water-leaving radiance dataset stored in the system over a long period of time. The spectral correction slope factor and intercept factor between bands are established and stored after the first generation. The water-free radiance provided by the reference payload is then converted to the water-free radiance at the central wavelength of the calibration payload. The gas absorption coefficient is obtained from the Hitran database, and the beam distribution function is accumulated through the absorption coefficient. The transmittance of each gas is calculated in combination with the gas density. An aerosol model containing aerosol LUTs is then used to establish it under different zenith angle conditions, covering a dynamic range of 0° to 80°, to obtain the diffuse transmittance of each gas. The total atmospheric absorption transmittance is then obtained based on the product of the absorption coefficients of each gas. For the coupled ocean-atmosphere system, in the absence of solar flash and whitecaps on the sea surface, the radiation transfer model is used to combine it with the atmospheric reflected radiation to simulate and calculate the total radiation brightness.

[0140] In this embodiment, the specific calculation formula for the total radiant brightness is as follows:

[0141]

[0142] Where, represents the total radiance; Represents the Rayleigh scattered radiation produced by scattering of air molecules; represents aerosol scattered radiation; represents the diffuse transmittance from the water surface to the sensor; Represents the total atmospheric absorption and transmittance.

[0143] The specific calculation formula for the water-off radiance at the central wavelength of the calibration load described in S2.3 is as follows:

[0144]

[0145] Where, Represents the water-leaving radiance of the reference source in the corresponding band; represents the reference source spectral response function; Represents the spectral response function of the load to be calibrated;

[0146] The specific calculation formula for the total atmospheric absorption and transmittance is as follows:

[0147]

[0148]

[0149]

[0150]

[0151]

[0152]

[0153]

[0154] Where, From the absorption coefficient Calculated in represents the gas density;

[0155] The representative interval can be selected according to the bandwidth of the calibration payload band; and Respectively Gas permeability from the Sun to the ocean surface and from the ocean surface to the satellite; and represent the solar zenith angle and satellite observation angle respectively;

[0156] and They represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite, respectively.

[0157] In addition, it should be noted that the configuration parameters of the hyperspectral Rayleigh LUTs are shown in Table 1 below:

[0158]

[0159] The dynamic parameters of the aerosol lookup table of the aerosol model are shown in Table 2 below:

[0160]

[0161] Step 6: Summarize the calibration coefficients and generate corresponding daily or periodic reports and calibration coefficient adoption recommendations through the calibration coefficient comprehensive analysis module. Then, the working status recording module records system working information, tracks the working status of other modules, and generates a work log.

[0162] It should be further explained that the calibration coefficient comprehensive analysis module is divided into two modes: daily and periodic. In the daily mode, the calibration report displays the process quantity of each calibration algorithm on that day. In the periodic mode, different calibration coefficients obtained by different calibration methods at different times are evaluated to achieve a comprehensive analysis of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within a specified time, and provide a set of calibration coefficients that are finally adopted.

[0163] In addition, it should be noted that the calibration coefficient finally adopted is specifically obtained by comparing the relative percentage differences of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within a specified time. If the relative percentage difference is less than 5%, all calibration gain coefficients are arranged in ascending order, and the average value of the semi-quartile range is calculated as the final calibration coefficient. If the relative percentage difference is greater than 5%, abnormal data points with a deviation of one standard deviation from the mean of the calibration gain coefficient are eliminated based on the stability of the trend over time, and the relative percentage difference is calculated again to obtain the final calibration coefficient.

[0164] To verify the reliability of the hyperspectral aerosol and Rayleigh lookup tables established above, we verify them with the MODIS dedicated lookup table in NASA's operational system. First, we show the three parameters of Rayleigh reflectance. The following comparison shows the consistency of the corresponding parameters with MODIS. The spectral curves are compared under different conditions of wind speed of 2 m / s, solar zenith angle (solz), and satellite zenith angle (senz). It can be seen that the Rayleigh reflectance spectral similarity is 1 under the same observation angle, and the mean absolute relative error (MARD) of each parameter is better than 1%, as shown in Table 3.

[0165] Table 3 Different observation geometries Spectral curve statistical parameters (wind speed 2m / s)

[0166]

[0167] At the same time, in order to verify the feasibility of aerosol optical depth, the distribution of the extinction coefficient and scattering phase function (scattering angle is 2°) spectral curves of eight relative humidity models at the median fine mode fraction are shown here. Figure 3 It can be seen that the spectral shape of the extinction coefficient and scattering phase function obtained by hyperspectral convolution is consistent with that of the MODIS lookup table. The spectral similarity of the extinction coefficient and scattering phase function at various relative humidity levels is close to 1, and the MARD of each parameter is better than 0.5%, indicating a high correlation.

[0168] Table 4 Statistical parameters of the extinction coefficient and scattering phase function (scattering angle is 2°) spectral curve under relative humidity conditions

[0169]

[0170] To validate the algorithm for establishing spectral correction factors, we used the HY3A COCTS2 and Sentinel-3A OLCI as examples. Given the relatively low water-leaving radiance values ​​at wavelengths above 620 nm, only the 360-565 nm water-leaving radiance comparison scatter plot is presented here. Although the Sentinel-3A OLCI's minimum wavelength is 400 nm, the water-leaving radiance at 400 nm, after correction by the correction factor, still represents the water-leaving radiance of the HY3A COCTS2 in the UV band well. The correlation coefficient is greater than 0.995, and the correlation for other visible bands is better than 0.999. The MAPD is better than 1% when the center band difference is less than 20 nm, and good consistency is also shown when the center band deviation is 40 nm, thus verifying the effectiveness of the algorithm.

[0171] Table 5 Statistical parameters of water-offset radiance of the central band of Sentinel-3A OLCI after spectral response correction and HY3A COCTS2

[0172]

[0173] Subsequently, the onboard calibration and cross-calibration results of COCTS-HY1C / D will be realized based on this system, and the test data time range is from August 1, 2022 to August 1, 2023. First, the calibration coefficients are obtained by cross-calibrating the onboard calibration spectrometer and COCTS.

[0174] Depend on Figure 6 The onboard calibration coefficients of COCTS-HY1C at 490, 520, and 670 nm during this period were 0.99, 0.98, and 0.94, with standard deviations of 0.009, 0.010, and 0.019, and coefficients of variation of 0.88%, 1.06%, and 2.00%, respectively. The onboard calibration coefficients of COCTS-HY1D at 490, 520, and 670 nm were 0.97, 0.97, and 0.95, with coefficients of variation of 1.94%, 1.95%, and 2.44%, respectively.

[0175] In order to verify the stability of the calibration coefficients, considering that the sensor attenuates over time, the above-mentioned on-board calibration coefficients are used to establish stable radiation performance for COCTS-HY1C / D. On this basis, cross-calibration is performed. The coefficient of variation of the obtained calibration coefficients in each band within one year is basically better than 0.5%. The results are as follows: Figure 7 and as shown in Table 6. This verifies the applicability of the Grace-OC system in the application of the two sensors COCTS-HY1C / D, and also verifies the effectiveness and stability of the on-board calibration and cross-calibration algorithms;

[0176] Table 6 COCTS-HY1C / D cross-calibration statistical parameters

[0177]

[0178] The above is a detailed description of an embodiment of the present invention, but the content is only a preferred embodiment of the present invention and should not be considered to limit the scope of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the scope of the patent coverage of the present invention.

Claims

1. A universal ocean color satellite on-orbit radiometric calibration method, characterized in that: It includes data analysis and quality control module, radiation calibration module, calibration coefficient comprehensive analysis module and working status recording module; The specific working steps of the universal ocean color satellite on-orbit radiometric calibration method based on the above module are as follows: Step 1: Input the data into the data analysis and quality control module to analyze the input calibration load data and various auxiliary data; Step 2: Select the corresponding calibration algorithm based on the input data; Step 3: When it is cross calibration, the radiation calibration module is scheduled to dynamically select the calibration area; Step 4: Perform quality control screening to extract pixels that are suitable for the calibration method and meet the matching criteria, and generate the corresponding matching data set; Step 5: The data analysis and quality control module outputs the scheduling instructions and matching data sets, and the radiation calibration module calculates the calibration coefficients according to different calibration algorithms; Step 6: Summarize the calibration coefficients and generate corresponding daily or periodic reports and calibration coefficient adoption recommendations through the calibration coefficient comprehensive analysis module. Then, the working status recording module records system working information, tracks the working status of other modules, and generates a work log. The specific steps for calculating the calibration coefficients using the radiation calibration module according to different calibration algorithms in step 5 are as follows: S2.1: Hyperspectral Rayleigh LUTs are established based on OSOAA. By simulating different solar zenith angles, satellite zenith angles, and wind speeds, the three Stokes vector coefficients of Rayleigh reflectance are obtained. The relative azimuth angle is then Fourier expanded based on each set of Stokes vector coefficients. S2.2: Based on the 80 aerosol models classified by Ahmad, the single scattering albedo and scattering phase function of the aerosol are calculated using the Mie scattering model. The quadratic coefficient fitting values ​​of the aerosol single and multiple scattering albedo are obtained by simulation using the OSOAA radiative transfer model at different aerosol optical depths. S2.3: Using the spectral response correction method, the spectral correction factors generated by the spectral response functions of the two satellites based on the long-term hyperspectral MOBY water-leaving radiance dataset stored in the system are used to construct the spectral correction slope factors and intercept factors between the bands. After the initial generation, these factors are stored. The water-leaving radiance provided by the reference payload is then converted to the water-leaving radiance at the central wavelength of the calibration payload. S2.4: Obtain gas absorption coefficients from the Hitran database. Calculate the transmittance of each gas by accumulating the beam spread function based on the absorption coefficients and combining them with the gas density. Then, use an aerosol model containing aerosol lookup tables (LUTs) to build an aerosol model at different zenith angles, covering a dynamic range of 0° to 80°, to obtain the diffuse transmittance of each gas. The total atmospheric absorption transmittance is then calculated as the product of the absorption coefficients of each gas. S2.5: For the coupled ocean-atmosphere system, in the absence of solar flare and sea whitecaps, the total radiance is calculated by combining it with atmospheric reflected radiation using a radiative transfer model.

2. The universal ocean color satellite on-orbit radiometric calibration method according to claim 1, characterized in that: The quality control principles of the data analysis and quality control module specifically include: For field observation data, the QA scoring system is used for scoring. If the score generated by the QA scoring system is greater than 0.8, the data is directly adopted; otherwise, it is converted into remote sensing reflectance; For satellite image data, when L2 level data is provided, it is screened according to the provided internationally accepted 32-bit flag to eliminate the influence of land, solar flares, saturated radiance, and cloud factors; otherwise, the satellite impact data is screened based on the principles of wind speed less than 6m / s, 5×5 grid top of atmosphere reflectivity variation coefficient within 0.15, solar zenith angle less than 70 degrees, satellite observation angle less than 56 degrees, and the longest near-infrared band top of atmosphere remote sensing reflectivity less than 0.

027.

3. The universal ocean color satellite on-orbit radiometric calibration method according to claim 1, characterized in that: The specific steps of the quality control screening in the fourth step are as follows: S1.1: If the calibration is for system replacement, exclude data with a time difference greater than 1 hour and pixels in satellite images whose distance from the on-site observation site exceeds the spatial resolution of the calibration payload; S1.2: If cross-calibration is used, data with a time difference greater than 30 minutes will be excluded, as will data where the spatial distance between the nearest pixels of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel, and the overlap range of the two payload images is less than 15×15 pixels.

4. The universal ocean color satellite on-orbit radiometric calibration method according to claim 1, characterized in that: The specific calculation formula for the total radiant brightness described in S2.5 is as follows: L t (λ)=(L r (λ)+L a (λ)+L ra (λ)+t(λ)L w (λ))×T g (l) Where, L t (λ) represents the total radiance; L r (λ) represents the Rayleigh scattered radiation generated by air molecules; L a (λ)+L ra (λ) represents aerosol scattered radiation; t(λ) represents the diffuse transmittance from the water surface to the sensor; T g (λ) represents the total atmospheric absorption and transmittance, L w (λ) represents the water-leaving radiance at the central wavelength of the calibration payload.

5. The universal ocean color satellite on-orbit radiometric calibration method according to claim 1, characterized in that: The calibration coefficient comprehensive analysis module is divided into two modes: daily and periodic. In the daily mode, the calibration report displays the process quantity of each calibration algorithm on that day. In the periodic mode, different calibration coefficients obtained by different calibration methods at different times are evaluated to achieve a comprehensive analysis of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within a specified time, and provide a set of calibration coefficients finally adopted.

6. The universal ocean color satellite on-orbit radiometric calibration method according to claim 5, characterized in that: The specific calculation steps of the final calibration coefficient are as follows: S3.1: Compare the relative percentage differences between the calibration coefficients obtained by different calibration methods or the calibration coefficients obtained by the same method within a specified time period. If the relative percentage difference is less than 5%, arrange all calibration gain coefficients in ascending order and calculate the average of the semi-quartile range as the final calibration coefficient; S3.2: If the relative percentage difference is greater than 5%, then based on the stability of the trend over time, remove the abnormal data points that deviate from the mean of the calibration gain coefficient by one standard deviation, and recalculate the relative percentage difference to obtain the final calibration coefficient.

Images