A method for inverting the spectral response function in the visible light band

Abstract

The present invention discloses a method for inverting the spectral response function in the visible light band, including: obtaining the spectral response function of the satellite in orbit and the simulated data of the hyperspectral reflectance at the top of the atmosphere through data simulation, and calculating the reflected data of the satellite in orbit through a forward model; determining the objective function of the inversion model according to the reflected data of the satellite in orbit; determining the constraint conditions of the inversion model; training the inversion model through a linear programming solver to obtain a trained inversion model; inputting the hyperspectral reflectance data at the top of the atmosphere and the satellite reflectance value into the trained inversion model to obtain the corresponding spectral response function. The present invention uses an inversion algorithm based on linear programming for the visible light band of the satellite and the characteristics of the satellite, makes full use of the sufficient prior information of the target satellite to formulate reasonable constraints, and obtains a method for inverting the spectral response function in the visible light band with universality and robustness.

Description

Technical Field

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

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

[0003] At present, considerable progress has been made in radiometric calibration, and the accuracy of spectral calibration has become increasingly important. It can be said that the accuracy of on-orbit spectral calibration largely affects the quality of remote sensing inversion products. When using an inaccurate instrument spectral response function, the wavelength scale may differ from the true scale by several tenths 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 spectrometers cannot calibrate the spectrum during on-orbit operation due to the lack of airborne spectral calibration equipment. Usually, these instruments use the pre-launch spectral response function, that is, the spectral response function obtained from laboratory spectral calibration. During the entire service life, it is also assumed that the performance of the spectral response function will not change. In addition, the phenomenon of degradation of the spectral response function has been found in the monitoring and flight calibration of many optical remote sensing devices. And this phenomenon is very likely caused by certain factors such as instrument aging, chemical molecular contamination in the space environment, and temperature changes of satellite instruments, as a large number of studies have shown.

[0005] In summary, due to limitations such as pre-launch laboratory equipment and environment, the spectral response function obtained from laboratory spectral calibration may be inaccurate. At the same time, during on-orbit operation, the spectral response function may degrade. These problems result in the unreliability of the spectral response function and it cannot be directly applied. Summary of the Invention

[0006] Aiming at the above deficiencies in the prior art, the method for inverting the spectral response function in the visible light band provided by the present invention solves the problem that due to limitations such as pre-launch laboratory equipment and environment, the spectral response function obtained from laboratory spectral calibration may be inaccurate; at the same time, during on-orbit operation, the spectral response function may degrade. These problems result in the unreliability of the spectral response function and it cannot be directly applied.

[0007] To achieve the above invention objective, the technical solution adopted by the present invention is: A method for inverting the spectral response function in the visible light band, including:

[0008] Through data simulation, obtain the satellite on-orbit spectral response function and the simulated data of the hyperspectral reflectance at the top of the atmosphere, and calculate the satellite reflection data through the forward model;

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

[0010] Determine the constraint conditions of the inversion model;

[0011] Train the inversion model through a linear programming solver to obtain a trained inversion model;

[0012] Input the hyperspectral reflectance data at the top of the atmosphere and the satellite reflectance value 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 simulated data of the hyperspectral reflectance at the top of the atmosphere is obtained through simulation using the libRadtran radiative transfer model.

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

[0015]

[0016] where refm is the satellite observed reflectance, S is the hyperspectral reflectance at the top of the atmosphere, rsr represents the satellite on-orbit spectral response function, and λ represents the band.

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

[0018]

[0019]

[0020] where δ i represents the absolute value of the i-th satellite reflectance residual, the value range of i is 1 to m, m represents the number of groups of the input hyperspectral reflectance data at the top of the atmosphere, refm i represents the reflectance of the i-th band, rsr λ represents the spectral response function obtained through inversion, rsr λ 0 represents the spectral response function measured in the laboratory before the satellite is launched, and v represents the energy proportionality factor.

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

[0022] Furthermore: The expressions for the constraint conditions of the inversion model are:

[0023]

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

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

[0026] For the visible light band of the satellite, according to the characteristics of the satellite, for the top-of-atmosphere hyperspectral reflectance data simulated for multi-ground object target observations, using the inversion algorithm based on linear programming, making reasonable constraints by fully utilizing the sufficient prior information of the target satellite, and having good universality and robustness. Description of the Drawings

[0027] Figure 1 It is a flowchart of the method for inverting the spectral response function in the visible light band.

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

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

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

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

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

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

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

[0035] Figure 9 It is a schematic diagram of the degradation situation of the spectral response function of each channel.

