Quantitative inversion method of airborne mid-infrared hyperspectral remote sensing data in mining areas

Abstract

The invention belongs to the technical field of remote sensing image processing, and discloses a method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data in a mining area, including: step one: simulation analysis to select the best inversion band; step two: derive the mid-infrared radiation transmission equation; step three: call the radiation transmission model to simulate the atmospheric parameters required for inversion; step four: perform atmospheric correction on the original image; step five: estimate the initial surface temperature using the NEM module; step six: recalculate the emissivity of the remaining bands according to the initial temperature and linearly fit; step seven: construct a cost function model by fitting parameters and initial temperature and iterate to find the optimal solution for temperature; step eight: determine whether the temperature converges according to the optimal solution temperature and the initial surface temperature, and determine the inverted surface temperature and surface emissivity. The invention realizes high-precision extraction of surface temperature and emissivity in a field mining area, and can be used for field mining area detection using mid-infrared hyperspectral remote sensing technology.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing image processing, and in particular relates to a quantitative inversion method for airborne mid-infrared hyperspectral remote sensing data of a mining area. Background Art

[0002] The exploration and development of mineral resources plays a vital role in social and economic construction and development. The spectra of some ground objects such as minerals and rocks have obvious characteristics in the mid-infrared spectral region. Obtaining the emissivity spectrum of ground objects with high spatial resolution and high spectral resolution is of great significance for target identification, mineral mapping and other applications. Surface temperature is a key parameter in the study of physical changes on the earth's surface. In the infrared region, it is often coupled with emissivity to emit thermal radiation outward. How to separate surface temperature from emissivity from the radiance data collected by the remote sensor is a key issue in mid-infrared remote sensing.

[0003] Due to the coupling problem between the temperature and emissivity of the ground objects, for remote sensing observation data of N bands, it is necessary to invert 1 surface temperature value and N band emissivity values, that is, N radiation transfer equations correspond to N+1 unknowns. The solution of this underdetermined set of equations requires the introduction of other reasonable assumptions, which have been solved in the field of thermal infrared remote sensing. However, during the day, the transmission process of mid-infrared and thermal infrared radiation is essentially different. According to the blackbody radiation law, the earth's surface and the sun both emit mid-infrared radiation, and the radiation values ​​on the ground are of the same magnitude. Therefore, it is difficult to eliminate the influence of the environment and obtain accurate surface information, which leads to great difficulties in the inversion of surface parameters of mid-infrared remote sensing data.

[0004] In view of the technical defects of separating surface temperature and emissivity information from mid-infrared hyperspectral remote sensing data, it is urgent to develop an inversion method for separating mid-infrared emissivity and surface temperature, which can effectively overcome the current technical defects. Summary of the invention

[0005] The purpose of the present invention is to provide an inversion method to solve the problem of difficulty in separating mid-infrared emissivity and surface temperature in mining areas. By reasonably simplifying the complex mid-infrared radiation transmission equation, introducing a mature temperature emissivity separation algorithm and supplemented by accurate atmospheric correction, the surface temperature and mid-infrared emissivity information can be inverted to provide necessary data support for large-scale mapping and mineral identification in field mining areas, and fully explore the potential and advantages of aerial mid-infrared hyperspectral remote sensing.

[0006] The technical solution to achieve the purpose of the present invention is:

[0007] A method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data in a mining area, the method comprising the following steps:

[0008] Step 1: Perform transmittance simulation analysis on the mid-infrared region spectrum and select the inversion band combination;

[0009] Step 2: Analyze the mid-infrared radiation transmission process and establish a parameterized model of the ground-level radiation signal;

[0010] Step 3: Obtain atmospheric profile data from the National Center for Environmental Prediction of the United States, and input it into MODTRAN radiation transfer software to calculate the atmospheric parameters required for inversion;

[0011] Step 4: Perform atmospheric correction on the original image according to the atmospheric parameters to obtain the ground radiation of the target;

[0012] Step 5: Use the normalized emissivity module to set the maximum emissivity estimate and calculate the initial surface temperature for all bands;

[0013] Step 6: Recalculate the band emissivity based on the initial surface temperature, and perform piecewise linear fitting on the emissivity of all mid-infrared bands to obtain fitting parameters;

[0014] Step 7: Fitting parameters and initial temperature to construct a cost function model, and iteratively solve the local optimal solution and fitting parameters of the surface temperature corresponding to the minimum value of the cost function;

