Visible light wave band spectral response function inversion method

Abstract

The invention discloses a visible light wave band spectral response function inversion method which comprises the following steps: acquiring a satellite on-orbit spectral response function and atmospheric top hyperspectral reflectivity simulation data through data simulation, and calculating satellite on-orbit reflection data through a forward model; determining a target function of the inversion model according to the satellite in-orbit reflection data; determining constraint conditions of the inversion model; training the inversion model through a linear programming solver to obtain a trained inversion model; and inputting the atmospheric top hyperspectral reflectance data and the satellite reflectance value into the trained inversion model to obtain a corresponding spectral response function. According to the satellite visible light wave band and the satellite characteristics, an inversion algorithm based on linear programming is used, and sufficient prior information of a target satellite is fully utilized to formulate reasonable constraints; and a visible light wave band spectral response function inversion method with universality and robustness is obtained.

Description

Technical Field

[0001] The present application relates to the field of optical remote sensing science and technology, and in particular to a method for inverting a spectral response function in a visible light band. Background Art

[0002] The spectral response function is an important parameter of the remote sensor imaging system, also known as the bandpass function, slit function, etc. It characterizes the response characteristics of the sensor to radiation in different wavelength ranges. It is closely related to the remote sensor radiation calibration coefficient and directly affects the observed target radiation.

[0003] At present, radiometric calibration has made considerable progress, and the accuracy of spectral calibration has become increasingly important. It can be said that the accuracy of on-orbit spectral calibration greatly affects the quality of remote sensing inversion products. When an inaccurate instrument spectral response function is used, the wavelength scale may differ from the true scale by a few hundredths of a nanometer. This uncertainty in the wavelength scale may introduce pseudo-noise, which in turn affects algorithms based on spectral information, thereby affecting the accuracy of the acquired data.

[0004] Most spectroradiometers cannot be calibrated on-orbit due to the lack of onboard spectral calibration equipment. Typically, these instruments use the pre-launch spectral response function, which is the spectral response function obtained by laboratory spectral calibration. It is also assumed that the performance of the spectral response function will not change throughout the service life. In addition, the degradation of the spectral response function has been found in the monitoring and flight calibration of many optical remote sensing equipment. This phenomenon is likely caused by certain factors such as instrument aging, chemical molecular contamination in the space environment, and temperature changes in satellite instruments, as shown by a large number of studies.

[0005] In summary, due to the limitations of laboratory equipment and environment before launch, the spectral response function obtained by laboratory spectral calibration may be inaccurate. At the same time, the spectral response function may degrade during on-orbit operation. These problems make the spectral response function unreliable and cannot be directly applied. Summary of the invention

[0006] In view of the above shortcomings in the prior art, the present invention provides a visible light band spectral response function inversion method to solve the problem that the spectral response function obtained by laboratory spectral calibration may be inaccurate due to the limitations of laboratory equipment and environment before launch; at the same time, the spectral response function may be degraded during orbital operation. These problems lead to the problem that the spectral response function is unreliable and cannot be directly applied.

[0007] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is: a visible light band spectral response function inversion method, comprising:

[0008] Through data simulation, the satellite on-orbit spectral response function and the simulated data of the high spectral reflectance at the top of the atmosphere are obtained, and the satellite reflectance data is calculated through the forward model;

[0009] Determine the objective function of the inversion model based on satellite reflection data;

[0010] Determine the constraints of the inversion model;

[0011] The inversion model is trained by a linear programming solver to obtain a trained inversion model;

[0012] The top-of-atmosphere hyperspectral reflectance data and satellite reflectance values ​​are input into the trained inversion model to obtain the corresponding spectral response function.

[0013] Furthermore: the satellite on-orbit spectral response function is obtained through degradation simulation, and the high spectral reflectance simulation data of the top of the atmosphere is obtained through libRadtran radiation transfer model simulation.

[0014] Further: The expression for calculating satellite reflection data is:

[0015]

[0016] Among them, refm is the satellite observed reflectivity, S is the high spectral reflectivity of the top of the atmosphere, rsr represents the satellite on-orbit spectral response function, and λ represents the band.

[0017] Further: The objective function Z of the inversion model is:

[0018]

[0019]

[0020] Among them, δ i Represents the absolute value of the residual reflectance of the i-th satellite. The value range of i is 1 to m. m represents the number of groups of input atmospheric top hyperspectral reflectance data. i represents the reflectivity of the i-th band, rsr λ represents the spectral response function obtained by inversion, rsr λ 0 represents the spectral response function measured in the laboratory before the satellite was launched, and v represents the energy scaling factor.