[0036] Figure 10 It is a distribution diagram of the custom Sigmoid function.

[0037] Figure 11 It is a diagram showing the mutation of the spectral response function.

[0038] Figure 12 It is a diagram showing the relationship between adjacent points under mutation constraints.

[0039] Figure 13 It is a diagram showing the relationship between adjacent points under the combined action of mutation constraints and difference constraints.

[0040] Figure 14 It is a diagram of the inversion result of the accuracy experiment.

[0041] Figure 15 It is a diagram of the inversion result of the universality experiment.

[0042] Figure 16 It is the overall flow chart of the simulation deviation analysis of the top-of-atmosphere hyperspectral reflectance.

[0043] Figure 17 It is the correlation distribution diagram of each channel among the three types of ground objects selected for the robustness experiment.

[0044] Figure 18 It is the mean value diagram of the relative deviation of each channel in the robustness experiment.

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

[0046] Figure 20 It is the inversion result diagram of the robustness and universality experiment. Detailed implementation manners

[0047] The following describes the detailed implementation manners of the present invention to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the detailed implementation manners. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions made using the concept of the present invention are within the scope of protection.

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

[0049] Through data simulation, obtain the on-orbit spectral response function of the satellite and the simulated data of the top-of-atmosphere hyperspectral reflectance, and calculate the on-orbit reflection data of the satellite through the forward model;

[0050] Determine the objective function of the inversion model according to the on-orbit reflection data of the satellite;

[0051] Determine the constraint conditions of the inversion model;

[0052] Train the inversion model through a linear programming solver to obtain a trained inversion model;

[0053] Input the top-of-atmosphere hyperspectral reflectance data and the satellite reflectance value 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 remote sensing physics inversion, the forward problem is always considered known. In this embodiment, we start with the forward problem, first describe its basic principle, and then describe in detail the inversion modeling method and technology. The expression for calculating the satellite on-orbit reflection data is as shown in Equation (1):

[0055]

[0056] where refm is the satellite observed reflectance, S is the top-of-atmosphere hyperspectral reflectance, 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 takes 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 generated by integrating and simulating the hyperspectral reflectance data and the on-orbit spectral response function. 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, when conducting the inversion experiment in this embodiment, data preparation needs to be carried out first. In the data preparation stage, we will obtain the on-orbit spectral response function through degradation simulation, and at the same time use the llibRadtran radiative transfer model to simulate and obtain the top-of-atmosphere hyperspectral reflectance data. 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 during the inversion process of the spectral response function.

[0058] Preferably, as Figures 2 - 8 shown, in this embodiment, in the visible light band, the spectral response function distributions of the MERSI / VIRR sensors of FY3 A / B / C satellites and the MODIS sensor of AQUA satellite in the visible light band are respectively shown. They are the 1st, 2nd, 3rd, 8th, 9th, 10th, 11th, 12th, 13th channels of FY3 A / B / C MERSI, the 1st, 7th, 8th, 9th channels of FY3 A / B / C VIRR, and the 1st, 3rd, 4th, 8th, 9th, 10th, 11th, 12th channels of AQUA MOIDS.

[0059] It can also be seen through Figures 2 - 8 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 (unit: nm) and response value (range: 0 to 1). The wavelength is in ascending order, and the value range is selected from 390 nm to 800 nm in the visible light range.

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

[0061] Table 1 Selected Satellite and Its Sensor Information

[0062]

[0063]

[0064] The degradation of the spectral response function of each channel is as Figure 9 shown, Figure 9 where 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 channel 12 of the MODIS sensor of the AUQA satellite 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 usage process of uvspec is divided into three parts:

[0066] A1. Preparation of input files;

[0067] A2. Invocation of the uvspec model;

[0068] A3. Reading of output files.

[0069] Among them, the core link is the preparation of input files. uvspec reads the data of each parameter according to the input files. Therefore, while preparing the input files, it is also necessary to prepare and process the input data corresponding to each parameter. In addition, after processing the output files of uvspec, the obtained result is the radiance value, and the simulated top-of-atmosphere hyperspectral reflectance result data can be obtained through subsequent calculations.

[0070] In step A1, the input of the model includes the following multiple parts. For example, 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, provided by the user of the transfer model himself; the boundary conditions can set the reflecting surfaces at the top and bottom of the atmosphere; the relevant settings for solving the radiative transfer equation solver, etc. Uvspec provides a variety of radiative transfer equation solvers; the user specifies the heliocentric distance correction parameter, etc. to provide parameters for post-processing tasks. The uvspec input file contains multiple lines of instruction information, and each line represents an input of uvspec. Each line contains an instruction name followed by the corresponding parameter value. Information can be commented through "#". The following example is an input file set with the ocean as the ground object target:

