An on-orbit automatic updating method for a remote sensing instrument reflection band calibration parameter

Abstract

The present application relates to the field of communication technology, provide a kind of remote sensing instrument reflection band calibration parameter on-orbit automatic updating method, the deviation of observed reflectivity is calculated using the simulated reflectivity of multi-stable target, the reliability of processing result can be improved by two layers of data quality control strategy, and for the annual cycle variation of observed reflectivity relative deviation, the corresponding annual daily correction model is established, and relatively accurate observed reflectivity deviation can be obtained;For the annual cycle variation of calibration slope, the corresponding annual daily correction model is established, and relatively accurate calibration slope and decay rate can be obtained.The present application provides a kind of time continuous, stable on-orbit observation deviation calculation, and on-orbit automatic updating method and device based on observation deviation and attenuation calibration parameter, can obtain more stable and accurate instrument observation deviation change and update suggestion, meet the demand of on-orbit change and timely discovery and correction.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of communication technology, in particular to a method for automatically updating on-orbit calibration parameters of a remote sensing instrument in a reflective band. BACKGROUND

[0002] The performance of a remote sensing instrument is affected by the space working environment and the change of device characteristics, and the performance on-orbit changes. The reflective band, especially the short-wave part, has a significant radiation response time decay problem. It is necessary to find the changes in time and provide corresponding calibration correction parameters to compensate for the influence of instrument response decay on radiation data and maintain the quality and stability of radiation data. For instruments with reliable on-board calibration devices, such as MODIS and VIIRS, the relative change information of the instrument radiation response can be obtained in time by using periodic solar observation, and calibration correction can be implemented accordingly. For instruments without reliable on-board calibration devices, the calibration parameter information of a specific time is usually obtained by simultaneous nadir observation (SNO) with high-performance reference instruments and atmospheric top observation (field calibration) based on ground synchronous measurement.

[0003] The technical solution of simultaneous nadir observation (SNO) based on high-performance reference instruments is limited by the orbits of the reference instrument and the target instrument, and the comparison opportunity is limited, such as the cross-comparison samples of two sun-synchronous near-polar orbit satellites distributed in the polar region. The technical solution of atmospheric top observation (field calibration) based on ground synchronous measurement is limited by the implementation and site clear sky conditions of ground synchronous measurement, and there are problems of manpower and material resources cost. Also, the effective comparison opportunity is limited. At the same time, due to the influence of the potential nonlinear response characteristics of the instrument, the spatial and radiation distribution differences of the sample target, the observation geometric conditions, etc., the obtained radiation deviation or calibration parameter has fluctuation, which affects the reliability of the monitoring result. Since it cannot provide time-continuous and stable monitoring results, it cannot meet the demand of timely discovery and reliable correction of the on-orbit change of the instrument response. Therefore, a method for automatically updating on-orbit calibration parameters of a remote sensing instrument in a reflective band is needed. SUMMARY

[0004] In view of the deficiencies of the prior art, the present application provides a method for automatically updating on-orbit calibration parameters of a remote sensing instrument in a reflective band, which solves the problem of not being able to provide time-continuous and stable monitoring results and not being able to meet the demand of timely discovery and reliable correction of the on-orbit change of the instrument response.

[0005] To achieve the above object, the present application is implemented by the following technical scheme:

[0006] A method for automatically updating on-orbit calibration parameters of a remote sensing instrument in a reflective band, comprising the following steps:

[0007] Step one, select stable target sites distributed around the world, extract target instrument observation data, according to the set scheme, based on environmental field data and radiation transfer model, calculate the instantaneous simulated reflectivity of the top of the atmosphere of the stable target site, and obtain the simulated reflectivity factor RefF after distance correction and solar zenith angle correction sim , form the matching sample data pair of instrument observation reflectivity factor RefF mea and simulated reflectivity factor RefF sim , which contains the earth observation count value EV, cold air observation count value SV, satellite zenith angle (SenZ), solar zenith angle (SolZ), and intra-annual day count Julday, etc.

[0008] Step two, according to the set scheme, the matching sample data of observation and simulation reflectivity is subjected to the first quality control screening, and the unreliable factors such as large simulation error, saturated observation data, cloud / dust and other non-clear sky effects are removed, and the matching sample data set after the first quality control is formed.

