Universal water color satellite in-orbit radiometric calibration method

Abstract

The invention discloses a universal water color satellite in-orbit radiometric calibration method, which belongs to the field of radiometric calibration and comprises a data analysis and quality control module, a radiometric calibration module, a calibration coefficient comprehensive analysis module and a working state recording module. According to the invention, the algorithm of the hyperspectral atmospheric absorption coefficient is optimized, and the compatibility of the on-orbit radiometric calibration algorithm to various types of water color satellite loads is realized.

Description

Technical Field

[0001] The present invention relates to the field of radiometric calibration, and particularly to a general water color satellite on-orbit radiometric calibration method. 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 important scientific and economic benefits to our in-depth understanding of the ocean. In order to obtain these high-quality optically active component products, the uncertainty of ocean radiation measurement (such as water-leaving radiance Lw or remote sensing reflectance Rrs) needs to be controlled within 5% through atmospheric correction (Hooker et al., 1992). However, since Lw in the visible (VIS) band only accounts for 10% or less of the total radiance (Lt) measured by ocean color sensors at the top-of-atmosphere (TOA), the uncertainty of Lt is required to be less than 0.5%. To meet this requirement, ocean color sensors need to perform a series of radiometric calibration tasks to couple and reduce the uncertainty brought by atmospheric correction and ocean color inversion algorithms. Although existing publicly available OCRS on-orbit radiometric calibration systems already include the above-mentioned multiple calibration methods, these systems are usually only applicable to specific sensors, are difficult to be directly applied to other similar payloads, and are even less able to meet the requirements of future hyperspectral and high spatial resolution sensors. Given the operating conditions of numerous OCRS systems internationally and the upcoming launch of more advanced ocean color sensors, this study aims to develop a general on-orbit radiometric calibration system, namely the general ocean color radiometric calibration software. Summary of the Invention

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

[0004] To achieve the above purpose, the present invention adopts the following technical solutions: A general water color satellite on-orbit radiometric calibration method includes a data parsing and quality control module, a radiometric calibration module, a calibration coefficient comprehensive analysis module, and a working state recording module; The specific working steps of the above-mentioned general water color satellite on-orbit radiometric calibration method based on the above modules are as follows: The first step: Input the data into the data parsing and quality control module, and parse the input calibration payload data and various types of auxiliary data; The second step: Select the corresponding calibration algorithm according to the input data; Step 3: When it is cross-calibration, the scheduling radiometric calibration module dynamically selects the calibration area; Step 4: Conduct quality control screening, extract pixel points applicable to the calibration method and meeting the matching criteria, and generate the corresponding matching data set; Step 5: The data parsing and quality control module outputs scheduling instructions and the matching data set, and calculates the calibration coefficients according to different calibration algorithms through the radiometric calibration module; Step 6: Summarize the calibration coefficients, generate the corresponding daily or periodic report and suggestions for adopting the calibration coefficients through the calibration coefficient comprehensive analysis module, and then the working status recording module records the system working information, tracks the working status of other modules, and generates a work log.

[0005] As a further solution of the present invention, the quality control principles of the data parsing and quality control module specifically include: For on-site observation data, use the QA scoring system to score. 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, screen according to the provided international general 32-bit flag, and eliminate the influence of various factors such as land, solar flare, saturated radiance, and cloud; otherwise, according to the principles that the wind speed is less than 6 m / s, the coefficient of variation of the top-of-atmosphere reflectance of the 5×5 grid is within 0.15, the solar zenith angle is less than 70 degrees, the satellite observation angle is less than 56 degrees, and the top-of-atmosphere remote sensing reflectance of the longest near-infrared band is less than 0.027, screen the satellite impact data.

[0006] As a further solution of the present invention, the specific steps of the quality control screening in Step 4 are as follows: S1.1: If it is system substitution calibration, eliminate the data with a time difference greater than 1 hour, and the pixel points in the satellite image whose distance from the on-site observation site exceeds the spatial resolution of the calibration payload; S1.2: If it is cross-calibration, eliminate the data with a time difference greater than 30 minutes, and the spatial distance between the nearest pixel points of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel point, and eliminate it when the overlapping range of the two payload images is less than 15×15 pixels.