[0071]

[0072]

[0073] To ensure that the uvspec model can run correctly and calculate accurately, it is necessary to preprocess the model input data before use to meet its requirements and specifications for the input data format, and to ensure that simulation failures and errors will not occur due to non-compliant data.

[0074] The processing of the input data mainly involves extracting, converting, and saving the 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 specifications for use 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, and only the data provided by the user needs to be processed. Among them, the solar zenith angle, total ozone column (general unit: DU), altitude, heliocentric distance (days), and aerosol optical depth at 550 nm can be extracted directly from the relevant product data according to time and longitude and latitude positions, and no additional processing is required. The remaining data that requires additional processing is shown in Table 2 below:

[0075] Table 2 Description of data that requires additional processing

[0076]

[0077]

[0078] In step A2, the uvspec radiative transfer model is called using the uvspec command on the terminal. The program can be run by providing the address of the pre-prepared input file and the preset output file address to the command. After the run is completed, the target data can be found in the output file.

[0079] In step A3, extracting the radiance value from the output file, calculating the hyperspectral reflectance, and finally storing the data can be summarized as output file processing.

[0080] There are two factors affecting the output file format: the selection of the uvspec radiative transfer equation solver and whether the output azimuth angle (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)) ....

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

[0088] The meanings represented by the symbols in the output file are shown in Table 3 below:

[0089] Table 3 Symbols and Their Meaning Descriptions

[0090]

[0091]

[0092] Among them, uu in the output file is the target data radiance value to be extracted in the simulation process. Then, the output value is further calculated and the corresponding Earth-Sun distance correction is performed. The uu data in the output file is the Earth scene radiance at the sensor entrance pupil, with the unit of mW / (m2 nm sr), denoted as L. Since the Earth-Sun distance correction is performed using the day_of_year instruction in the input file, the formula for the top-of-atmosphere hyperspectral reflectance (2) is as follows:

[0093]

[0094] where ρ is the top-of-atmosphere hyperspectral reflectance; F S is the extraterrestrial instantaneous solar irradiance, with the unit of mW / (m2nm); θ is the solar zenith angle; λ is the wavelength value, with the unit of nm.

[0095] The formula for the instantaneous solar irradiance (3):

[0096]

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

[0098] Finally, the top-of-atmosphere hyperspectral reflectance is obtained through the following formula (4):

[0099]

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

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

[0102] Table 4 Simulation schemes for the top-of-atmosphere hyperspectral reflectance of ocean, desert, and glacier ground objects

[0103]

[0104]

[0105] Among the input data with the "*" symbol in Table 4 are all provided by llibRadtran, and the rest need to be prepared by oneself. For the introduction of the data sources, see Table 5 below:

[0106] Table 5 Information on the data sources to be prepared by oneself

[0107]

[0108]

[0109] Among them, the data to be prepared by oneself needs to be preprocessed before use. First, data conversion is performed according to the scale factor and added offset provided in the product data. Then, the data corresponding to the required spatio-temporal conditions is extracted. Since the research is aimed at the scenario of clear sky during the day without clouds, the satellite observation product data requires more stringent quality control screening to ensure the quality of the observation data. The filtering strategies include: day-night filtering, water depth and altitude filtering, coordinate filtering, angle filtering, center stability filtering, and cloud filtering. For desert and glacier ground objects, there are differences in filtering. For example, in water depth and altitude filtering, non-land areas are filtered out; in coordinate filtering, the target research area such as a mature observation station is locked through coordinate filtering, and more stringent filtering of the stability standard is performed for center stability, etc.

[0110] The linear programming method is one of the constrained optimization methods, which is a mathematical method and theory for finding the extreme value of a linear objective function under given linear constraint conditions. The linear programming model includes an objective function and a set of constraint conditions. The solution objective of the inversion method based on linear programming is to minimize the absolute value of the satellite reflectance residual. In addition, according to the prior information, the degradation degree of the spectral response function of each channel of the satellite studied in the present application is not too large, that is, the difference between the on-orbit spectral response function obtained by model inversion and the spectral response function before launch is not too large. To avoid the ill-posed inversion problem and improve the robustness of the model, an L1-norm regularization term is added, which is actually a background constraint. The specific representation form is shown in Equation (5):