[0009] Step three, based on the matching sample data set after the first quality control, for the ith day, in the set time window [i-N1, i+N2] days, (EV-SV) and RefF sim data fitting is carried out, and abnormal data is removed for the time window, and the fitting slope Slope of the ith day is obtained, and the first Slope time sequence is formed.

[0010] Step four, based on the Slope time sequence processed in step three, taking the selected date i0 such as launch day and start day as the reference time, the time variation model Func_Slope of Slope is obtained according to the set scheme, and the daily SlopeM time sequence data estimated based on the time variation model Func_Slope is averaged to obtain the average value SlopeA.

[0011] Step five, based on SlopeM, the new observation reflectivity RefF mea2 (RefF mea2 =SlopeM*(EV-SV)) is calculated, and steps two to four are repeated to form the matching sample data set after the second quality control, the second Slope time sequence data set, and the second daily SlopeM time sequence data set.

[0012] Step six, based on the matching sample data set after the second quality control, PDif (PDif=RefF mea / RefF sim -1) and PDif2 (PDif2=RefF mea2 / RefF sim -1) are calculated, and the PDif and PDif2 time sequence data sets are formed.

[0013] Step seven, based on the PDif2 time series data, according to the set scheme, through the multi-year average of daily data processing, the PDif daily correction model Func_PDif is established for the annual change correction of PDif data, and the PDifCor after the annual change correction is obtained;

[0014] Step eight, based on the second Slope time series data set, according to the set scheme, through the detrended and multi-year average of daily data processing, the Slope daily correction model Func_Slope is established for the annual change correction of Slope data, and the SlopeCor after the annual change correction is obtained;

[0015] Step nine, based on the PDifCor and SlopeCor results, according to the set scheme, it is judged whether the scaling coefficient needs to be updated;

[0016] Step ten, when the scaling coefficient needs to be updated in step nine, according to the set scheme, the scaling coefficient to be updated is calculated for the first and non-first update respectively.

[0017] Preferably, the observation reflectivity calculation step in step one is: inputting the target instrument L1 data, searching for the observation pixel point closest to the stable target station according to the latitude and longitude position of the stable target station, extracting the data of n*n window with the observation pixel point as the center, and calculating the window mean (Mean), standard deviation (Std), and variation coefficient CV (CV=Std / Mean), obtaining the observation reflectivity factor RefF according to the radiation scaling conversion relationship, and forming the stable target observation data. mea

[0018] Preferably, the simulation reflectivity calculation step in step one is: inputting the stable target observation data, according to the latitude and longitude and time information, spatiotemporally matching with the pre-constructed atmospheric aerosol / water vapor / ozone, surface BRDF climate state data, and atmospheric temperature and humidity pressure wind environment field data, and using any one of the radiation transfer modes such as MODTRAN and 6SV to calculate the atmospheric top instantaneous reflectivity Ref sim , and performing the solar zenith angle and the distance correction to obtain the simulation reflectivity factor RefF sim .

[0019] Preferably, the specific steps of the quality control screening in step two are:

[0020] S1, according to the set angle quality control threshold, using the solar zenith angle M° and the satellite zenith angle N° to perform SolZ and SenZ identification, and removing the data exceeding the threshold;

[0021] ​S2, according to the set observation saturation threshold, each channel EV identification is carried out, and the data exceeding the threshold is removed;

[0022] S3, according to the set space uniformity quality control reference channel F CV and threshold, CV EV identification is carried out, and the data exceeding the threshold is removed;

[0023] S4, according to the set cloud / sand detection channel F1 and F2 and threshold, channel observation reflectivity feature identification is carried out, and the data exceeding the threshold is removed;

[0024] S5, according to the set observation bias quality control reference channel F RE and threshold, observation and simulation reflectivity relative deviation (RefF mea / RefF sim -1) identification is carried out, and the data exceeding the threshold is removed.

[0025] Preferably, the specific steps in step three for each time window are as follows:

[0026] S1, remove the data less than 0 (EV-SV) or RefF sim ;

[0027] S2, according to the set model, fit (EV-SV) as x and RefF sim as y;