[0007] As a further solution of the present invention, the specific steps of calculating the calibration coefficients according to different calibration algorithms through the radiometric calibration module in Step 5 are as follows: S2.1: Establish hyperspectral Rayleigh LUTs based on OSOAA, simulate different solar zenith angles, satellite zenith angles and wind speeds, obtain the three coefficients of the Stokes vector of Rayleigh reflection, and perform Fourier expansion on the relative azimuth angle according to the obtained groups of Stokes vector coefficients; S2.2: Based on the 80 groups of aerosol models divided by Ahmad, calculate and obtain the single-scattering albedo and scattering phase function of aerosols through the Mie scattering model. Through the OSOAA radiative transfer model, simulate under different aerosol optical depths to obtain the fitting values of the quadratic term coefficients of the single and multiple scattering albedos of aerosols; S2.3: Adopt the spectral response correction method to construct the spectral correction slope factor and intercept factor between bands through the spectral correction factor generated by two satellite spectral response functions based on the hyperspectral MOBY water-leaving radiance dataset stored in the system over a long time range, and store it after the first generation. Then convert the water-leaving radiance provided by the reference payload into the water-leaving radiance at the central wavelength of the calibration payload; S2.4: Obtain the gas absorption coefficients from the Hitran database, calculate the transmittance of each gas through the absorption coefficient cumulative beam distribution function and in combination with the gas density. Then adopt the aerosol model containing aerosol LUTs, established under different zenith angle conditions, covering the dynamic range from 0° to 80°, to obtain the diffuse transmittance of each gas. Then obtain the total absorption transmittance of the atmosphere according to the product of the absorption coefficients of each gas; S2.5: For the coupled ocean-atmosphere system, in the case of no solar glint and sea surface whitecaps, simulate and calculate by combining with the atmospheric reflected radiation through the radiative transfer model to obtain the total radiance.

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

[0009] In the formula, represents the total radiance; represents the Rayleigh scattering radiation generated by the scattering of air molecules; represents the aerosol scattering radiation; represents the diffuse transmittance from the water surface to the sensor; represents the total absorption transmittance of the atmosphere.

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

[0011] In the formula, represents the water-leaving radiance of the corresponding band of the reference source; represents the spectral response function of the reference source; represents the spectral response function of the payload to be calibrated; The specific calculation formula of the total absorption transmittance of the atmosphere described in S2.4 is as follows:

[0012]

[0013]

[0014]

[0015]

[0016]

[0017]

[0018] wherein, is calculated from the of the absorption coefficient; represents the gas density; represents the interval, which can be selected according to the bandwidth of the calibration load band; and respectively represent the gas transmittance from the sun to the sea surface and from the sea surface to the satellite; and respectively represent the solar zenith angle and the satellite observation angle; and respectively represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite.

[0019] 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 shows the process quantities of each calibration algorithm on the current day. In the periodic mode, evaluation is performed based on the different calibration coefficients obtained by different calibration methods within different time periods, so as to realize the comprehensive analysis of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within a specified time period, and a set of finally adopted calibration coefficients is given.

[0020] As a further solution of the present invention, the specific calculation steps of the finally adopted calibration coefficient are as follows: S3.1: Compare 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 period. If the relative percentage difference is less than 5%, then arrange all the calibration gain coefficients in ascending order, and calculate the average value of the semi-interquartile range as the final calibration coefficient; S3.2: If the relative percentage difference is greater than 5%, then eliminate the abnormal data points that are one standard deviation away from the mean of the calibration gain coefficient deviation according to the stability of the trend over time, and calculate the relative percentage difference again to obtain the final calibration coefficient.

[0021] Compared with the prior art, the beneficial effects of the present invention are as follows: The present invention optimizes the algorithm for the atmospheric absorption coefficient of hyperspectral by replacing the single look-up table in various calibration algorithms with a hyperspectral Rayleigh and aerosol look-up table at 1nm intervals established based on the OSOAA radiative transfer model and long-term observational data of multiple types of water body sites in Ahmad et al.'s AERONET-OC. At the same time, a spectral correction factor for the spectral response difference of the water-leaving radiance of each sensor is constructed for the cross-calibration algorithm, realizing the compatibility of the on-orbit calibration algorithm for multiple types of ocean color satellite payloads. Description of the Drawings

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

[0023] Figure 1 It is a system block diagram of a general ocean color satellite on-orbit calibration method proposed by the present invention; Figure 2 It is an operation flowchart of a general ocean color satellite on-orbit calibration method proposed by the present invention.

[0024] Figure 3 It is a 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 look-up table; circles represent the MODIS look-up table); Figure 4 It is a spectral variation of the extinction coefficient and scattering phase function (scattering angle is 2°) under different relative humidities (triangles represent the convolution results of the hyperspectral Rayleigh look-up table; circles represent the MODIS look-up table); Figure 5 It is a scatter plot comparison of the water-leaving radiance of the central bands after spectral response correction of Sentinel-3A OLCI and HY3A COCTS2; Figure 6 It is a time series diagram of the on-orbit calibration coefficients of COCTS-HY1C / D; Figure 7 It is a spectral diagram of the cross-calibration coefficients of COCTS-HY1C / D. Detailed Embodiments

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

[0026] Referring to Figure 1-2 , a general ocean color satellite on-orbit calibration method includes a data parsing and quality control module, a radiation calibration module, a calibration coefficient comprehensive analysis module, and a working status recording module.

[0027] Among them, the data parsing and quality control module reads and organizes according to the dedicated parsing units of each data source. At the same time, it determines the calibration method based on the input data source, and then performs quality control on the unified format data. The quality control criteria adopt corresponding data matching criteria for different data sources.