[0111]

[0112] where 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, and hereinafter it is collectively referred to as the energy proportion factor.

[0113] In the actual application scenario, multiple groups of top-of-atmosphere hyperspectral reflectance and satellite observation reflectance will be used for inversion, and the specific representation form is shown in Equation (6):

[0114]

[0115] where, δ 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 top-of-atmosphere hyperspectral reflectance data; refm i is the reflectance of the i-th band, which is obtained by convolving the top-of-atmosphere hyperspectral reflectance simulated by llibRadtran and the spectral response function in the simulation experiment, and is represented by refsh i and the physical form is shown in Equation (7):

[0116]

[0117] where, S λ(i) refers to the top-of-atmosphere hyperspectral reflectance of the i-th group in the m groups of simulated top-of-atmosphere hyperspectral reflectance data files S; is the spectral response function after the degradation process 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, and λ b is the rightmost wavelength value of the band. Since only the rsr value in the discrete spectral channels can be actually obtained, the continuous form of refsh is converted into the discrete form as follows in Equation (8):

[0118]

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

[0120]

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

[0122] Non - negative constraint: Mathematically, the spectral response function represents a ratio. Its value range is [0, 1]. The mathematical formula is expressed as (12):

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

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

[0125]

[0126] Variation constraint: Given the degradation trend of the spectral response function, the smaller its response value, the smaller the degradation degree. Use the Sigmoid function, the pre - emission spectral response function value corresponding to the wavelength and the maximum response difference between adjacent wavelengths in the pre - emission spectral response function value to constrain the change range of the response value. The custom Sigmoid function is represented by as follows. The following equations (14) and (15):

[0127]

[0128]

[0129] The mathematical formula of the variation constraint is expressed as (16):

[0130]

[0131] Among them, represents the j - th sampling point of the spectral response function at the wavelength λ obtained in the pre - emission laboratory, rsr λ(j) is similar. rsr λ(j) is the value after variation (i.e., the target spectral response function value), and the variation size Δrsr λ(j) is as shown in the above equation (4 - 11). When is 0.2, the custom Sigmoid function is distributed as Figure 10 shown, Figure 10In it, the abscissa is the spectral response function before emission, and the ordinate represents the variation size between the inversion result and the spectral response function before emission. The variation situation is as Figure 11 shown.

[0132] Difference constraint: The difference between adjacent response values of the spectral response function before emission and the difference factor α are used to limit the variation size of the spectral response function values between adjacent wavelengths. Among them, the difference factor α can be adjusted according to the experience of the spectral response function before emission. The mathematical formula is expressed as (17):

[0133]

[0134] The difference constraint can effectively limit the variation degree of adjacent points. If there is only the variation constraint, the effect is as Figure 12 shown.

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

[0136]

[0137] This means that a huge gap in the response values of adjacent points in the search results is allowed. The combined action of the difference constraint and the variation constraint can 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.

[0138] This means that a huge gap in the response values of adjacent points in the search results is allowed. The combined action of the difference constraint and the variation constraint can 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 Figure 13 shown.

[0139] To sum up, the objective function of the inversion method based on linear programming is expressed as the following formula (19):

[0140]

[0141] Z is the final objective function value, which is the sum of i deltas. The constraint conditions are summarized as the following formula (20):

[0142]

[0143] 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, and rsr λ(k) represents the kth sampling point of the spectral response function. is the variation size, and n represents that the spectral response function has n sampling points.

[0144] Preferably, when setting the numerical values of the energy ratio factor v and the difference factor α parameter of the inversion method, the energy ratio factor v is usually set to be smaller when the noise in the research scenario is ideal; the magnitude of the difference factor α is set according to the maximum value of the difference between the response values of adjacent wavelengths of the spectral response function in the inversion target band.

[0145] After training the model by calling the linear programming solver, the model corresponding to the optimal solution of the spectral response function obtained by automatic search and output is used as the trained model. When obtaining the trained model, the quality of the inversion result will be judged according to the evaluation index, and the curve comparison diagram of the spectral response function before and after inversion will be visually displayed. In the research, the quality of the inversion result is judged by three typical measurement indexes: Mean Absolute Error (MAE), Root Mean Square Error (RMSE), and Coefficient of Determination (R-Square, R 2 ). 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 the goodness of fit, is the square of the correlation coefficient, and its magnitude determines the closeness of the relationship between the inversion result and the initial value. The calculation methods of the three are shown in formulas (21) - (23) respectively:

[0146]