[0028] S3, calculate the standard deviation of the deviation between the simulation reflectivity and the model estimated reflectivity, and remove the data exceeding the threshold with the threshold being 2 times the standard deviation;

[0029] S4, repeat S2, S3 and S2 to obtain the fitting slope Slope of the i-th day.

[0030] Preferably, the specific steps of processing the Slope time sequence in step three in step four are as follows:

[0031] S1, taking the selected date i0 as the reference time, based on the Slope time sequence, according to the set model such as linear, quadratic polynomial, etc., time change trend fitting is carried out;

[0032] S2, calculate the standard deviation of the deviation between the actual Slope value and the model estimated value, and remove the data exceeding the threshold with the threshold being 2 times the standard deviation;

[0033] S3, repeat steps S1, S2 and S1 to obtain the time change model Func_Slope (DSL) of Slope, wherein DSL is the day number of distance from the reference time;

[0034] S4, based on the Func_Slope(DSL) model estimation, obtain the daily SlopeM time series;

[0035] S5, average the daily SlopeM time series data to obtain SlopeA.

[0036] Preferably, the specific steps of step seven for establishing the PDif annual correction model and performing change correction are:

[0037] S1, based on the PDif2 time series data, average the PDif2 data of the same year within the day count Julday to obtain the multi-year average daily PDif2 within the year ave-daily data;

[0038] S2, calculate the mean of the PDif2 ave-daily data set to obtain PDif2 ave ;

[0039] S3, according to a set model such as a Sin function, perform time variation fitting of PDif2 ave-daily -PDif2 ave to obtain the PDif annual daily correction model Func_PDif(Julday), Julday = 1, 2... 366;

[0040] S4, obtain the observation bias PDifCor = PDif - Func_PDif(Julday) after correction.

[0041] Preferably, the specific steps of step eight for establishing the Slope annual correction model and performing change correction are:

[0042] S1, based on the second Slope and the daily SlopeM time series data, remove the trend of the data set over time, Slope2 = Slope / SlopeM;

[0043] S2, average the Slope2 data of the same year within the day count Julday to obtain the multi-year average daily Slope2 ave-daily data within the year;

[0044] S3, calculate the mean of the Slope2 ave-daily data set, Slope2 ave ;

[0045] S4, according to a set model such as a Sin function, perform Slope2 ave / Slope2 ave-dailySlopeCor = Slope * Func_Slope(Julday) of the SlopeCor time series data of the i day, and the SlopeCor time series data of the i day is obtained by the following steps:

[0046] SlopeCor = Slope * Func_Slope(Julday) of the SlopeCor time series data of the i day, and the SlopeCor time series data of the i day is obtained by the following steps:

[0047] The specific steps in the ninth step are preferably determined as follows:

[0048] S1, for the i day, the PDifCor time series data is processed by the moving average method in the set time window, and the PDifCor time series data of the i day is obtained by the following steps: i ;

[0049] S2, in the time sequence, the PDifCor time series data is processed by the moving average method in the set time window, and the PDifCor time series data of the i day is obtained by the following steps: i When the PDifCor time series data exceeds the set threshold value, it is determined that the calibration coefficient needs to be updated.

[0050] S3, for the non-first update, the decay rate scheme is used for judgment.

[0051] S4, for the i day, the decay rate Decay is calculated relative to the last update: Decay = 1 - Func_NRes(i) / Func_NRes(T1), where T1 is the date of the first update:

[0052] The selected date i0 (such as the launch date, the start-up date) is taken as the reference time, and the normalized response NRes is calculated: the SlopeCor results from the reference time to the i day are extracted, and the time variation trend fitting is performed according to the set model, such as linear, quadratic polynomial, etc., to obtain the time variation model Func_SlopeCor(DSL) of SlopeCor; NRes(i) = Func_SlopeCor(i0) / Func_SlopeCor(i);

[0053] Based on the NRes time series, the time variation trend fitting is performed according to the set model, such as linear, quadratic polynomial, etc., to obtain the time variation model Func_NRes(DSL) of NRes;

[0054] The decay rate Decay is calculated: Decay = 1 - Func_NRes(i) / Func_NRes(T1), where T1 is the date of the first update:

[0055] S5, when the Decay exceeds the set threshold value, it is determined that the calibration coefficient needs to be updated.