[0028] For in-situ observation data (MOBY), data with the highest quality assessment level is directly adopted. Otherwise, it is converted into remote sensing reflectance, and a QA scoring system (Wei, J., Lee, Z.-P., and Shang, S. (2016). Asystem to measure the data quality of spectral remote sensing reflectance ofaquatic environments. J. Geophys. R., 121, 8189–8207.) is used for scoring, and data with a score greater than 0.8 is adopted. For satellite image data, when providing L2-level data, it is screened according to the provided international standard 32-bit flag, and factors such as land, solar flare, saturated radiance, and cloud are excluded; otherwise, it is screened according to the following principles: 1) wind speed is less than 6 m / s; 2) the coefficient of variation of the top-of-atmosphere reflectance within a 5×5 grid is within 0.15; 3) solar zenith angle is less than 70°; 4) satellite observation angle is less than 56°; 5) the top-of-atmosphere remote sensing reflectance of the longest near-infrared band is less than 0.027.

[0029] For pixel points that meet the quality control conditions, the corresponding spatio-temporal matching criteria are selected according to the calibration algorithm.

[0030] System substitution calibration: 1) Data with a time difference greater than 1 hour is excluded; 2) Pixel points in the satellite image that are more than the spatial resolution of the calibration payload away from the in-situ observation site are excluded.

[0031] Cross calibration: 1) Data with a time difference greater than 30 minutes is excluded; 2) Exclude when the overlapping range of the two payload images is less than 15×15 pixels; 3) Exclude when the spatial distance between the nearest pixel points of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel point.

[0032] The radiometric calibration module includes four algorithms: on-orbit calibration, system substitution calibration, cross calibration, and land calibration. The principles of each algorithm are introduced as follows: a) On-orbit calibration The on-orbit calibration method of this system is solar calibration. The main calibration devices of the on-orbit calibration system include a solar diffuser (SD), a solar diffuser stability monitor (SDSM), a solar attenuator, a rotating telescope assembly (RTA), and a half-angle mirror (HAM), etc.

[0033] Sunlight is attenuated by the solar attenuation screen and then enters the solar diffuser. After that, it passes through the rotating telescopic device and is reflected by the half-angle mirror into the calibration sensor. At this time, the solar irradiance can be expressed as:

[0034] where: represents the solar irradiance at the entrance pupil when using the half-angle mirror surface; is the extraterrestrial solar irradiance; represents the sun-earth correction factor; and are the solar zenith angle and azimuth angle in the SD coordinate system, respectively; and are the observed zenith angle and azimuth angle of the calibration payload observing SD, respectively; is the bidirectional reflectance distribution function of the solar diffuser; and are the solar zenith angle and azimuth angle incident on the solar attenuation screen, respectively; is the transmittance of the solar attenuation screen; is the BRDF correction coefficient of the solar diffuser.

[0035] The irradiance value entering the calibration payload can be calculated by the above formula (1). On this basis, the solar calibration gain coefficient in each band can be finally calculated by the following formula .

[0036]

[0037] where: is the DN value when observing the solar diffuser using the half-angle mirror surface, is the DN value when observing the cold sky using the half-angle mirror surface.

[0038] b) System substitution calibration The system substitution calibration and cross-calibration gain coefficient are defined as the ratio of the theoretical radiance at the top of the atmosphere to the actual observed radiance at the top of the atmosphere, as shown in formula (3).

[0039]

[0040] The actual observed radiance at the top of the atmosphere is obtained by combining the DN value observed by the calibration payload data and the absolute radiometric calibration coefficient calculated in the laboratory before launch according to formula (3), and then converted into the reflectance at the top of the atmosphere through formula (4). The theoretical observed radiance at the top of the atmosphere is derived to the top of the atmosphere through the system substitution calibration algorithm in combination with the on-site measured water body information and meteorological auxiliary data.

[0041]

[0042] Wherein: is the solar irradiance; is the sun-earth corrected distance; is the solar zenith angle.

[0043] The system substitution calibration algorithm is usually expressed as the inverse process of the atmospheric correction algorithm. The atmospheric correction algorithm refers to the process of removing the atmospheric path reflectance from the top-of-atmosphere reflectance received by the water color satellite sensor to obtain the water body signal. However, the atmospheric correction algorithm itself is affected by various factors such as the optical properties of the water body, the atmosphere, and the sun-satellite observation geometry, which causes errors, making the uncertainty of the retrieved water-leaving radiance increase beyond the limit range. By introducing the atmospheric correction algorithm into the system substitution calibration algorithm, the uncertainty of the atmospheric correction can be coupled during the calibration process, minimizing the deviation between the remote sensing reflectance retrieved from the top-of-atmosphere reflectance and the measured data. The cross-calibration algorithm process is the same as the system substitution calibration algorithm process.

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