[0015] Step 8: Compare the local optimal solution of the surface temperature with the initial surface temperature, determine whether the temperature converges, and determine the inverted surface temperature and surface emissivity.

[0016] In the step 1, a continuous band combination with an atmospheric transmittance higher than 90% is selected for quantitative inversion of surface temperature and emissivity.

[0017] The transmission process of mid-infrared radiation in the atmosphere in step 2 is:

[0018]

[0019] Among them, L i is the entrance pupil radiance of sensor band i, L g_i is the ground radiation of band i at the ground observation point, τ i is the atmospheric transmittance of band i in the path from the observation point to the sensor, is the atmospheric upward radiation in band i, is the atmospheric scattered upward solar radiation in band i.

[0020] The parameterized model of the ground surface radiation signal in step 2 is:

[0021]

[0022] Among them, ε i is the emissivity of the surface, T s is the surface temperature, is the atmospheric downward radiation, is the downward solar radiation scattered by the atmosphere, is the direct irradiance of the sun on the surface, and B is the Planck function.

[0023] The atmospheric parameters required for inversion in step 3 include: atmospheric transmittance τ i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Atmospheric downward radiation Downward solar radiation scattered by the atmosphere Direct solar irradiance

[0024] The off-ground radiation of the target on the ground in step 4 is based on the atmospheric transmittance τ obtained in step 3 i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Formula (1) is used to calculate the original data L i Atmospheric correction is performed.

[0025] The calculation formula of the initial surface temperature in step 5 is:

[0026]

[0027] Where T0 is the initial surface temperature, ε max is the maximum emissivity, B -1 It is the Planck inverse function, and min means taking the minimum value of the vector.

[0028] The step six is ​​specifically as follows: based on the surface radiation equation (2), the initial surface temperature is used to recalculate the emissivity of the remaining bands, based on the piecewise linear emissivity constraint assumption, the emissivity spectrum is divided into M groups and respectively fitted by least squares linear fitting, each group has N bands (N≥3), and the linear fitting formula of the band emissivity is:

[0029] ε i =a k λ i} k ,k=1,…,M,i∈[(k-1)N+1,kN] (4)

[0030] Among them, a k , b k is the kth group of fitting parameters, λ i is the central wavelength of band i.

[0031] The calculation formula of the cost function in step seven is:

[0032]

[0033] Where E is the cost function, ‖·‖2 represents the second norm, is the off-ground radiation in band i estimated from the temperature and fitting parameters.

[0034] The above-ground radiation of band i estimated by temperature and fitting parameters in step 7 The calculation formula is:

[0035]

[0036] Among them, T is the local optimal solution of temperature when the cost function E takes the minimum value.

[0037] The step eight is specifically as follows: if the difference between the local optimal solution of the temperature and the initial surface temperature is less than a threshold, the iteration converges, and the temperature is the inverted surface temperature. The emissivity restored by the fitting parameters is the inverted emissivity. If the difference between the local optimal solution of the temperature and the initial surface temperature is greater than or equal to a threshold, and the temperature does not converge, the temperature is taken as the initial surface temperature, and steps six to eight are repeated.

[0038] The beneficial technical effects of the present invention are:

[0039] 1. The quantitative inversion method for airborne mid-infrared hyperspectral remote sensing data of mining areas provided by the present invention further explores and realizes the temperature emissivity separation of mid-infrared remote sensing data based on the research on thermal infrared inversion temperature emissivity.

[0040] 2. The quantitative inversion method of airborne mid-infrared hyperspectral remote sensing data of mining areas provided by the present invention adopts NCEP atmospheric profile data and MODTRAN radiation transmission model to accurately simulate the real atmospheric conditions at the time and location of imaging and obtain the necessary atmospheric parameters required for inversion.

[0041] 3. Mid-infrared hyperspectral remote sensing data will be interfered by solar radiation during the day, and the existing hyperspectral inversion method is not applicable. The present invention realizes the accurate inversion of surface temperature and mid-infrared emissivity during the day by accurately simulating solar radiation values.

[0042] 4. The wavelength range covered by the existing mid-infrared emissivity has been expanded. The mid-infrared emissivity information inverted by the present invention can provide a new data source for mineral identification technology in mining areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 The present invention provides a flow chart of the method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of mining areas.

[0044] FIG2 is a mid-infrared band atmospheric parameter curve calculated in Example 1 of the present invention, Figure 2a is the atmospheric transmittance curve, Figure 2bis the atmospheric upward radiation curve, Figure 2c is the atmospheric downward radiation curve, Figure 2d is the upward solar radiation curve scattered by the atmosphere, Figure 2e is the downward solar radiation curve scattered by the atmosphere, Figure 2f Direct solar radiation.