[0021] Furthermore: the constraints of the inversion model include non-negativity constraints, normalization constraints, variation constraints and difference constraints.

[0022] Further: The expression of the constraint condition of the inversion model is:

[0023]

[0024] Among them, α represents the difference factor, rsr λ(j) represents the jth sampling point of the spectral response function, rsr λ(k) represents the kth sampling point of the spectral response function, is the variation size, and n means that the spectral response function has n sampling points.

[0025] The beneficial effects of the present invention are:

[0026] For the satellite visible light band, according to the characteristics of the satellite, for the hyperspectral reflectance data of the top of the atmosphere simulated by observations of multiple ground objects, an inversion algorithm based on linear programming is used to make full use of sufficient prior information of the target satellite to formulate reasonable constraints, which has good universality and robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 This is a flow chart of the inversion method of the spectral response function in the visible light band.

[0028] Figure 2 This is the distribution diagram of FY3A MERSI spectral response function.

[0029] Figure 3 This is the distribution diagram of FY3A VIRR spectral response function.

[0030] Figure 4 This is the distribution diagram of FY3B MERSI spectral response function.

[0031] Figure 5 This is the distribution diagram of FY3B VIRR spectral response function.

[0032] Figure 6 This is the distribution diagram of FY3C MERSI spectral response function.

[0033] Figure 7 This is the distribution diagram of FY3C VIRR spectral response function.

[0034] Figure 8 This is the distribution diagram of AQUA MODIS spectral response function.

[0035] Fig. 9 Schematic diagram of the degradation of the spectral response function of each channel.

[0036] Fig.10 It is a custom Sigmoid function distribution graph.

[0037] Fig.11 This is a graph showing the variation of the spectral response function.

[0038] Fig.12 This is a graph showing the relationship between adjacent points under mutation constraints.

[0039] Fig.13 This is a diagram showing the relationship between adjacent points under the combined effect of variation constraints and difference constraints.

[0040] Fig.14 This is the inversion result diagram of the accuracy experiment.

[0041] Fig.15 This is the inversion result diagram of the universal experiment.

[0042] Fig.16 The overall flow chart for the bias analysis of top-of-atmosphere hyperspectral reflectance simulation.

[0043] Fig.17 Correlation distribution plots of each channel in the three features selected for the robustness experiment.

[0044] Fig.18 This is the relative deviation mean of each channel in the robustness experiment.

[0045] Fig.19 This is the relative deviation standard deviation diagram of each channel in the robustness experiment.

[0046] Fig. 20 This is the inversion result diagram of the robust universality experiment. DETAILED DESCRIPTION

[0047] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0048] like Figure 1 As shown, in one embodiment of the present invention, a method for inverting a spectral response function in a visible light band is provided, comprising:

[0049] Through data simulation, the satellite's on-orbit spectral response function and the simulated data of the top-of-atmosphere high spectral reflectance are obtained, and the satellite's on-orbit reflectance data is calculated through the forward model;

[0050] Determine the objective function of the inversion model based on the satellite's on-orbit reflection data;

[0051] Determine the constraints of the inversion model;

[0052] The inversion model is trained by a linear programming solver to obtain a trained inversion model;

[0053] The top-of-atmosphere hyperspectral reflectance data and satellite reflectance values ​​are input into the trained inversion model to obtain the corresponding spectral response function.

[0054] The essence of satellite remote sensing is inversion. The inversion problem exists relative to the forward problem. In the inversion of remote sensing physics, the forward problem is always considered to be known. This embodiment starts with the forward problem, first describes its basic principles, and then describes the inversion modeling method and technology in detail. The expression for calculating the satellite's on-orbit reflection data is as follows:

[0055]

[0056] Among them, refm is the satellite observed reflectivity, S is the high spectral reflectivity of the top of the atmosphere, rsr represents the satellite on-orbit spectral response function, and λ represents the band.