[0056] For the first update, the Func_SlopeCor(T1) is taken as the calibration coefficient to be updated, and for the non-first update, the Func_SlopeCor(T1) * (1 + Decay) is taken as the calibration coefficient to be updated.

[0057] The application provides a method for automatically updating a remote sensing instrument reflection band calibration parameter on orbit.

[0058] 1. The application uses a multi-stable target to simulate reflectivity to calculate the deviation of observed reflectivity, and through two-layer data quality control strategies, the reliability of the processing result can be improved, and relatively accurate observed reflectivity deviation change conditions can be obtained; through normalization response calculation, the instrument performance change conditions can be relatively accurately reflected; in addition, the instrument observation deviation and the normalization response both consider the annual periodic change, and a corresponding annual daily correction model is established, so that the influence of the annual periodic change can be removed, and more stable and accurate long-term change results can be obtained.

[0059] 2. In the application, whether to need calibration update is judged by using two schemes of reflectivity relative deviation and decay rate, the relative deviation scheme is used for judgment first, the update judgment based on the relative deviation uses a time window moving average scheme, the influence of data randomness can be removed; after implementing a calibration update, the decay rate scheme is used for judgment, the influence of data randomness of the relative deviation scheme can be further avoided, and the judgment result is more stable.

[0060] 3. Based on the method, the application can realize daily real-time monitoring of instrument observation deviation and performance change, timely provide business calibration update judgment, support calibration correction compensation, avoid data quality deterioration, and keep the data precision stable in a predetermined range. BRIEF DESCRIPTION OF DRAWINGS

[0061] Figure 1 is a flow processing diagram of the application;

[0062] Figure 2 is a PDif annual correction model result graph of the application taking band 1 of FY-3D / MERSI as an experimental example;

[0063] Figure 3 is a calibration coefficient update process PDif result graph of the application taking band 1 of FY-3D / MERSI as an experimental example. DETAILED DESCRIPTION

[0064] The technical solutions in the embodiments of the application will be clearly and completely described below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, rather than all the embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the protection scope of the application.

[0065] Embodiment:

[0066] Please refer to the drawingsFigure 1 The embodiment of the present application provides a kind of remote sensing instrument reflection band calibration parameter on-orbit automatic updating method, comprising the following steps:

[0067] Step one, select the stable target site distributed in the world, extract target instrument observation data, according to the set scheme, based on environmental field data and radiation transfer model, the instantaneous simulation reflectivity of the atmosphere top of stable target site is calculated, and the matching sample data pair of instrument observation and simulation reflectivity is formed, wherein it contains earth observation count value EV, cold observation count value SV, satellite zenith angle (SenZ), solar zenith angle (SolZ), day count Julday in a year, etc.;

[0068] Wherein, the observation reflectivity calculation step is: input target instrument L1 data, according to the latitude and longitude position of stable target site, find the observation pixel point closest to it, extract the data of n*n window with it as center, and calculate window mean (Mean), standard deviation (Std), variation coefficient CV (CV=Std / Mean), form stable target observation data, including earth observation count value (EV), cold observation count value (SV), earth observation count value after dark background deduction (DN, DN=EV-SV), observation reflectivity of solar located at zenith average day-earth distance (RefF mea ), satellite zenith angle (SenZ), solar zenith angle (SolZ);N takes 3, (k0, k1, k2) is calibration coefficient.

[0069] The formula is: RefF mea =k0+DN*k1+DN*DN*k2

[0070] Wherein, the simulation reflectivity calculation step is: input stable target observation data, according to latitude and longitude and time information, with the atmospheric aerosol / water vapor / ozone, ground BRDF climatological data constructed in advance and atmospheric temperature and humidity pressure wind and other environmental field data are matched in space and time, and for target observation angle, any one of radiation transfer mode such as MODTRAN, 6SV is used to calculate atmospheric top instantaneous reflectivity Ref sim , and solar zenith angle and day-earth distance correction is carried out, and solar located at zenith average day-earth distance simulation reflectivity factor (RefF sim =Ref sim *cos (SolZ)*dcorF) is obtained.

[0071] Day-earth distance correction factor formula is: dcorF=1. / (1-0.01673*cos (0.9856*(Julday-4)*3.1415926 / 180)) 2