[0045] Figure 3 It is the off-ground radiation image calculated in Example 1 of the present invention.

[0046] Figure 4 are the initial values ​​of the ground surface temperature and emissivity calculated in Example 1 of the present invention.

[0047] Figure 5 It is the surface temperature and emissivity inversion value calculated in Example 1 of the present invention. DETAILED DESCRIPTION

[0048] The present invention is further described in detail below with reference to the accompanying drawings and embodiments.

[0049] like Figure 1 As shown, the quantitative inversion method of airborne mid-infrared hyperspectral remote sensing data of a mining area provided by the present invention first selects a band with a higher transmittance in the atmospheric window, performs atmospheric correction on the data, and then performs surface temperature emissivity separation processing to generate an emissivity image of the surface temperature and the mid-infrared band, which specifically includes the following steps:

[0050] Step 1: Perform transmittance simulation analysis on the mid-infrared region spectrum and select the inversion band combination;

[0051] When simulating the atmospheric transmittance in the mid-infrared region, the bands with low transmittance are discarded, and high-quality bands with atmospheric transmittance higher than 90% are selected as data for surface inversion.

[0052] Step 2: Analyze the mid-infrared radiation transmission process and establish a parameterized model of the ground-off signal;

[0053] According to the radiation transfer theory of thermodynamic equilibrium cloudless atmosphere, the radiance of the remote sensor in the mid-infrared band under cloudless daytime conditions can be expressed as:

[0054]

[0055] Among them, L i is the entrance pupil radiance of sensor band i, L g_i is the ground radiation of the ground observation point, τ i is the atmospheric transmittance of band i in the path from the observation point to the sensor, is the atmospheric upward radiation in band i, is the atmospheric scattered upward solar radiation in band i.

[0056] The surface of the mining area can be approximately regarded as a Lambertian body, and a parameterized model of the surface radiation signal is established. The radiation L g_i for:

[0057]

[0058] Among them, ε i is the emissivity of the surface, T s is the surface temperature, is the atmospheric downward radiation, is the downward solar radiation scattered by the atmosphere, is the direct irradiance of the sun on the surface, and B is the Planck function.

[0059] Step 3: Obtain the atmospheric profile data from the National Center for Environmental Prediction (NCEP) of the United States, and input it into the MODTRAN radiation transfer software to calculate the atmospheric parameters required for inversion;

[0060] The atmospheric profile data provided by the National Center for Environmental Prediction (NCEP) of the United States is downloaded using the image imaging latitude and longitude and international time. The observation geometric information such as the observation zenith angle, flight altitude, and solar zenith angle of the remote sensor are used as input data. The required atmospheric parameters are output using the MODTRAN atmospheric radiation transmission software, including: atmospheric transmittance τ i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Atmospheric downward radiation Downward solar radiation scattered by the atmosphere Direct solar irradiance

[0061] Step 4: Perform atmospheric correction on the original image based on atmospheric parameters such as atmospheric transmittance and uplink radiation value to obtain the ground radiation of the target on the ground;

[0062] The radiation information received by the sensor is affected by atmospheric attenuation and path radiation. Based on the atmospheric transmittance τ obtained in step 3 i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Formula (1) is used to calculate the original data L i Perform atmospheric correction to obtain the ground-surface radiation.

[0063] Step 5: Use the Normalized Emissivity Module (NEM) to set the maximum emissivity estimate and calculate the initial surface temperature for all bands;

[0064] Use the Normalized Emissivity Module (NEM) to set the maximum emissivity estimate ε max =0.97, used to replace the emissivity ε in formula (2)i , the surface temperature of all bands can be calculated based on formula (2), and the minimum temperature value in the surface temperature of each band is taken as the initial surface temperature T0; the calculation formula of the initial surface temperature T0 is as follows (3):

[0065]

[0066] In the formula, B -1 It is the Planck inverse function, and min means taking the minimum value of the vector.

[0067] Step 6: Recalculate the emissivity of the remaining bands according to the initial surface temperature, and perform linear fitting on the emissivity of all mid-infrared bands to obtain fitting parameters;

[0068] The emissivity of the remaining bands is recalculated according to the initial surface temperature. Based on the piecewise linear emissivity constraint assumption, all mid-infrared bands are grouped and the fitting parameters are obtained by linearly fitting each group of emissivity using the least squares method.