[0057] The inversion model infers the spectral response function of the on-orbit instrument based on the given top-of-atmosphere hyperspectral reflectance data and the given satellite reflectance value. Specifically, the inversion model uses the top-of-atmosphere hyperspectral reflectance data and the satellite reflectance value as input data to achieve the inversion of the on-orbit spectral response function. However, the focus of this study is to first solve the inversion algorithm problem. In this process, the satellite reflectance value used is calculated by the forward model, and is generated by the hyperspectral reflectance data and the on-orbit spectral response function through integral simulation. Therefore, the data required for the inversion experiment include the on-orbit spectral response function and the top-of-atmosphere hyperspectral reflectance data. For this reason, data preparation is required before the inversion experiment is performed in this embodiment. In the data preparation stage, we will obtain the on-orbit spectral response function through degradation simulation, and use the llibRadtran radiation transfer model to simulate and obtain the top-of-atmosphere hyperspectral reflectance data. In order to verify the robustness of the inversion, the study will analyze the simulation deviation when simulating the top-of-atmosphere hyperspectral reflectance, and take it into account in the inversion process of the spectral response function.

[0058] Preferably, if Figure 2-Figure 8 As shown, in this embodiment, the visible light bands respectively show the spectral response function distribution of the visible light bands of the FY3 A / B / C MERSI / VIRR sensors and the AQUA MODIS sensor, which are channels 1, 2, 3, 8, 9, 10, 11, 12, and 13 of FY3 A / B / C MERSI, channels 1, 7, 8, and 9 of FY3 A / B / C VIRR, and channels 1, 3, 4, 8, 9, 10, 11, and 12 of AQUA MODIS.

[0059] pass Figure 2-Figure 8 It can also be seen that the common shapes of the spectral response function are single-peak shape and double-peak shape, and the storage form of the spectral response function contains two columns, namely wavelength (in nm) and response value (ranging from 0 to 1). The wavelength is in ascending order, and the value range is selected from 390nm to 800nm ​​in the visible light range.

[0060] The spectral response function may degrade in many forms, such as central wavelength drift, response value change, bandwidth scaling, etc. Among them, the change of response value is the most important form of degradation, and other forms of change can be regarded as a special case of responsivity change. In this embodiment, within the visible light band, spectral response functions are selected for experiments based on the characteristics of wide bandwidth, narrow bandwidth, single peak and double peak, and the degradation of the spectral response function on orbit is simulated by adjusting the response value. The selected satellite and its sensor information are shown in Table 1 below:

[0061] Table 1 Information of selected satellites and their sensors

[0062]

[0063]

[0064] The degradation of the spectral response function of each channel is as follows: Fig. 9 As shown, Fig. 9 Here, ori_rsr represents the spectral response function measured before launch, deg_rsr represents the spectral response function after simulated degradation, and the naming rule of the spectral response function is "satellite name_sensor name_channel number". For example, the spectral response function corresponding to the 12th channel of the AUQA satellite MODIS sensor is AUQA_MODIS_ch12.

[0065] Preferably, the simulation process of the top-of-atmosphere hyperspectral reflectance can also be simulated by the uvspec model. The basic use process of uvspec is divided into three parts:

[0066] A1. Prepare input files;

[0067] A2, uvspec model call;

[0068] A3. Output file reading.

[0069] The core part is the preparation of the input file. uvspec will read the data of each parameter according to the input file. Therefore, when preparing the input file, you must also prepare and process the input data corresponding to each parameter. In addition, after processing the output file of uvspec, the radiance value is obtained, and the simulated atmospheric top hyperspectral reflectance result data can be obtained through subsequent calculation.

[0070] In step A1, the input of the model includes the following parts, such as the composition setting of the atmosphere. The absorption and scattering characteristics of these components have two sources, one is the data provided by the llibRadtran software package, and the other is the custom method, that is, the user of the transmission model provides it; the boundary conditions can set the reflection surfaces at the top and bottom of the atmosphere; the radiation transfer equation solver and other solution-related settings, uvspec provides a variety of radiation transfer equation solvers; the user specifies the sun-earth distance correction parameters and other parameters for post-processing tasks. The uvspec input file contains multiple lines of command information, each line represents a uvspec input. Each line contains a command name, followed by the corresponding parameter value. Information annotations can be used with "#". The following example is an input file set with the ocean as the ground object target:

[0071]

[0072]

[0073] In order to ensure that the uvspec model can run correctly and calculate accurately, the model input data needs to be preprocessed before use to meet its requirements and specifications for the input data format, ensuring that there will be no simulation failures and errors due to data not meeting the specifications.

[0074] The processing of input data mainly involves extracting, converting and saving data according to the requirements of the llibRadtran manual. Specifically, it is necessary to extract atmospheric parameters, surface parameters, spectral data, etc., and perform necessary conversions on them. Finally, the processed data is named and stored according to the requirements of the model specification so that it can be used as input data when using the uvspec model. The data requirements have been introduced in the previous section. Many specified data can directly use the data provided by the software package. Only the data provided by the user needs to be processed. Among them, the solar zenith angle, total ozone column (universal unit DU), altitude, sun-earth distance (days) and 550nm aerosol optical thickness can be extracted from the relevant product data according to the time and longitude and latitude positions, and no additional processing is required. The remaining data that require additional processing are shown in Table 2 below:

[0075] Table 2 Description of data requiring additional processing

[0076]

[0077]

[0078] In step A2, the uvspec command is used in the terminal to call the uvspec radiation transfer model. The prepared input file address and the preset output file address are provided to the command to run the program. After the operation is completed, the target data can be found in the output file.

[0079] In step A3, extracting the radiance value in the output file, calculating the hyperspectral reflectance and the final storage of the data can be summarized as output file processing.

[0080] There are two factors that affect the output file format: the selection of the uvspec radiation transfer equation solver and whether the output azimuth (phi) and the cosine value of the output polar angle (umu) are defined. In the example, the disort solver is selected and phi and umu are set. Therefore, the output file format is as follows:

[0081] lambda edir edn eup uavgdir uavgdn uavgup

[0082] phi(0)...phi(m)

[0083] umu(0)u0u(umu(0))uu(umu(0),phi(0))...uu(umu(0),phi(m))

[0084] umu(1)u0u(umu(1))uu(umu(1),phi(0))...uu(umu(1),phi(m))

[0085] umu(2)u0u(umu(2))uu(umu(2),phi(0))...uu(umu(2),phi(m)) .... ....

[0088] umu(n)u0u(umu(n))uu(umu(n),phi(0))...uu(umu(n),phi(m))

[0089] The meanings of the symbols in the output file are shown in Table 3 below:

[0090] Table 3 Symbols and their meanings

[0091]

[0092]

[0093] The uu in the output file is the radiance value of the target data to be extracted in the simulation process. Then the output value is further calculated and the corresponding sun-earth distance correction is performed. The uu data in the output file is the radiance of the earth scene at the sensor entrance pupil, in units of mW / (m2 nm sr), denoted as L. Since the day_of_year command is used in the input file to perform the sun-earth distance correction, the calculation formula (2) for the high spectral reflectance at the top of the atmosphere is as follows:

[0094]

[0095] Where ρ is the high spectral reflectance at the top of the atmosphere; F S is the instantaneous solar irradiance outside the earth, the unit is mW / (m2nm); θ is the solar zenith angle; λ is the wavelength value, the unit is nm.

[0096] Instantaneous solar irradiance calculation formula (3):

[0097]

[0098] where F is the extraterrestrial solar irradiance at a distance of one astronomical unit (AU) and d is the distance between the Sun and the Earth in AU.

[0099] Finally, the high spectral reflectance of the top of atmosphere is obtained by the following calculation formula (4):

[0100]

[0101] Experiments usually simulate in batches according to bands to obtain ρ corresponding to multiple wavelength values, which constitutes a hyperspectral reflectance data of the top of the atmosphere, which is a vector. After the hyperspectral reflectance data of the top of the atmosphere is calculated, it must be stored in a certain way for the convenience of subsequent experiments. The hyperspectral reflectance data of the top of the atmosphere obtained by llibRadtran simulation has a high resolution. The hyperspectral reflectance value of the top of the atmosphere corresponds to the wavelength. The data is stored in a .csv format file at intervals of 1nm. The target satellite sensor has multiple channels, and different spectral response functions are located in different bands. The storage method is separated by the type of land object and the band. One band corresponds to one land object and is stored in a .csv file. One line is a hyperspectral reflectance. The size of the entire file is the number of data items × the band length. The extracted data is used as inversion model input and is named in a certain format when stored. For example, the file name is S_20Seawater30Sand10Ice_401_423.csv, which means that this file consists of 20 marine feature types, 30 desert feature types and 10 glacier feature type data, and its band range is 401nm-423nm.

[0102] In summary, the present application can also determine the simulation scheme through multi-ground object simulation data, which involves data provided by llibRadtran and data that needs to be prepared by the user. The simulation schemes for three types of ground objects are shown in Table 4 below:

[0103] Table 4 Simulation scheme of top-of-atmosphere hyperspectral reflectance for ocean, desert and glacier features

[0104]

[0105]

[0106] The input data marked with "*" in Table 4 are provided by llibRadtran. The rest need to be prepared by yourself. The data sources are introduced in Table 5 below:

[0107] Table 5 Self-prepared data source information

[0108]

[0109]