[0072] Julday is the day count within the year relative to January 1

[0073] Step two, the matching sample data of observed and simulated reflectivity is screened according to a set scheme for the first time, and unreliable factors such as large simulation error, saturated observation data, cloud / dust and other non-clear sky effects are removed to form a matching sample data set after the first quality control.

[0074] S1, according to the set angle quality control threshold, such as the solar zenith angle 60° and the satellite zenith angle 50°, SolZ and SenZ are identified, and data exceeding the threshold are removed;

[0075] S2, according to the set observation saturation threshold, such as 4000, EV of each channel is identified, and data exceeding the threshold are removed;

[0076] S3, according to the set spatial uniformity quality control reference channel and threshold, such as the second channel (550 nm) and 0.1, CV EV is identified, and data exceeding the threshold are removed;

[0077] S4, according to the set cloud / dust detection reference channel and threshold, such as the first channel (470 nm) and the second channel (550 nm) and 0.65, the sum of the observation reflectivity factors of the two channels is identified, and data exceeding the threshold are removed;

[0078] S5, according to the set observation deviation quality control reference channel and threshold, such as the fourth channel (865 nm) and 30%, the relative deviation (RefF mea / RefF sim -1) of observed and simulated reflectivity is identified, and data exceeding the threshold are removed;

[0079] Step three, based on the matching sample data set after the first quality control, for the ith day, a set time window is used, such as [i-4, i+5] days, (EV-SV) and RefF sim data fitting is performed, and the fitting slope forms the first Slope time series, wherein, for each time window:

[0080] S1, remove data less than 0 of (EV-SV) or RefF sim ;

[0081] S2, taking (EV-SV) as x and RefF sim as y, fitting is performed according to a set model, such as linear;

[0082] S3, the standard deviation of the deviation between the simulated reflectivity and the model estimated reflectivity is calculated, and data exceeding the threshold of 2 times the standard deviation are removed;

[0083] S4, repeat S2, S3 and S2, obtain the i-th day of the fitting slope Slope.

[0084] Step four, based on the Slope time series processed in step three, obtain the daily SlopeM time series and the average SlopeA, the specific steps are as follows:

[0085] S1, taking the selected date i0 such as launch day, start day, etc. as the reference time, based on the Slope time series, according to the set model such as linear, quadratic polynomial, etc., the time change trend fitting is carried out;

[0086] S2, calculate the standard deviation of the deviation between the actual Slope value and the model estimated value, and remove the data exceeding the threshold value of 2 times the standard deviation;

[0087] S3, repeat steps S1, S2 and S1, obtain the time change model Func_Slope(DSL) of Slope, wherein DSL is the day number of distance from the reference time;

[0088] S4, based on the Func_Slope(DSL) model estimation, obtain the daily SlopeM time series;

[0089] S5, average the daily SlopeM time series data to obtain SlopeA.

[0090] Step five, based on SlopeM, calculate the new observation reflectivity RefF mea2 (RefF mea2 =SlopeM*(EV-SV)), repeat steps two to four to form the second time quality controlled matching sample data set, the second time Slope time series data set, and the second time daily SlopeM time series data set.

[0091] Step six, based on the second time quality controlled matching sample data set, calculate PDif (PDif = RefF mea / RefF sim -1) and PDif2 (PDif2 = RefF mea2 / RefF sim -1), form the PDif and PDif2 time series data set.

[0092] Step seven, based on the PDif2 time series data, establish the PDif intra-year correction model for the intra-year change correction of PDif data, the detailed steps are as follows:

[0093] S1, based on the PDif2 time series data, average the PDif2 data of the same intra-year day count Julday to obtain the intra-year daily PDif2ave-daily data;

[0094] S2, calculate the mean of PDif2 data set, obtain PDif2 ave-daily data set, obtain PDif2 ave ;

[0095] S3, according to the set model, such as Sin function, PDif2 ave-daily -PDif2 ave time change fitting, obtain the daily correction model of PDif in the year Func_PDif(Julday), Julday=1, 2…366;

[0096] S4, obtain the observation bias after correction PDifCor=PDif-Func_PDif(Julday).