[0045] Wherein: and represent the Rayleigh scattering and aerosol (including the interaction between molecules and aerosols) reflectances respectively; represents the white cap reflectance; represents the remote sensing reflectance; and represent the atmospheric diffuse transmittances from the sun to the water body and from the water body to the sensor respectively; represents the atmospheric transmittance.

[0046] It can be seen from formula (5) that the top-of-atmosphere reflectance includes the contributions of physical quantities such as Rayleigh scattering, aerosol scattering, water body remote sensing reflectance, and sea surface white caps. The contribution generated by solar flares is ignored in the above calculation formula. To avoid errors, pixels with values greater than 0 are removed by normalizing the solar flare coefficient. Next, for each parameter.

[0047] First of all, the sea surface white cap is formed by the breaking of sea waves caused by the wind on the sea surface, and affects the top-of-atmosphere reflectance by reflecting sunlight and skylight. It is mainly manifested as the relative relationship with the sea surface wind speed. When the wind speed is less than 6.33 m / s, the white cap reflectance of the water body is 0, and the specific calculation formula is shown in (6). Among them, is the normalized white cap reflectance, is the sea surface wind speed.

[0048]

[0049] To simplify the calculation process of the algorithm and improve the calculation efficiency, both the system substitution calibration and the cross-calibration adopt the lookup table method for lookup calculation. The Rayleigh lookup table is divided into multiple files by wavelength bands. It contains parameters such as Rayleigh optical depth, solar zenith angle, satellite observation angle, depolarization factor, wind speed, and Rayleigh reflectance coefficient. According to the wind speed provided by the auxiliary data, the atmospheric pressure, and the observation geometry provided by the satellite data to be calibrated, the calculation can be carried out. .

[0050] The calculation of atmospheric transmittance requires atmospheric pressure, Rayleigh optical depth, and observation geometry. Among them, the atmospheric pressure and the observation geometry are provided by the auxiliary data and the data to be calibrated respectively. The Rayleigh optical depth The calculation is shown in formula (7). The transmittance of a single atmospheric gas is expressed as the product of the atmospheric gas content and the gas absorption coefficient. To correct the influence of the observation geometry on the atmospheric transmittance, the specific calculation formula is shown in (8).

[0051]

[0052]

[0053] Where: are the current atmospheric pressure and the standard atmospheric pressure respectively; is the satellite observation angle; represents the gas type.

[0054] The system substitution calibration separates the calculation processes of the near-infrared to short-wave infrared and ultraviolet to visible light bands and calculates them in two steps: The first step is based on the "dark pixel" hypothesis theory and combines the known water aerosol types to achieve near-infrared to short-wave infrared calibration; the second step is based on the atmospheric correction algorithm and extrapolates the calibration gain coefficient of the visible light band through two near-infrared / short-wave infrared bands. The following specifically describes the calculation of the aerosol scattering contribution based on the above process.

[0055] Based on the two prerequisites of the system substitution calibration algorithm, the South Pacific Circulation Region and the South Indian Ocean Circulation Region waters with cleaner, more uniform, and stable non - coastal marine - type atmospheric conditions on Earth are selected in this paper for near - infrared and short - wave infrared calibration. It is considered that the aerosol types in the two circulation regions are a mixed aerosol model with 50% each of marine - type aerosol M70 with a relative humidity of 70% and marine - type aerosol model M90 with a relative humidity of 90%. First, assuming that the initial system calibration gain coefficient of the longer short - wave infrared band (1640 nm) is 1 in the waters with known aerosol types, the aerosol reflectance of other short - wave infrared bands (1245 nm) can be estimated by combining the aerosol lookup table. Subsequently, the aerosol optical properties in the near - infrared band are estimated through the "dark pixel" atmospheric correction algorithm based on the band infrared as a reference. Finally, the near - infrared - short - wave infrared calibration gain coefficient can be obtained based on the calculation results of the above other components. At this time, the system substitution calibration gain coefficient of the 1640 - nm band is 1. Since the system substitution calibration algorithm assumes that the longer near - infrared band (865 nm) is 1, the calibration gain coefficient of the 1640 - nm band is adjusted according to the system substitution calibration gain coefficient of 865 nm at this time. The above process is repeated for multiple iterations until the system substitution calibration gain coefficient of the 865 - nm band is 1. The finally obtained near - infrared - short - wave infrared system substitution calibration gain coefficient is the final value.

[0056] Subsequently, the near - infrared - band infrared system substitution calibration gain coefficient is used for the waters of the in - situ observation site. At this time, based on the calibration gain coefficients of the two near - infrared / short - wave infrared bands, the aerosol type and the multiple - scattering reflectance of the aerosol in this area are determined, and then the single - scattering reflectance of the near - infrared / short - wave infrared aerosol is calculated by combining the lookup table. . Finally, based on the aerosol single - scattering reflectance, the aerosol optical depth conversion formula, the extinction coefficients of different bands and the aerosol single - multiple conversion, the ultraviolet - visible aerosol scattering reflectance can be calculated.