[0110] Among them, the self-prepared data needs to be preprocessed before use. First, the data is converted according to the scale factor and added offset provided in the product data. Then the data of the corresponding field is extracted according to the required spatiotemporal conditions. Since the research is aimed at scenes with clear and cloudless skies during the day, satellite observation product data requires more stringent quality control screening to ensure the quality of the observation data. The filtering strategies include: day and night filtering, water depth and altitude filtering, coordinate filtering, angle filtering, center stability filtering and cloud filtering. For desert and glacial features, the filtering will be different, such as water depth and altitude filtering, filtering out non-land areas; coordinate filtering, through coordinate filtering to lock in target research areas such as mature observation stations, and the center stability will be filtered with stricter stability standards.

[0111] The linear programming method is a type of constrained optimization method. It is a mathematical method and theory for finding the extreme value of a linear objective function under given linear constraints. The linear programming model includes an objective function and a set of constraints. The solution goal of the inversion method based on linear programming is to minimize the absolute value of the satellite reflectivity residual. In addition, according to prior information, the degree of degradation of the spectral response function of each channel of the satellite studied in this application will not be too large, that is, the difference between the on-orbit spectral response function obtained by model inversion and the spectral response function before launch will not be too large. In order to avoid the problem of ill-conditioned inversion and improve the robustness of the model, the L1 normal form regularization term is added, which is actually a background constraint. The specific representation is shown in formula (5):

[0112]

[0113] Among them, rsrλ is the spectral response function obtained by inversion, that is, the optimal solution when the objective function is minimized; rsr λ 0 represents the spectral response function measured in the laboratory before launch; the setting strategy of v is to divide the objective function into two parts, allocate the energy of the overall objective function to the original residual term and the regularization term, so that v represents the proportion of the regularization term in the overall energy, which is collectively referred to as the energy proportional factor below.

[0114] In practical applications, multiple sets of top-of-atmosphere hyperspectral reflectance and satellite-observed reflectance are used for inversion. The specific expression is shown in formula (6):

[0115]

[0116] Among them, δ i is the absolute value of the i-th reflectance residual; the value range of i is 1 to m, and m refers to the number of groups of input atmospheric top hyperspectral reflectance data; refm i is the reflectance of the ith band. In the simulation experiment, the atmospheric top high spectral reflectance and spectral response simulated by llibRadtran are convolved with this function and refsh is used. i The physical form is shown in formula (7):

[0117]

[0118] Among them, S λ(i) It refers to the top-of-atmosphere hyperspectral reflectance of the i-th group in the m-group simulated top-of-atmosphere hyperspectral reflectance data file S; is the spectral response function after degradation described above. In the simulation experiment, is regarded as the true value of the inversion solution, that is, the true value of the on-orbit spectral response function; here λ is the band, λ a is the leftmost wavelength value of the band, λ b is the rightmost wavelength value of the band. Since only the RSR value in the discrete spectrum channel can be obtained in practice, the continuous form of Refsh is converted into a discrete form as shown in the following formula (8):

[0119]

[0120] At this time S λ is an m×n matrix, is an n×1 vector, and refsh is an m×1 vector. It can be expressed by formulas (9) to (11):

[0121]

[0122] Based on sufficient prior information of the target satellite, the following multiple constraints are set to optimize and narrow the search space so that the inversion method can quickly find the optimal solution without iterative search. The constraints of the inversion model include non-negative constraints, normalization constraints, variation constraints, and difference constraints.

[0123] Non-negative constraint: The spectral response function is mathematically expressed as a ratio. Its value space is [0,1]. The mathematical formula is (12):

[0124] 0≤rsr λ(j) ≤1; 1≤j≤n (12)

[0125] Normalization constraint: The spectral response function is normalized before being input into the inversion model, that is, the integral of the spectral response function after processing is 1. Therefore, the inversion result should also follow this constraint. The mathematical formula is expressed as (13):

[0126]

[0127] Variation constraint: Given the degradation trend of the spectral response function, the smaller the response value, the smaller the degree of degradation. Use the sigmoid function and the wavelength corresponding to the spectral response function value before emission. And the maximum response difference of adjacent wavelengths in the spectral response function value before emission To constrain the range of response values, the custom Sigmoid function is used As shown in equations (14) and (15):

[0128]

[0129]

[0130] The mathematical formula of the mutation constraint is expressed as (16):

[0131]