[0097] Step eight, based on the second Slope time series data set, establish the Slope annual correction model for the annual change correction of Slope data, the detailed steps are as follows:

[0098] S1, based on the second Slope and daily SlopeM time series data, remove the trend of data set changing with time, Slope2=Slope / SlopeM;

[0099] S2, average the Slope2 data of the same day count Julday in the year, obtain the annual daily Slope2 ave-daily data;

[0100] S3, calculate the mean of Slope2 ave-daily data set, Slope2 ave ;

[0101] S4, according to the set model, such as Sin function, Slope2 ave / Slope2 ave-daily time change fitting, obtain the daily correction model of Slope in the year Func_Slope(Julday), Julday=1, 2…366;

[0102] S5, obtain the calibrated slope SlopeCor=Slope*Func_Slope(Julday) after correction.

[0103] Step nine, based on the results of PDifCor and SlopeCor, according to the set scheme, judge whether the calibration coefficient needs to be updated, the detailed steps are as follows:

[0104] S1, for the i-th day, a moving average processing is performed on the PDifCor time series data in a set time window (e.g., [i-5, i] days) to obtain PDif i ;

[0105] S2, in time sequence, PDif i exceeds a set threshold value (e.g., 0.02), it is determined that the scaling coefficient needs to be updated;

[0106] S3, for non-first-time updates, a decay rate scheme is used for judgment;

[0107] S4, for the i-th day, a decay rate Decay relative to the last update is calculated:

[0108] A selected date i0(e.g., launch day, start-up day) is taken as a reference time, and a normalized response NRes is calculated: the SlopeCor results from the reference time to the i-th day are extracted, a time change trend fitting is performed according to a set model, such as linear, quadratic polynomial, etc., to obtain a time change model Func_SlopeCor(DSL) of SlopeCor; NRes(i) = Func_SlopeCor(i0) / Func_SlopeCor(i);

[0109] Based on the NRes time series, a time change trend fitting is performed according to a set model, such as linear, quadratic polynomial, etc., to obtain a time change model Func_NRes(DSL) of NRes;

[0110] A decay rate Decay is calculated: Decay = 1-Func_NRes(i) / Func_NRes(T1), where T1 is the date of the first update.

[0111] S5, when Decay exceeds a set threshold value (e.g., -0.02), it is determined that the scaling coefficient needs to be updated.

[0112] Step ten, it is determined that the scaling coefficient needs to be updated, for the first-time update, Func_SlopeCor(T1) is taken as the scaling coefficient to be updated; for non-first-time updates, Func_SlopeCor(T1)*(1+Decay) is taken as the scaling coefficient to be updated.

[0113] The present application is based on the simulated reflectivity of multi-stable targets, uses the deviation between the simulated reflectivity and the observed reflectivity, and carries out the annual cycle (daily) change correction of the deviation, which can continuously and accurately reflect the change of the instrument observation deviation; through the instrument normalization response and decay rate based on the calibration slope, and the annual cycle (daily) change correction, more stable and accurate instrument change results are obtained. When judging whether to update the calibration, the present application uses two schemes of reflectivity relative deviation and decay rate, first uses the time window moving average method to reduce the influence of data randomness, and then further improves the stability of the judgment through the decay rate scheme. This method realizes the daily real-time monitoring of the instrument observation deviation and performance change, provides timely judgment of calibration update, supports calibration correction compensation, ensures the data quality does not deteriorate, and maintains the data precision within the predetermined range.

[0114] Experimental example:

[0115] Please refer to the attached Figure 2 After the data is processed by the Func_PDif correction model, the volatility is obviously reduced, please refer to the attached Figure 3 There are three times of calibration coefficient update in the whole process, which are July 18, 2019, November 18, 2020 and April 25, 2022.