[0057]

[0058]

[0059] Where: is the aerosol optical depth; is the scattering phase function; is the aerosol single - scattering albedo.

[0060] The remote - sensing reflectance data is provided by the in - situ measurement site, and the calculation formula is as follows.

[0061]

[0062] Where: Provide the remote sensing reflectance obtained by converting the water-leaving radiance for the ocean optical buoy, which is the spectral response function of the calibration payload. which is the remote sensing reflectance corresponding to the central wavelength of the calibration payload.

[0063] c) Cross-calibration The cross-calibration algorithm is different from the system substitution calibration algorithm for calculating the aerosol scattering reflectance. It is necessary to use a reference payload to provide the aerosol optical depth and aerosol type of the reference band in the cross region. The estimation of the aerosol reflectance in the ultraviolet-visible band is based on the estimation of the aerosol scattering reflectance in the near-infrared band based on the "dark pixel" hypothesis. The process is as follows: The aerosol single-scattering reflectance is estimated by using the aerosol lookup table and the top-of-atmosphere reflectance of two near-infrared bands observed by the reference payload. Since the aerosol single-scattering reflectance is related to the inherent properties of the aerosol, such as: scattering phase function, single-scattering albedo, and attenuation coefficient of the aerosol. Therefore, the aerosol single-scattering reflectance of each aerosol model can be calculated, and the most suitable aerosol model can be selected to make the calculated radiation amount match the measured value best. Based on directly determining the aerosol type through statistical analysis, the estimation of the aerosol reflectance in the ultraviolet-visible band is directly extrapolated from the two known aerosol models.

[0064] The remote sensing reflectance is directly provided by the reference payload and obtained by performing an atmospheric correction program on the reference payload data. At this time, considering the deviation of the reference band from the central wavelength of the payload to be calibrated and the difference in spectral response, the band offset algorithm based on the bio-optical model can be used to select the closest bands of the calibration payload to the reference payload bands for conversion.

[0065] d) Land calibration The 6S radiative transfer model is selected for simulation in land calibration. Assuming cloud-free atmospheric conditions, the top-of-atmosphere radiance can be expressed as:

[0066]

[0067] where 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 atmospheric spherical albedo, is the transmittance formed by the atmospheric molecular absorption, are the diffuse transmittances on the two paths from the sun and the sensor to the ground target respectively, and are the top-of-atmosphere radiance and the top-of-atmosphere reflectance respectively, is the cosine of the solar zenith angle, is the top-of-atmosphere solar irradiance, is the Earth-Sun distance (unit: AU).

[0068] The top-of-atmosphere radiance in the visible and ultraviolet bands can be calculated and obtained through the 6S radiative transfer model . The land substitution calibration coefficient involves the differences between different detectors. The calculation formula is:

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

[0070] Among them: and are the look-up table parameters of the hyperspectral and calibration payload center bands respectively; is the spectral response function.

[0071] The atmospheric transmittance uses the gas absorption spectral lines provided by the Hitran database and calculates the gas transmittance line by line using the Line-by-line (LBL) method. The K-distribution method is introduced to improve the calculation efficiency of LBL. First, obtain the gas absorption coefficient from the Hitran database. Calculate the transmittance of each gas by combining the absorption coefficient cumulative beam distribution function with the gas density. The calculation formula is as follows.

[0072]

[0073]

[0074] is calculated from the of the absorption coefficient. is the gas density (provided by auxiliary data). The interval can be selected according to the bandwidth of the calibration payload band. From the above formula, the absorption transmittances of different gases can be obtained, and combined with the observation geometry, the total gas transmittance can be obtained:

[0075]

[0076]

[0077]

[0078]

[0079]

[0080] wherein: and respectively represent the gas transmittance of the atmosphere from the sun to the sea surface and from the sea surface to the satellite; and respectively represent the solar zenith angle and the satellite viewing angle; and respectively represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite.

[0081] In the cross - calibration algorithm, due to the spectral response differences between the calibration payload and the reference payload, considering the uncertainties and computational efficiency brought about by the uncorrected spectral function differences in the band - shift algorithm based on the bio - optical model, therefore, based on the spectral response correction method proposed by Doelling et al. (Doelling, D. R., et al. (2011). Spectral Reflectance Corrections for Satellite Intercalibrations Using SCIAMACHY Data. IEEE Geoscience and Remote Sensing Letters, 9(1), 119 - 123.), long - term high - spectral data measured in - situ by the Marine Optical Buoy MOBY is used to construct the spectral correction slope factor and the intercept factor , so as to convert the water - leaving radiance provided by the reference payload into the water - leaving 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 two satellite spectral response functions to generate spectral correction factors based on the long - term high - spectral water - leaving radiance dataset of MOBY stored in the system, and the generated factors will be stored after the first generation. Formula (23) shows the usage of the spectral correction slope factor and the intercept factor .