[0069] The initial surface temperature T0 is taken as the surface temperature T s Substituting into equation (2), we can easily obtain the emissivity spectrum ε i , the emissivity spectrum is divided into M groups and fitted linearly using the least squares method, each group has N bands (N≥3), that is, the emissivity of all bands can be obtained by piecewise fitting with fewer linear parameters. The linear fitting formula of the band emissivity is:

[0070] ε i =a k λ i +b k ,k=1,…,M,i∈[(k-1)N+1,kN] (4)

[0071] In the formula, a k , b k is the kth group of fitting parameters, λ i is the central wavelength of band i.

[0072] Step 7: Fitting parameters and initial temperature to build a cost function model, iteratively solve the temperature and fitting parameters corresponding to the minimum value of the cost function;

[0073] The cost function E is the binary norm of the error vector between the actual off-ground radiation and the estimated off-ground radiation:

[0074]

[0075] Among them, ‖·‖2 represents the two-norm, is the ground radiation of band i estimated from the temperature and fitting parameters:

[0076]

[0077] Among them, T is the local optimal solution temperature when the cost function E takes the minimum value.

[0078] The introduction of spectral piecewise linearization transforms the underdetermined system of equations in formula (2) into a system of equations with stable analytical solutions in formula (6). However, this system of equations is nonlinear and requires the use of the Newton iteration method of multivariate functions to find the root formula and gradually solve the approximate solution of the system of equations when the cost function E takes the minimum value. This approximate solution is the local optimal solution T of the surface temperature and the fitting parameter a. k , b k .

[0079] Step 8: Compare the local optimal solution of surface temperature with the initial surface temperature to determine whether the temperature converges and determine the inverted surface temperature and surface emissivity

[0080] If the difference between the local optimal solution of the surface temperature and the initial surface temperature is less than a threshold, the iteration converges and the temperature is the inverted surface temperature. The emissivity restored by the fitting parameters is the inverted emissivity. If the difference between the local optimal solution of the surface temperature and the initial surface temperature is greater than or equal to a threshold and the temperature does not converge, the temperature is taken as the initial surface temperature and steps six to eight are repeated.

[0081] Including: calculating the absolute value of the local optimal solution temperature T and the initial temperature T0, as follows:

[0082] |T-T0| <threshold (7)

[0083] Wherein, threshold is the threshold value. When the above formula is not satisfied, the iteration has not converged. Let T = T0 and repeat steps 6 to 8. When the above formula is satisfied, the iteration converges. The temperature T is the inverted surface temperature, and the surface emissivity is restored by formula (4).

[0084] Example 1: Taking a certain aerial mid-infrared image as an example, temperature and emissivity separation processing is performed. The specific processing flow is as follows: Figure 1 As shown, the following steps are included:

[0085] Step 1: According to the atmospheric transmittance at the remote sensor band response position, a band within the wavelength range of 3.4μm to 4.1μm is selected as input data for inverting surface temperature and emissivity.

[0086] Step 2: Download the atmospheric profile data from the official website of the National Center for Environmental Prediction (NCEP) of the United States, and linearly interpolate the atmospheric profile data to the time and position of image imaging, and input it into the MODTRAN atmospheric radiation transfer software in combination with the observed geometric data to calculate the required atmospheric parameters, as shown in Figure 2.

[0087] Step 3: Use atmospheric transmittance, atmospheric up-radiation and atmospheric scattered up-radiation to perform atmospheric correction on the original image to obtain the off-ground radiation data, such as Figure 3 shown.

[0088] Step 4: According to the derived mid-infrared radiation transfer equation, the maximum emissivity is set to 0.97, the normalized emissivity module (NEM) is solved, and the initial values ​​of the surface temperature and emissivity are calculated, as follows: Figure 4 shown.

[0089] Step 5: Introduce the prior knowledge of spectral piecewise linearity, iteratively solve and obtain the optimal temperature and piecewise linear parameters that minimize the cost function. This temperature is used as the inversion temperature, and the inverted emissivity is reconstructed by the piecewise linear parameters, such as Figure 5 shown.

[0090] The present invention is described in detail above with reference to the accompanying drawings and embodiments, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge of ordinary technicians in the field without departing from the purpose of the present invention. The contents not described in detail in the present invention can adopt the existing technology.

Claims