[0116] Although the embodiments of the present application have been shown and described, it can be understood by those skilled in the art that various changes, modifications, replacements and variations can be made to these embodiments without departing from the principles and spirits of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for on-orbit automatic updating of reflectance band calibration parameters of a remote sensing instrument, characterized in that, The method comprises the following steps: Step one, select stable target sites distributed around the world, extract target instrument observation data, according to the set scheme, based on environmental field data and radiation transfer model to calculate the atmospheric top instantaneous simulated reflectivity of stable target sites, and through the distance between the earth and the sun and the solar zenith angle correction to get the simulated reflectivity factor RefF sim , form the matching sample data pair of instrument observation reflectivity factor RefF mea and simulated reflectivity factor RefF sim , which contains the earth observation count value EV, cold air observation count value SV, satellite zenith angle SenZ, solar zenith angle SolZ and annual day count Julday; Step two, the matching sample data of observation and simulated reflectivity is subjected to first quality control screening according to a set scheme, and unreliable factors including large simulation error, saturated observation data and non-clear sky influence of cloud / dust are removed to form a first quality controlled matching sample data set; Step three, based on the matched sample data set after the first quality control, for the ith day, with the set time window [i-N1, i+N2] days, EV-SV and RefF sim Data fitting, and abnormal data elimination for the time window, obtain the fitting slope Slope of the ith day, the fitting slope forms the first Slope time series; Step four, based on the slope time sequence processed in step three, a time variation model Func_Slope of slope is obtained according to a set scheme with the launch date and the selected date i0 as the reference time, and the daily slopeM time sequence data estimated based on the time variation model Func_Slope is averaged to obtain the average value SlopeA; Step five, calculate new observed reflectivity RefF based on SlopeM mea2 =SlopeM*(EV-SV), repeat steps two-four to form a second quality controlled matched sample dataset, a second Slope time series dataset, a second daily SlopeM time series dataset; Step six, based on the matched sample dataset after the second quality control, calculate PDif = RefF mea / RefF sim - 1 and PDif2 = RefF mea2 / RefF sim - 1, forming the PDif and PDif2 time series dataset; Step seven, based on the PDif2 time sequence data, a PDif daily subscription model Func_PDif is established through the calculation and processing of the multi-year average daily data within a year according to a set scheme, which is used for the intra-annual variation correction of PDif data to obtain the PDifCor after the intra-annual variation correction; Step eight, based on the second slope time sequence data, a slope daily subscription model Func_Slope is established through the calculation and processing of the multi-year average daily data within a year according to a set scheme, which is used for the intra-annual variation correction of slope data to obtain the SlopeCor after the intra-annual variation correction; Step nine, based on the PDifCor and SlopeCor results, it is judged according to a set scheme whether the scaling coefficient needs to be updated; Step ten, when the scaling coefficient needs to be updated according to the judgment in step nine, the scaling coefficient to be updated is calculated according to a set scheme for the first update and non-first update respectively.

2. The method according to claim 1, wherein the method is characterized by, The observation reflectivity calculation step in the step one is: input target instrument L1 data, find the observation pixel point closest to the stable target station according to the latitude and longitude position of the stable target station, extract n*n window data with the observation pixel point as the center, and calculate the window mean Mean, standard deviation Std, and variation coefficient CV=Std / Mean, and obtain the observation reflectivity factor RefF according to the radiation scaling conversion relationship mea , to form stable target observation data.

3. The method according to claim 1, wherein the method is characterized by: The simulation reflectivity calculation step in step one is: input stable target observation data, according to the latitude, longitude and time information, time and space matching is performed with the atmospheric aerosol / water vapor / ozone, surface BRDF climatic state data and environmental field data including atmospheric temperature, humidity, pressure and wind which are constructed in advance, and for the target observation angle, the instantaneous reflectivity Ref of the top of the atmosphere is calculated by using any one of the radiation transfer modes in MODTRAN and 6SV sim , and the solar zenith angle and the distance between the sun and the earth are corrected to obtain the simulation reflectivity factor RefF of the sun located at the average distance between the sun and the earth sim .

4. The method according to claim 1, wherein, The specific steps of the quality control screening in step two are: S1, according to the set angle quality control threshold, the SolZ and SenZ are identified by using the solar zenith angle M° and the satellite zenith angle N°, and the data exceeding the threshold is removed; S2, according to the set observation saturation threshold, the EV of each channel is identified, and the data exceeding the threshold is removed; S3. According to the set spatial uniformity quality control reference channel F CV and threshold, CV EV discrimination, eliminate data exceeding threshold; S4, according to the set cloud / dust detection channel F1 and F2 and the threshold, the channel observation reflectivity feature is identified, and the data exceeding the threshold is removed; S5. According to the set observation bias quality control reference channel F RE and threshold, the observation and simulated reflectivity relative bias RefF mea / RefF sim -1 identify, reject data exceeding threshold.