[0082]

[0083] where is the radiance leaving water corresponding to the reference load band, reference load spectral response function, is the spectral response function of the load to be calibrated.

[0084] A general method for on-orbit radiometric calibration of water color satellites based on the above modules has the following specific working steps: The first step: Input the data into the data parsing and quality control module to parse the input calibration load data and various auxiliary data.

[0085] It should be further noted that the quality control principles of the data parsing and quality control module specifically include: For in-situ observation data, a 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 international common 32-bit flag, and the influences of land, solar flare, saturated radiance, cloud and other factors are eliminated; otherwise, according to the principles that the wind speed is less than 6 m / s, the coefficient of variation of the top-of-atmosphere reflectance within a 5×5 grid is less than 0.15, the solar zenith angle is less than 70 degrees, the satellite observation angle is less than 56 degrees, and the top-of-atmosphere remote sensing reflectance of the longest near-infrared band is less than 0.027, the satellite image data is screened.

[0086] The second step: Select the corresponding calibration algorithm according to the input data.

[0087] The third step: When it is cross-calibration, dispatch the radiometric calibration module to dynamically select the calibration area.

[0088] The fourth step: Conduct quality control screening, extract the pixel points applicable to the calibration method and meeting the matching criteria, and generate the corresponding matching data set.

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

[0090] The fifth step: The data parsing and quality control module outputs the scheduling instructions and the matching data set, and calculates the calibration coefficients through the radiometric calibration module according to different calibration algorithms.

[0091] Specifically, based on OSOAA, hyperspectral Rayleigh LUTs are established. By simulating different solar zenith angles, satellite zenith angles, and wind speeds, the three coefficients of the Stokes vector of Rayleigh reflection are obtained. Based on the obtained Stokes vector coefficients for each group, Fourier expansion is performed on the relative azimuth angle. Based on the 80 groups of aerosol models divided by Ahmad, the single-scattering albedo and scattering phase function of aerosols are calculated through the Mie scattering model. Through the OSOAA radiative transfer model, simulations are carried out at different aerosol optical depths to obtain the fitting values of the quadratic term coefficients of the single and multiple scattering albedos of aerosols. The spectral response correction method is adopted. Based on the spectral correction factors generated from the hyperspectral MOBY water-leaving radiance dataset stored in the system over a long time range by two satellite spectral response functions, the spectral correction slope factor and intercept factor between bands are constructed and stored after the first generation. Then, the water-leaving radiance provided by the reference payload is converted into the water-leaving radiance at the central wavelength of the calibration payload. The gas absorption coefficients are obtained from the Hitran database, and the transmittance of each gas is calculated by accumulating the beam distribution function of the absorption coefficient and combining with the gas density. Then, an aerosol model containing aerosol LUTs is adopted, established under different zenith angle conditions, covering a dynamic range from 0° to 80°, to obtain the diffuse transmittance of each gas. Then, the total absorption transmittance of the atmosphere is obtained based on the product of the absorption coefficients of each gas. For the coupled ocean-atmosphere system, in the absence of solar glint and sea surface whitecaps, the total radiance is simulated and calculated by combining with the atmospheric reflected radiation through the radiative transfer model.

[0092] In this embodiment, the specific calculation formula for the total radiance is as follows:

[0093] In the formula, represents the total radiance; represents the Rayleigh scattering radiation generated by the scattering of air molecules; represents the aerosol scattering radiation; represents the diffuse transmittance from the water surface to the sensor; represents the total absorption transmittance of the atmosphere.

[0094] The specific calculation formula for the water-leaving radiance at the central wavelength of the calibration payload described in S2.3 is as follows

[0095] In the formula, represents the water-leaving radiance of the corresponding band of the reference source; represents the spectral response function of the reference source; represents the spectral response function of the payload to be calibrated; The specific calculation formula for the total absorption transmittance of the atmosphere is as follows:

[0096]

[0097]

[0098]

[0099]

[0100]

[0101]

[0102] In the formula, is calculated from of the absorption coefficient; represents the gas density; represents the interval, which can be selected according to the bandwidth of the calibration load band; and respectively represent the gas transmittance from the sun to the sea surface and from the sea surface to the satellite; and respectively represent the solar zenith angle and the satellite observation angle; and respectively represent the total gas transmittance from the sun to the sea surface and from the sea surface to the satellite.

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

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

[0105] Step 6: Summarize the calibration coefficients. Through the calibration coefficient comprehensive analysis module, generate the corresponding daily or periodic reports, as well as suggestions for adopting the calibration coefficients. After that, the working status recording module records the system working information, tracks the working status of other modules, and generates a work log.