[0132] in, represents the jth sampling point of the spectral response function on band λ obtained in the laboratory before launch, rsr λ(j) Similar. λ(j) is the value after mutation (i.e., the target spectral response function value), and the mutation size Δrsr λ(j) That is See above formula (4-11). When it is 0.2, the custom Sigmoid function Distribution Fig.10 As shown, Fig.10In the figure, the horizontal axis is the spectral response function before launch, and the vertical axis represents the variation between the inversion result and the spectral response function before launch. Fig.11 shown.

[0133] Differential constraint: The difference between adjacent response values ​​of the spectral response function before emission and the differential factor α are used to limit the difference in the spectral response function values ​​between adjacent wavelengths. The differential factor α can be adjusted based on the experience of the spectral response function before emission. The mathematical formula is (17):

[0134]

[0135] The difference constraint can effectively limit the degree of variation of adjacent points. If there is only a variation constraint, the effect is as follows: Fig.12 shown.

[0136] Obviously, the orange lines in the above figure may be the relative positions of two adjacent points in the search results after degradation. Taking the figure as an example, the maximum difference reaches Δ max , as shown in formula (18).

[0137]

[0138] This means that a large gap between the response values ​​of adjacent points in the search results is allowed. The difference constraint and the variation constraint work together to more effectively ensure that the variation of the spectral response function is not too large, effectively reducing the search space and improving the search efficiency.

[0139] This means that a large gap between the response values ​​of adjacent points in the search results is allowed. The difference constraint and the variation constraint work together to more effectively ensure that the variation of the spectral response function is not too large, effectively reduce the search space, and improve the search efficiency. The effect is as follows: Fig.13 shown.

[0140] In summary, the objective function of the inversion method based on linear programming is expressed as follows (19):

[0141]

[0142] Z is the final objective function value, which is the sum of i δ. The constraint fusion is summarized as follows (20):

[0143]

[0144] Among them, α represents the difference factor, which can be adjusted according to different inversion target channels, rsr λ(j) represents the jth sampling point of the spectral response function, rsr λ(k) represents the kth sampling point of the spectral response function, is the variation size, and n means that the spectral response function has n sampling points.

[0145] Preferably, when setting the values ​​of the energy scale factor v and the differential factor α parameters of the inversion method, the energy scale factor v is usually set to be approximately as small as possible when the research scene noise is ideal; the size of the differential factor α is set according to the maximum value of the difference between adjacent wavelength response values ​​of the spectral response function of the inversion target band.

[0146] After calling the linear programming solver to train the model, the model corresponding to the optimal solution of the spectral response function is automatically searched and output as the trained model. When the trained model is obtained, the quality of the inversion result will be judged according to the evaluation indicators, and the curve comparison of the spectral response function before and after inversion will be visualized. The mean absolute error (MAE), root mean square error (RMSE) and determination coefficient (R-Square, R) are used to judge the quality of the inversion results in the study. 2 ) are the three typical measurement indicators. Among them, MAE can better reflect the actual situation of the inversion result error; RMSE can measure the deviation between the inversion result and the initial value; R 2 Also known as goodness of fit, it is the square of the correlation coefficient, and its size determines the closeness of the inversion result to the initial value. The calculation methods of the three are shown in formulas (21) to (23):

[0147]

[0148] where y i represents the target spectral response function, that is, the spectral response function after on-orbit degradation processing, represents the spectral response function of the inversion result, i represents the i-th value of the spectral response function, and when RMSE and MAE are smaller, the inversion effect is better, and R 2 The closer it is to 1, the better the fit of the model.

[0149] In order to prove the effectiveness of the proposed method, this application sets up three experiments: accuracy experiment, universality experiment, and robustness experiment. All experimental codes are implemented by python3.8 programming, run in windows11 system, and the cpu is AMD Ryzen 7 5800H.

[0150] Accuracy experiment:

[0151] In order to verify the accuracy of the linear programming method proposed in this application, this experiment uses Figure 2-Figure 8 The experiment was conducted using the selected spectral response function and 50 pieces of top-of-atmosphere hyperspectral reflectance data of ocean features.

[0152] The inversion results of the inversion method based on linear programming are as follows: Fig.14 As shown, the evaluation indexes of the results obtained by inverting 20 groups of spectral response functions based on the linear programming inversion method are shown in Table 6 below:

[0153] Table 6 Evaluation index results of inversion method based on linear programming

[0154]

[0155]