5. The method according to claim 1, wherein, The specific steps for each time window in step three are: S1, reject EV-SV vs. RefF sim any one of the data is less than 0; S2, with EV-SV as x, RefF sim as y, according to the set model S3, the standard deviation of the deviation between the simulated reflectivity and the model estimated reflectivity is calculated, and the data exceeding the threshold of 2 times the standard deviation is removed; S4, S2, S3 and S2 are repeated to obtain the fitting slope Slope of the i-th day.

6. The method according to claim 1, wherein, The specific steps of processing the slope time sequence in step three in step four are: S1, based on the slope time sequence, the time variation trend is fitted according to any one of the set models of linear and quadratic polynomial with the selected date i0 as the reference time; S2, the standard deviation of the deviation between the actual slope value and the model estimated value is calculated, and the data exceeding the threshold of 2 times the standard deviation is removed; S3, steps S1, S2 and S1 are repeated to obtain the time variation model Func_Slope (DSL) of slope, wherein DSL is the daily record number from the reference time. S4, based on the Func_Slope(DSL) model estimation, obtain the daily SlopeM time series; S5, average the daily SlopeM time series data to obtain SlopeA.

7. The method according to claim 1, wherein the method is characterized by: The specific steps of step seven are to establish a PDif annual correction model and to correct the changes, and are as follows: S1, based on the PDif2 time series data, average the PDif2 data of the same year within the day count Julday to obtain the multi-year average daily PDif2 within the year ave-daily data; S2, compute PDif2 ave-daily mean of dataset, obtain PDif2 ave ; S3. Obtain PDif2 from the fitted model ave-daily -PDif2 ave time variation fitting, obtain the daily corrected model Func_PDif(Julday), Julday = 1,2...366; S4, obtain the observation bias PDifCor after correction PDif-Func_PDif(Julday).

8. The method according to claim 1, wherein, The specific steps of step eight are to establish a Slope annual correction model and to correct the changes, and are as follows: S1, based on the second Slope and the daily SlopeM time series data, remove the trend of the data set changing with time, Slope2=Slope / SlopeM; S2, average the Slope2 data for the same year's day count, Julday, to obtain a multi-year average of the daily Slope2 within a year ave-daily data; S3, compute Slope2 ave-daily Mean of dataset, Slope2 ave ; S4, according to the set model, Slope2 ave Slope2 ave-daily time variation fitting, obtain Slope daily correction model Func_Slope(Julday), Julday=1, 2…366; S5, obtain the calibrated slope SlopeCor after correction Slope*Func_Slope(Julday).

9. The method according to claim 1, wherein, The specific steps in step nine are as follows: S1, for the ith day, with a set time window, the PDifCor time series data is processed by moving average, and the PDif i ; S2, in time sequence, PDif i If the value exceeds a set threshold, it is determined that the scaling factor update is needed. S3, for non-first-time updates, use the decay rate scheme for judgment; S4, for the i-th day, calculate the decay rate relative to the last update Decay Take the selected date i0 as the reference time, calculate the normalized response NRes: extract the SlopeCor results from the reference time to the i-th day, set the model according to any one of the linear and quadratic polynomials, and perform time trend fitting to obtain the time variation model Func_SlopeCor(DSL) of SlopeCor; NRes(i)=Func_SlopeCor(i0) / Func_SlopeCor(i); Based on the NRes time series, set the model according to any one of the linear and quadratic polynomials, and perform time trend fitting to obtain the time variation model Func_NRes(DSL) of NRes; Calculate the decay rate Decay: Decay=1-Func_NRes(i) / Func_NRes(T1), where T1 is the date of the first update: S5, when Decay exceeds the set threshold, it is determined that the calibration coefficient needs to be updated.

10. The method according to claim 9, wherein the method is characterized by: In step ten, for the first update, take Func_SlopeCor(T1) as the updated calibration coefficient; for non-first-time updates, take Func_SlopeCor(T1)*(1+Decay) as the updated calibration coefficient.

Images