1. A method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data in mining areas, characterized in that: The method comprises the following steps: Step 1: Perform transmittance simulation analysis on the mid-infrared region spectrum and select the inversion band combination; Step 2: Analyze the mid-infrared radiation transmission process and establish a parameterized model of the ground-level radiation signal; Step 3: Obtain atmospheric profile data from the National Center for Environmental Prediction of the United States, and input it into MODTRAN radiation transfer software to calculate the atmospheric parameters required for inversion; Step 4: Perform atmospheric correction on the original image according to the atmospheric parameters to obtain the ground radiation of the target; Step 5: Use the normalized emissivity module to set the maximum emissivity estimate and calculate the initial surface temperature for all bands; Step 6: Recalculate the band emissivity based on the initial surface temperature, and perform piecewise linear fitting on the emissivity of all mid-infrared bands to obtain fitting parameters; Step 7: Fitting parameters and initial temperature to construct a cost function model, and iteratively solve the local optimal solution and fitting parameters of the surface temperature corresponding to the minimum value of the cost function; Step 8: Compare the local optimal solution of the surface temperature with the initial surface temperature, determine whether the temperature converges, and determine the inverted surface temperature and surface emissivity; The calculation formula of the cost function in step seven is: And=||(L g_i -THE' g_i )||2 (5) Among them, E is the cost function, ||·||2 represents the second norm, and L g_i is the ground radiation of band i at the ground observation point, L′ g_i is the off-ground radiation of band i estimated from the temperature and fitting parameters; The above-ground radiation L′ of band i estimated by temperature and fitting parameters in step 7 g_i The calculation formula is: Where T is the local optimal solution of the temperature when the cost function E takes the minimum value, a k 、b k is the kth group of fitting parameters, λ i is the central wavelength of band i, is the atmospheric downward radiation, is the downward solar radiation scattered by the atmosphere, is the direct irradiance of the sun on the surface, and B is the Planck function.

2. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 1, characterized in that: In the step 1, a continuous band combination with an atmospheric transmittance higher than 90% is selected for quantitative inversion of surface temperature and emissivity.

3. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 2, characterized in that: The mid-infrared radiation transmission process in step 2 is: Among them, L i is the entrance pupil radiance of sensor band i, L g_i is the ground radiation of band i at the ground observation point, τ i is the atmospheric transmittance of band i in the path from the observation point to the sensor, is the atmospheric upward radiation in band i, is the atmospheric scattered upward solar radiation in band i.

4. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 3 is characterized in that: The parameterized model of the ground surface radiation signal in step 2 is: Among them, ε i is the surface emissivity, T s is the surface temperature, is the atmospheric downward radiation, is the downward solar radiation scattered by the atmosphere, is the direct irradiance of the sun on the surface, and B is the Planck function.

5. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 4, characterized in that: The atmospheric parameters required for inversion in step 3 include: atmospheric transmittance τ i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Atmospheric downward radiation Downward solar radiation scattered by the atmosphere Direct solar irradiance 6. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 5, characterized in that: The off-ground radiation of the target on the ground in step 4 is based on the atmospheric transmittance τ obtained in step 3 i , atmospheric upward radiation Upward solar radiation scattered by the atmosphere Formula (1) is used to calculate the original data L i Atmospheric correction is performed.

7. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data in mining areas according to claim 6, characterized in that: The calculation formula of the initial surface temperature in step 5 is: Where T0 is the initial surface temperature, ε max is the maximum emissivity, B -1 It is the Planck inverse function, and min means taking the minimum value of the vector.

8. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data of a mining area according to claim 7, characterized in that: The step six is ​​specifically as follows: based on the off-ground radiation equation (2), the initial surface temperature is used to recalculate the emissivity of the remaining bands, based on the piecewise linear emissivity constraint assumption, the emissivity spectrum is divided into M groups and respectively fitted by least squares linear fitting, each group has N bands (N≥3), and the linear fitting formula of the band emissivity is: ε i =a k λ i +b k ,k=1,…,M,i∈[(k-1)N+1,kN] (4) Among them, a k 、b k is the kth group of fitting parameters, λ i is the central wavelength of band i.

9. The method for quantitative inversion of airborne mid-infrared hyperspectral remote sensing data in mining areas according to claim 8, characterized in that: The step eight is specifically as follows: if the difference between the local optimal solution of the temperature and the initial surface temperature is less than a threshold, the iteration converges, and the temperature is the inverted surface temperature. The emissivity restored by the fitting parameters is the inverted emissivity. If the difference between the local optimal solution of the temperature and the initial surface temperature is greater than or equal to a threshold, and the temperature does not converge, the temperature is used as the initial surface temperature, and steps six to eight are repeated.

Images