[0147] where y i represents the target spectral response function, that is, the spectral response function after in-orbit degradation processing, represents the inversion result spectral response function, i represents the i-th value of the spectral response function. When RMSE and MAE are smaller, it indicates a better inversion effect, while R 2 is closer to 1, indicating a better fitting degree of the model.

[0148] To prove the effectiveness of the proposed method, this application sets three experiments: accuracy experiment, universality experiment, and robustness experiment. All experimental codes are implemented by programming in Python 3.8, running on the Windows 11 system, and the CPU is AMD Ryzen 7 5800H.

[0149] Accuracy experiment:

[0150] To verify the accuracy of the method based on linear programming proposed in this application, this experiment uses Figures 2 - 8 the selected spectral response function and the top-of-atmosphere hyperspectral reflectance data of 50 ocean features for the experiment.

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

[0152] Table 6 Evaluation Index Results of Inversion by Inversion Method Based on Linear Programming

[0153]

[0154]

[0155] According to Figure 14 it can be shown that the results obtained by the inversion method based on linear programming proposed in this application have a high degree of fitting. Moreover, the spectral response functions used in the experiment all come from the original spectral response functions measured in the laboratory before satellite launch. The spectral response function lines measured by some satellite sensor channels are not absolutely smooth. For obvious fluctuation forms, the inversion method based on linear programming can well invert them. And from the numerical values of the evaluation indexes in Table 6, the effect is also very significant. The RMSE and MAE indexes are close to 0, and the fitting degree R2 is above 99%.

[0156] Universality Experiment:

[0157] To verify whether the inversion method based on linear programming can still obtain effective inversion results under the observation simulation of different or multiple ground object targets, 4 spectral response functions with typical forms are selected from the spectral response functions used in the accuracy experiment as the input data of the spectral response functions for this experiment. They are FY3B_MERSI_ch8 representing the narrow-band of single-peak form, FY3B_VIRR_ch8 representing the wide-band of single-peak form, AQUA_MODIS_ch8 representing the narrow-band of double-peak form, and FY3C_MERSI_ch3 representing the wide-band of double-peak form. The detailed information is shown in Table 7 below:

[0158] Table 7 Information Table of Spectral Response Functions for Universality Experiment of Inversion Method

[0159]

[0160] The top-of-atmosphere hyperspectral reflectance data are set according to different or multiple ground object targets. Each spectral response function corresponds to 6 groups of top-of-atmosphere hyperspectral reflectance data. The inversion effect using 50 sets of ocean ground object data has been described in the accuracy experiment and will not be involved in this experiment. Taking AQUA_MODIS_ch8 as an example, its corresponding data file is shown in Table 8 below:

[0161] Table 8 Sample Table of Top-of-atmosphere Hyperspectral Reflectance Data File

[0162]

[0163] The inversion result diagram is arranged according to the top-of-atmosphere hyperspectral reflectance data corresponding to the 6 numbered data in Table 6 for each spectral response function as Figure 15 shown. The evaluation index situation of each experimental group is as follows in Table 9:

[0164] Table 9 Results Table of Generalization Experiment Evaluation Index

[0165]

[0166]

[0167] From the spectral response function fitting diagram and the three evaluation indexes, whether using different types of ground objects or the input data simulated by multiple ground object targets for observation, the inversion algorithm based on linear programming can achieve good inversion results, and the fitting degree of the experimental groups is basically above 99%. The difference between groups lies in the different input data, and different values of the differential factor may need to be adjusted. Generally speaking, for the data input methods of different ground objects or combinations of multiple ground objects, the achieved fitting effects are very good, and the evaluation index results are similar, fully proving that the inversion method based on linear programming proposed in this application has excellent generalization.

[0168] Robustness experiment:

[0169] 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 errors of the input data, and the main source of the errors is the deviation generated during the simulation process of the top-of-atmosphere hyperspectral reflectance. Therefore, based on formula (7), the simulation deviation is added to correct the forward model. It is represented by formula (24):

[0170]

[0171] Among them, η represents the simulation deviation; refsh i is still regarded as the satellite observation reflectance. The robustness experiment is carried out around the setting of η.

[0172] First, the analysis of the deviation of the top-of-atmosphere hyperspectral reflectance simulation results needs to be carried out. This experiment designs an experiment taking AQUA MODIS as an example. The bands of its channels 1, 3, 4, 8, 9, and 10 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, and its data form is the same as that of 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, as Figure 16 shown.