[0156] according to Fig.14 It can be shown that the results obtained by the linear programming-based inversion method proposed in this application have a high degree of fitting, and the spectral response functions used in the experiment are all derived from the original spectral response functions measured in the laboratory before the satellite was launched. The spectral response function lines measured by some satellite sensor channels are not absolutely smooth. For more obvious fluctuation forms, the linear programming-based inversion method can be well inverted. And from the evaluation index values ​​in Table 6, the effect is also very significant, the RMSE and MAE indicators are close to 0, and the fitting degree R2 is above 99%.

[0157] Universal experiment:

[0158] In order to verify whether the inversion method based on linear programming can still obtain effective inversion results under different or multiple ground object target observation simulation conditions, four typical spectral response functions are selected from the spectral response functions used in the accuracy experiment as the spectral response function input data of this experiment. They are FY3B_MERSI_ch8 representing a single-peak narrow band, FY3B_VIRR_ch8 representing a single-peak wide band, AQUA_MODIS_ch8 representing a double-peak narrow band, and FY3C_MERSI_ch3 representing a double-peak wide band. See Table 7 for details:

[0159] Table 7 Inversion method universal experimental spectral response function information table

[0160]

[0161] The top of atmosphere hyperspectral reflectance data are set according to different types or multiple ground objects, and each spectral response function corresponds to 6 sets of top of atmosphere hyperspectral reflectance data. The accuracy experiment has described the inversion effect of 50 ocean ground object data, which is no longer involved in this experiment. If AQUA_MODIS_ch8 is taken as an example, its corresponding data file is shown in Table 8 below:

[0162] Table 8 Sample table of top of atmosphere hyperspectral reflectance data file

[0163]

[0164] The inversion result diagram is arranged according to the spectral response function corresponding to the six numbered atmospheric top spectral reflectance data in the above table as follows: Fig.15 As shown in Table 9, the evaluation indicators of the results of each experimental group are as follows:

[0165] Table 9 Results of universality experiment evaluation index

[0166]

[0167]

[0168] Judging from the spectral response function fitting diagram and the three evaluation indicators, whether the input data is different types of objects or multiple types of object target observation simulation, the inversion algorithm based on linear programming can achieve good inversion results, and the fitting degree of the experimental group is basically above 99%. The difference between the groups is that the input data is different, and the differential factor may need to adjust different values. In general, the data input method of different objects or multiple combinations of objects can achieve a good fitting effect, and the evaluation index results are similar, which fully proves that the inversion method based on linear programming proposed in this application has excellent universality.

[0169] Robustness experiment:

[0170] This experiment analyzes the robustness of the method. The robustness analysis is to evaluate whether the inversion method based on linear programming can tolerate the error of the input data. The main source of error is the deviation generated during the simulation of the hyperspectral reflectance at the top of the atmosphere. Therefore, on the basis of formula (7), the simulation deviation is added to correct the forward model. It is expressed using formula (24):

[0171]

[0172] Where η represents the simulation deviation; refsh i It is still regarded as the satellite observed reflectivity. The robustness experiment is carried out around the setting of η.

[0173] First, we need to analyze the deviation of the simulation results of the top of atmosphere hyperspectral reflectance. This experiment is designed using AQUA MODIS as an example. The bands of its 1, 3, 4, 8, 9, and 10 channels basically cover the visible light band range studied in this application. The satellite observation reflectance extracted from the MYD02 observation data product is used as refm. Its data format is the same as refsh. The simulation deviation is analyzed by analyzing refm and refsh. The experimental process is based on the top of atmosphere hyperspectral reflectance simulation process, such as Fig.16 shown.

[0174] like Fig.17As shown in the figure, this experiment selects the relevant distribution diagram of each channel in the three types of ground objects for display, where the horizontal axis is refsh and the vertical axis is refm; the title of the figure is a combination of "ground object" and "channel number"; the red line is the fitted straight line obtained by fitting using the regression equation; the right side is the quantity level of the sample point color, and the sample points are randomly extracted from the filtering results in the previous section. The relative deviation mean and relative deviation standard deviation of each channel of the three types of ground object targets are shown in the figure Fig.18 and Fig.19 ;

[0175] Fig.18 and Fig.19 In the figure, the horizontal axis represents each channel and its corresponding central wavelength, and the vertical axis rel_bias_mean and rel_bias_std represent the mean value of the relative deviation and the standard deviation of the relative deviation, which are calculated as follows:

[0176]

[0177] Depend on Fig.18 and Fig.19 It can be seen that the ocean is different from the desert and glacier, showing a negative average relative deviation relative to the true value; the standard deviation of the relative deviation of the desert is larger, which may be due to the large changes in the texture, optical properties and surface morphology of desert objects. Except for channels 1 and 10, the mean relative deviations of other bands are around or within 0.02.