[0106] It should be further noted that the calibration coefficient comprehensive analysis module is divided into two modes: daily and periodic. In the daily mode, the calibration report shows the process quantities of each calibration algorithm on the current day. In the periodic mode, evaluations are made based on the different calibration coefficients obtained by different calibration methods within different time periods, so as to achieve the comprehensive analysis of the calibration coefficients obtained by different calibration methods or the calibration coefficients of the same method within the specified time period, and give a set of finally adopted calibration coefficients.

[0107] In addition, it should be noted that the finally adopted calibration coefficient 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-interquartile range is calculated as the final calibration coefficient. If the relative percentage difference is greater than 5%, abnormal data points with a deviation mean of one standard deviation of the calibration gain coefficient are excluded based on the stability of the trend over time, and the relative percentage difference is calculated again to obtain the final calibration coefficient.

[0108] To verify the reliability of the above-established hyperspectral aerosol and Rayleigh lookup tables, they are verified here with the MODIS-specific lookup tables in the NASA operational system. First, the consistency of the three parameters of the Rayleigh reflectance with the corresponding parameters of MODIS is shown. The following shows the spectral curve comparison under different conditions of wind speed of 2 m / s, solar zenith angle (solz), and satellite zenith angle (senz). It can be found that the spectral similarity of the Rayleigh reflectance is 1 at the same observation angle, and the mean absolute relative error (MARD) of each parameter is better than 1%, as shown in Table 3; Table 3 Statistical parameters of spectral curves under different observation geometries (wind speed is 2 m / s)

[0109] At the same time, to verify the feasibility of the aerosol optical depth, the spectral curve distributions of the extinction coefficient and scattering phase function (scattering angle is 2°) at the median value of the fine mode fraction for eight relative humidity models are shown here. It can be Figure 3 seen that the extinction coefficient and scattering phase function obtained by hyperspectral convolution are consistent with those in the MODIS lookup table in terms of spectral shape. The spectral similarity of the extinction coefficient and scattering phase function at each relative humidity is close to 1, and the MARD of each parameter is better than 0.5%, showing a high correlation; Table 4 Statistical parameters of spectral curves of extinction coefficient and scattering phase function (scattering angle is 2°) under relative humidity conditions

[0110] In the verification of the algorithm for establishing the spectral correction factor, we took HY3A COCTS2 and Sentinel-3A OLCI as examples for inspection. Considering that the water-leaving radiance values in the bands after 620 nm are relatively small, only the scatter plots of the water-leaving radiance in the 360 nm - 565 nm band are shown here. Although the minimum band of Sentinel-3A OLCI is the 400 nm band, the water-leaving radiance in the ultraviolet band of HY3A COCTS2 can still be well characterized after the water-leaving radiance in the 400 nm band is corrected by the correction factor. The correlation coefficient is greater than 0.995, and the correlation in other visible light bands is better than 0.999. The MAPD is better than 1% when the deviation in the central band is less than 20 nm, and it also shows good consistency when the deviation in the central band is 40 nm, thus verifying the effectiveness of the algorithm.

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

[0112] Subsequently, based on this system, the on-orbit calibration and cross-calibration results of COCTS-HY1C / D were realized, and the time range of the test data was from August 1, 2022 to August 1, 2023. First, the calibration coefficients were obtained by cross-calibrating the on-orbit calibration spectrometer and COCTS.

[0113] From Figure 6 it can be obtained that the on-orbit calibration coefficients of COCTS-HY1C at 490, 520, and 670 nm during this period are 0.99, 0.98, and 0.94 respectively, the standard deviations are 0.009, 0.010, and 0.019 respectively, and the coefficients of variation are 0.88%, 1.06%, and 2.00% respectively; the on-orbit calibration coefficients of COCTS-HY1D at 490, 520, and 670 nm are 0.97, 0.97, and 0.95 respectively, and the coefficients of variation are 1.94%, 1.95%, and 2.44% respectively; To verify the stability of the calibration coefficients, considering the attenuation of the sensor over time, a stable radiation performance was established for COCTS-HY1C / D through the above on-orbit calibration coefficients. On this basis, cross-calibration was carried out, and the coefficients of variation of each band within one year were basically better than 0.5%. The obtained results are as Figure 7 shown in Table 6. Thus, the applicability of the Grace-OC system when applied to the two sensors of COCTS-HY1C / D was verified, and at the same time, the effectiveness and stability of the on-orbit calibration and cross-calibration algorithms were verified; Table 6 Statistical parameters of the cross-calibration of COCTS-HY1C / D

[0114] The above has described in detail an embodiment of the present invention. However, the content is only a preferred embodiment of the present invention and cannot be considered as limiting the scope of implementation of the present invention. All equivalent changes and improvements made according to the scope of the application of the present invention shall still fall within the scope covered by the patent of the present invention.

Claims