[0173] As Figure 17As shown in the figure, the relevant distribution diagrams of each channel of the three types of ground objects are selected for display in this experiment. In the figure, the abscissa is refsh and the ordinate is refm; the title of the figure is the combination of "ground object" and "channel number"; the red line is the fitting line obtained by fitting with the regression equation; on the right is the quantity level of the sample point color, and the sample points are randomly extracted from the filtering results of the previous section. The diagrams showing the mean relative deviation and the standard deviation of the relative deviation of each channel of the three types of ground object targets are as Figure 18 and Figure 19 ;

[0174] Figure 18 and Figure 19 In the figure, the abscissa is each channel and its corresponding central wavelength, and rel_bias_mean and rel_bias_std on the ordinate represent the mean relative deviation and the standard deviation of the relative deviation respectively. The calculation methods are as follows in equations (25) and (26):

[0175]

[0176] From Figure 18 and Figure 19 it can be seen that the ocean is different from the desert and the glacier, showing a negative average relative deviation with respect to the true value; the standard deviation of the relative deviation of the desert is larger, probably due to the large variations in the texture, optical properties, and surface morphology of the desert ground object. Except for channels 1 and 10, the mean relative deviation of other bands is about 0.02 or less.

[0177] According to the above analysis, the simulation deviation mainly comes from the systematic deviation of the radiative transfer model and the inaccuracy of the model input data. In order to simulate this deviation, three groups of Gaussian white noises with different simulation deviation levels are set, corresponding to the mean relative deviation of 0, 0.02, and 0.05 respectively, while the standard deviation of the relative deviation is 0.015. Different from the desert and glacier ground objects, the simulation results of the ocean ground object show a negative average relative deviation with respect to the real data. Therefore, when simulating the real observed reflectance, the noise mean of the desert and glacier ground objects is negative, while that of the ocean ground object is the opposite.

[0178] In this experiment, the model input data uses the four representative spectral response functions and the top-of-atmosphere hyperspectral reflectance of the three types of ground object combinations from the previous section (e.g., S_20Sand20Seawater10Ice_401_423csv). The results of the robustness experiment are as Figure 20 Arranged according to the deviation level, the evaluation index situations of the corresponding inversion results of each spectral response function under different deviation levels are as shown in Table 10 below:

[0179] Table 10 Results Table of Robustness Experiment Evaluation Index

[0180]

[0181]

[0182] From Figure 20 It can be clearly seen that as the added deviation increases, the fitting effect shows a deteriorating trend, 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, with the fitting degree of the single-peak spectral response function above 99% and the fitting degree of the double-peak spectral response function 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, and from the perspective of R2, it can basically be guaranteed to be above 90%, but the line shows obvious deviation.

[0183] Generally speaking, 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 as much as possible during application. At the same time, under the atmospheric top hyperspectral reflectance simulation scheme proposed in this application, the deviation of most channels is at the level of 0.02 / 0.015, and a good inversion effect can be obtained during the inversion of the spectral response function, which proves that the simulation scheme is reliable and basically applicable to this inversion method.

[0184] 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 foregoing embodiments, those of ordinary skill in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements on some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate 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; The expression for calculating satellite on-orbit reflection data is: in, represents the satellite observation reflectivity, represents the high spectral reflectance at the top of the atmosphere, represents the satellite on-orbit spectral response function, Indicates the band; The objective function of the inversion model is determined based on the satellite on-orbit reflection data; the objective function Z of the inversion model is: in, Indicates The absolute value of the satellite reflectivity residual, The value range is 1~ m , m Indicates the number of groups of input top-of-atmosphere hyperspectral reflectance data, Indicates The reflectivity of each band, represents the spectral response function obtained by inversion, represents the spectral response function measured in the laboratory before the satellite launch, v represents the energy scaling factor, Indicates that it contains m The first file of the simulated atmospheric top hyperspectral reflectance data file The top-of-atmosphere hyperspectral reflectance of the group; Determine the constraints of the inversion model: in, represents the difference factor, The spectral response function j sampling points, The spectral response function k sampling points, is the variation size, n This means that the spectral response function is n sampling points, Indicates the band obtained in the laboratory before launch The spectral response function of k sampling points, Indicates the band obtained in the laboratory before launch The spectral response function of j sampling points, Indicates the band obtained in the laboratory before launch The spectral response function of k +1 sampling point, It represents the maximum response difference between adjacent wavelengths in the spectral response function value before emission; 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 constraints of the inversion model include non-negativity constraints, normalization constraints, variation constraints and difference constraints.

Images