[0178] According to the above analysis, the simulation deviation mainly comes from the systematic deviation of the radiation transfer model and the inaccuracy of the model input data. In order to simulate this deviation, three groups of Gaussian white noise with simulated deviation levels were set, corresponding to the relative deviation means of 0, 0.02 and 0.05, and the standard deviation of the relative deviation was 0.015. Unlike desert and glacial features, the simulation results of marine features show a negative average relative deviation relative to the real data. Therefore, when simulating the real observed reflectivity, the noise mean of desert and glacial features is negative, while the opposite is true for marine features.

[0179] In this experiment, the model input data uses the four representative spectral response functions in the previous section and the top-of-atmosphere hyperspectral reflectance of three ground object combinations (e.g., S_20Sand20Seawater10Ice_401_423csv). The robustness experimental results are shown in Fig. 20 Arranged according to the deviation level, the corresponding inversion result evaluation indicators of each spectral response function at different deviation levels are as follows in Table 10:

[0180] Table 10 Robustness experiment evaluation index results

[0181]

[0182]

[0183] from Fig. 20 It can be clearly seen that with the increase of added deviation, the fitting effect has a trend of deterioration, and the corresponding three evaluation indicators also verify this trend. At the deviation level of 0 / 0.015, the inversion effect can be well guaranteed, the fitting degree of the single-peak spectral response function is above 99%, and the fitting degree of the double-peak spectral response function is above 97%; at the deviation level of 0.02 / 0.015, the fitting degree is still very obvious, and the line has only a slight deviation; when it is expanded to 0.05 / 0.015, although the basic shape of the line still maintains a certain fitting effect, it can basically be guaranteed to be above 90% from the perspective of R2, but the line shows obvious deviation.

[0184] In general, the experiment can prove that the inversion method based on linear programming has good robustness, and the input data should be controlled below the deviation level of 0.02 / 0.015 when applied. At the same time, under the simulation scheme of the atmospheric top hyperspectral reflectance proposed in this application, the deviation of most channels is at the level of 0.02 / 0.015, which can obtain good inversion effect when inverting the spectral response function, proving that the simulation scheme is reliable and basically applicable to this inversion method.

[0185] The above embodiments are only used to illustrate the technical solutions of the present application, rather than to limit them. Although the present application has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for inverting the spectral response function in the visible light band, characterized in that: include: Through data simulation, the satellite's on-orbit spectral response function and the simulated data of the top-of-atmosphere high spectral reflectance are obtained, and the satellite's on-orbit reflectance data is calculated through the forward model; Determine the objective function of the inversion model based on the satellite's on-orbit reflection data; Determine the constraints of the inversion model; The inversion model is trained by a linear programming solver to obtain a trained inversion model; The top-of-atmosphere hyperspectral reflectance data and satellite reflectance values ​​are input into the trained inversion model to obtain the corresponding spectral response function.

2. The visible light band spectral response function inversion method according to claim 1, characterized in that: The satellite's on-orbit spectral response function is obtained through degradation simulation, and the simulated data of the high spectral reflectance at the top of the atmosphere are obtained through the libRadtran radiation transfer model simulation.

3. The visible light band spectral response function inversion method according to claim 1, characterized in that: The expression for calculating satellite on-orbit reflection data is: Among them, refm is the satellite observed reflectivity, S is the high spectral reflectivity of the top of the atmosphere, rsr represents the satellite on-orbit spectral response function, and λ represents the band.

4. The visible light band spectral response function inversion method according to claim 3, characterized in that: The objective function Z of the inversion model is: Among them, δ i Represents the absolute value of the residual reflectance of the i-th satellite. The value range of i is 1 to m. m represents the number of groups of input atmospheric top hyperspectral reflectance data. i represents the reflectivity of the i-th band, rsr λ represents the spectral response function obtained by inversion, rsr λ 0 represents the spectral response function measured in the laboratory before the satellite was launched, and v represents the energy scaling factor.

5. The visible light band spectral response function inversion method according to claim 1, characterized in that: The constraints of the inversion model include non-negativity constraints, normalization constraints, variation constraints and difference constraints.

6. The visible light band spectral response function inversion method according to claim 5, characterized in that: The constraint condition of the inversion model is expressed as: Among them, α represents the difference factor, rsr λ(j) represents the jth sampling point of the spectral response function, rsr λ(k) represents the kth sampling point of the spectral response function, is the variation size, and n means that the spectral response function has n sampling points.

Images