1. A general method for on-orbit radiometric calibration of ocean color satellites, characterized in that, It includes a data parsing and quality control module, a radiometric calibration module, a calibration coefficient comprehensive analysis module, and a working status recording module; The specific working steps of the above-mentioned modules for the universal on-orbit radiometric calibration method of water color satellites are as follows: The first step: Input the data into the data parsing and quality control module to parse the input calibration payload data and various auxiliary data; The second step: Select the corresponding calibration algorithm according to the input data; The third step: When it is cross-calibration, schedule the radiometric calibration module to dynamically select the calibration area; The fourth step: Conduct quality control screening, extract the pixel points applicable to the calibration method and meeting the matching criteria, and generate the corresponding matching data set; The fifth step: The data parsing and quality control module outputs the scheduling instructions and the matching data set, and calculates the calibration coefficients through the radiometric calibration module according to different calibration algorithms; The sixth step: Summarize the calibration coefficients, generate the corresponding daily or periodic reports and calibration coefficient adoption suggestions through the calibration coefficient comprehensive analysis module, and then the working status recording module records the system working information, tracks the working status of other modules, and generates a work log.

2. The universal on-orbit radiometric calibration method for ocean color satellites according to claim 1, characterized in that, The specific quality control principles of the data parsing and quality control module are as follows: For field observation data, use the QA scoring system to score. If the score generated by the QA scoring system is greater than 0.8, directly adopt the data; otherwise, convert it into remote sensing reflectance; For satellite image data, when L2-level data is provided, screen according to the provided international general 32-bit flag, and eliminate the influence of various factors such as land, solar flare, saturated radiance, and cloud; otherwise, screen the satellite image data according to the principles that the wind speed is less than 6m / s, the coefficient of variation of the top-of-atmosphere reflectance in a 5×5 grid is within 0.15, the solar zenith angle is less than 70 degrees, the satellite viewing angle is less than 56 degrees, and the top-of-atmosphere remote sensing reflectance of the longest near-infrared band is less than 0.

027.

3. A general method for on-orbit radiometric calibration of ocean color satellites according to claim 1, characterized in that, The specific steps of the quality control screening described in the fourth step are as follows: S1.1: If it is system substitution calibration, eliminate the data with a time difference greater than 1 hour, and the pixel points in the satellite image whose distance from the field observation site exceeds the spatial resolution of the calibration payload; S1.2: If it is cross-calibration, eliminate the data with a time difference greater than 30 minutes, and the pixel points to be eliminated when the spatial distance between the nearest pixel points of the calibration payload and the reference payload is greater than the minimum spatial resolution pixel point and the overlapping range of the two payload images is less than 15×15 pixels.

4. A general water color satellite on-orbit radiometric calibration method according to claim 3, characterized in that The specific steps of calculating the calibration coefficients through the radiometric calibration module according to different calibration algorithms described in the fifth step are as follows: S2.1: Establish hyperspectral Rayleigh LUTs based on OSOAA, simulate different solar zenith angles, satellite zenith angles, and wind speeds to obtain the three coefficients of the Stokes vector of Rayleigh reflection, and perform Fourier expansion on the relative azimuth angle according to the obtained groups of Stokes vector coefficients; S2.2: Based on the 80 groups of aerosol models divided by Ahmad, calculate the single-scattering albedo and scattering phase function of aerosols through the Mie scattering model, and simulate through the OSOAA radiative transfer model at different aerosol optical depths to obtain the quadratic coefficient fitting values of the single and multiple scattering albedos of aerosols; S2.3: The spectral response correction method is adopted to construct the spectral correction slope factor and intercept factor between the bands by using the spectral response function of two satellites based on the spectral correction factor generated by the hyperspectral MOBY water-leaving radiance data set stored in the system under a long time range, and store them after the first generation. After that, the water-leaving radiance provided by the reference payload is converted into the water-leaving radiance at the central wavelength of the calibration payload. S2.4: Obtain gas absorption coefficients from the Hitran database, accumulate beam spread functions through absorption coefficients, and calculate the transmittance of each gas in combination with gas density. Then, an aerosol model including aerosol LUTs is established under different zenith angle conditions, covering a dynamic range of 0° to 80°, to obtain the diffuse transmittance of each gas. Then, the total absorption transmittance of the atmosphere is obtained based on 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 radiative transfer model is used to combine the atmospheric reflected radiation to obtain the total radiance.

5. A general water color satellite on-orbit radiometric calibration method according to claim 4, characterized in that, The specific calculation formula for the total radiant brightness described in S2.5 is as follows: In the formula, represents the total radiance; represents the Rayleigh scattering radiation generated by air molecule scattering; represents the aerosol scattering radiation; represents the diffuse transmittance from the water surface to the sensor; represents the total atmospheric absorption transmittance.

6. A general water 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.

7. A general water color satellite on-orbit radiometric calibration method according to claim 6, characterized in that The specific calculation steps of the calibration coefficient finally adopted are as follows: S3.1: Compare 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 order from small to large, and the average value of the semi-quartile range is calculated as the final calibration coefficient; S3.2: If the relative percentage difference is greater than 5%, then the 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.

Images