Radiometric correction method and system for airborne hyperspectral imager

Abstract

The application discloses an airborne hyperspectral imager radiation correction method and system, which is used for improving the radiation correction accuracy under complex weather. The method comprises the following steps: calculating the theoretical sun position based on the initial coordinates of the unmanned aerial vehicle, ephemeris data and atmospheric parameters; combining the IMU heading and wind speed and direction data to perform heading and differential compensation on the sun azimuth angle, and performing refraction correction on the elevation angle, and outputting high-precision sun position parameters; synchronously collecting hyperspectral images and irradiance data, using quaternion interpolation combined with IMU data to realize space-time alignment and eliminate attitude disturbance; constructing a hemispherical space irradiance model based on the aligned data, dividing the cloud layer mode (clear sky, thin cloud and thick cloud) according to the atmospheric optical thickness variance, and respectively adopting Lambert body proportional correction, discrete longitudinal scale method or deep convolution spectral network for radiation correction. The application significantly improves the correction accuracy and generalization ability, supports large-area image splicing, and is suitable for quantitative applications such as geological exploration and environmental monitoring.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of sunlight correction of spectral imagers, in particular to a radiation correction method and system for an airborne hyperspectral imager. BACKGROUND

[0002] An airborne hyperspectral imager is a remote sensing device that continuously samples the spectrum from visible light to short-wave infrared (400-2500 nm). Its push-broom working principle can achieve high-resolution ground spectral data acquisition. The airborne hyperspectral imager system includes a high-precision hyperspectral imager, a high-precision IMU / GPS system, and a high-performance data acquisition unit, and has formed typical application cases in the fields of geological exploration, agricultural monitoring, environmental protection, etc.

[0003] The airborne hyperspectral imager is easily disturbed by factors such as changes in solar elevation angle (daytime fluctuations up to 60°), atmospheric turbulence (causing ±35% fluctuations in light intensity), and cloud cover (causing a 40-80% drop in irradiance). This causes spectral line distortion. SUMMARY

[0004] In order to solve the problems of the prior art, the purpose of the present application is to provide a radiation correction method and system for an airborne hyperspectral imager, which jointly corrects the solar azimuth and elevation angles, corrects the atmospheric refraction based on atmospheric parameters, classifies the cloud state into three categories of clear sky, thin cloud and thick cloud based on the mode classification strategy of atmospheric optical thickness variance, and adaptively selects the correction model, thereby avoiding spectral line distortion of sunlight data.

[0005] To achieve the above purpose, the technical solution adopted by the present application is as follows:

[0006] The present application provides a radiation correction method for an airborne hyperspectral imager, which comprises:

[0007] S101, obtaining the initial coordinates of the unmanned aerial vehicle, pre-storing ephemeris data, and obtaining atmospheric parameters;

[0008] S102, calculating the theoretical solar azimuth and the theoretical solar elevation angle according to the ephemeris data and the initial coordinates of the unmanned aerial vehicle;

[0009] S103, obtaining the IMU heading data and the wind speed and direction data of the unmanned aerial vehicle, compensating the theoretical solar azimuth based on the IMU heading data, compensating the theoretical solar azimuth based on the wind speed and direction data, obtaining the corrected solar azimuth, correcting the theoretical solar elevation angle based on the atmospheric parameters, and obtaining the corrected solar elevation angle;

[0010] S104, collecting hyperspectral images by the airborne hyperspectral imager, and synchronously collecting irradiance data by the irradiance sensor;

[0011] S105, aligning the hyperspectral images and irradiance data collected at different times in space based on the IMU heading data through a quaternion interpolation compensation method;

[0012] S106, based on the aligned hyperspectral images and irradiance data and the corrected solar azimuth angle and the corrected solar elevation angle, performing hemispherical space irradiance modeling to obtain a hemispherical space irradiance model;

[0013] S107, calculating an atmospheric optical thickness value and a variance of the atmospheric optical thickness values of a plurality of consecutive frames based on the atmospheric parameters, dividing the cloud layer mode into a clear sky mode, a thin cloud mode and a thick cloud mode based on the atmospheric optical thickness value and the variance, if it is the clear sky mode, using a proportional correction hemispherical space irradiance model based on the Lambertian assumption; if it is the thin cloud mode, using a discrete ordinate method to correct the hemispherical space irradiance model; if it is the thick cloud mode, using a deep convolution spectral network to correct the hemispherical space irradiance model.

[0014] As a preferred technical solution, in step S103, the theoretical solar azimuth angle is compensated based on the IMU heading data and the theoretical solar azimuth angle is compensated based on the wind speed and direction data, including: calculating a solar vector IMU heading compensation term, calculating a solar vector wind speed differential compensation term, and calculating a corrected solar azimuth angle according to the theoretical solar azimuth angle, the solar vector IMU heading compensation term and the solar vector wind speed differential compensation term; in step S103, the theoretical solar elevation angle is corrected based on the atmospheric parameters to obtain a corrected solar elevation angle, including: calculating a solar vector atmospheric refraction correction term, correcting the solar vector atmospheric refraction correction term in combination with the atmospheric parameters, and calculating a corrected solar elevation angle according to the theoretical solar elevation angle and the corrected solar vector atmospheric refraction correction term.

[0015] As a preferred technical solution, in step S103, the airborne hyperspectral imager and the irradiance sensor are carried on the unmanned aerial vehicle gimbal, and the unmanned aerial vehicle gimbal can perform pitching motion and rotating motion; in step S104, the unmanned aerial vehicle gimbal simultaneously performs the fuzzy PID control method, the fuzzy PID control method comprising: initializing PID parameters; reading the current azimuth and pitch angle of the unmanned aerial vehicle gimbal as feedback values, reading the corrected sun azimuth and sun elevation angle as set points, and reading the current wind speed and wind direction data; calculating the control error and error change rate of the azimuth and pitch angle of the unmanned aerial vehicle gimbal; based on the wind speed and wind direction data, calculating the disturbance torque of the wind resistance on the unmanned aerial vehicle gimbal, and converting the disturbance torque into a compensation amount of the azimuth and pitch angle offset of the unmanned aerial vehicle gimbal; inputting the control error and error change rate into a fuzzy controller, and the fuzzy controller adjusting the PID parameters in real time according to a fuzzy rule base; and the fuzzy controller calculating a control amount based on the adjusted PID parameters and the compensation amount of the azimuth and pitch angle offset of the unmanned aerial vehicle gimbal, and converting the control amount into a driving signal for the pitching motion and rotating motion of the unmanned aerial vehicle gimbal.

[0016] As a preferred technical solution, in step S104, the irradiance sensor adopts a three-channel irradiance sensor, which collects visible light channel irradiance data, near-infrared channel irradiance data and short-wave infrared channel irradiance data; in step S107, the weights of the visible light channel irradiance data, the near-infrared channel irradiance data and the short-wave infrared channel irradiance data are adjusted according to the atmospheric optical thickness value and the cloud layer mode, a three-channel weighted fusion function is established, the weighted fused irradiance data is calculated, and the hemispherical space irradiance model is reconstructed according to the weighted fused irradiance data.

[0017] As a preferred technical solution, in step S105, the four-element number interpolation compensation method is used to align the hyperspectral images and the irradiance data collected at different times in space based on the IMU heading data, comprising: using the four-element number interpolation compensation method to solve the problem of IMU attitude jitter: IMU original data preprocessing, including accelerometer calibration and gyroscope denoising; real-time quaternion solving; key point interpolation compensation; through field of view projection transformation, realizing the spatial alignment of the hyperspectral images and the irradiance data: establishing a dynamic projection plane; performing pixel-level coordinate conversion of the hyperspectral images and the irradiance data and geographic coordinate projection; and realizing irradiance field resampling based on the bilinear interpolation method.

[0018] As a preferred technical solution, in step S106, the hemispherical space irradiance modeling includes: dividing the upper hemispherical space into several layers of zenith angle x several layers of azimuth angle; calculating the solid angle weight of each grid cell; independently processing the visible light, near-infrared and short-wave infrared three bands; subtracting the diffuse background radiation from the total irradiance based on the Lambertian assumption initial value, locking the direct light core area through the solar direction vector; using a 4th order spherical harmonic function to decompose the spatial radiation distribution; injecting the cloud optical thickness parameter to correct the scattering intensity; selecting the phase function model according to the cloud thickness; enabling the gradient smoothing transition algorithm in the cloud edge area; calculating the local incidence angle of the grid; scaling the irradiance intensity according to the slope; attenuating the radiation in the full shadow area to a certain proportion of the reference value; using a linear attenuation model in the half-shadow transition area.

[0019] As a preferred technical solution, in step S107, the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption includes: triggering the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption when the optical thickness is less than a preset threshold and the irradiance data variance is less than a preset threshold; synchronously collecting irradiance data in real time; calculating the top-of-atmosphere irradiance based on the solar constant correction model and compensating the atmospheric transmittance; generating a proportion correction coefficient to normalize the measured irradiance of the irradiance sensor to the ideal value at the top of the atmosphere, and triggering an abnormal alarm if the proportion correction coefficient exceeds the dynamic range constraint; calculating the pixel-level albedo; and realizing terrain effect compensation by correcting the slope / slope direction.

[0020] As a preferred technical solution, in step S107, the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption includes: triggering the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption when the optical thickness is less than a preset threshold and the irradiance data variance is less than a preset threshold; synchronously collecting irradiance data in real time; calculating the top-of-atmosphere irradiance based on the solar constant correction model and compensating the atmospheric transmittance; generating a proportion correction coefficient to normalize the measured irradiance of the irradiance sensor to the ideal value at the top of the atmosphere, and triggering an abnormal alarm if the proportion correction coefficient exceeds the dynamic range constraint; calculating the pixel-level albedo; and realizing terrain effect compensation by correcting the slope / slope direction.

[0021] As a preferred technical solution, in step S107, the adopting the deep convolution spectral network to correct the hemispherical space irradiance model comprises: constructing a deep convolution spectral network, the deep convolution spectral network comprising a multi-source input layer, a spectral-space double-flow encoder, a cross-modal attention fusion layer, a physical constraint decoder, and a reflectivity output layer; the multi-source input layer is used for inputting visible light RGB images, near-infrared band data, irradiance time series data, a solar position vector, and cloud layer optical thickness; the spectral-space double-flow encoder comprises a spectral flow and a spatial flow, the spectral flow outputs weighted coefficients of three channels of visible light / near-infrared / short-wave infrared and original observation values of the channels, and the spatial flow outputs spatial data; the spatial data and the weighted coefficients of the three channels of visible light / near-infrared / short-wave infrared and the original observation values of the channels are normalized through a cross-modal attention fusion function of the cross-modal attention fusion layer; a reflectivity physical boundary constraint layer is constructed based on the physical constraint decoder, and a band continuity loss function is established to calculate reflectivity; and the reflectivity output layer outputs the calculated reflectivity.

[0022] The application also provides an airborne hyperspectral imager radiation correction system based on the above-mentioned airborne hyperspectral imager radiation correction method, the system comprising: a positioning and ephemeris module for obtaining initial coordinates of a UAV and storing ephemeris data; an atmospheric parameter module for obtaining atmospheric parameters and calculating atmospheric optical thickness values and their continuous multi-frame variances; a solar position correction module connected to the positioning and ephemeris module and the atmospheric parameter module, for calculating a theoretical solar azimuth angle and a theoretical solar elevation angle according to the ephemeris data and the initial coordinates of the UAV, performing IMU heading compensation on the theoretical solar azimuth angle based on IMU heading data, performing wind speed differential compensation on the theoretical solar azimuth angle based on wind speed and direction data, performing atmospheric refraction correction on the theoretical solar elevation angle based on atmospheric parameters, and outputting the corrected solar azimuth angle and the solar elevation angle; a data acquisition module for acquiring hyperspectral images and irradiance data; a spatial alignment module connected to the data acquisition module, for aligning the hyperspectral images and the irradiance data based on the IMU heading data by using a quaternion interpolation compensation method; and a comprehensive correction module connected to the spatial alignment module, the solar position correction module, and the atmospheric parameter module, for constructing a hemispherical space irradiance model based on the aligned images, the irradiance data, and the corrected solar azimuth angle and the solar elevation angle, dividing a cloud layer mode according to the atmospheric optical thickness values and the variances, and calling corresponding correction models according to the cloud layer mode.

[0023] Compared with the prior art, the application has the following beneficial effects:

[0024] The airborne hyperspectral imager radiation correction method of the present application can jointly correct the solar azimuth angle and the solar elevation angle by fusing multiple data sources (such as ephemeris, IMU attitude, wind speed and direction, and real-time atmospheric parameters), significantly reducing the angle calculation error, and realizing the spatio-temporal accurate alignment of the hyperspectral image and the irradiance sensor through quaternion interpolation, solving the platform jitter and data misplacement problem. For cloud interference, a mode classification strategy based on the variance of the atmospheric optical thickness divides the cloud state into three categories: clear sky, thin cloud and thick cloud, and adaptively selects the correction model: the clear sky mode adopts the Lambertian body proportion correction, the thin cloud mode introduces the discrete ordinate method to calculate the multiple scattering effect, and the thick cloud mode uses the deep convolution spectral network to learn the nonlinear mapping relationship of the irradiance drop, effectively recovering the shielded band characteristics. The present application systematically optimizes the irradiance modeling accuracy, and provides high-fidelity data basis for agricultural vegetation monitoring, shallow water environment detection, mineral identification and other application scenarios. BRIEF DESCRIPTION OF DRAWINGS

[0025] Figure 1 A step flow chart of the airborne hyperspectral imager radiation correction method of the present application;

[0026] Figure 2 An example diagram of the present application in which the unmanned aerial vehicle collects hyperspectral images through the airborne hyperspectral imager and synchronously collects irradiance data through the irradiance sensor. DETAILED DESCRIPTION

[0027] In order to enable those skilled in the art to better understand the present application scheme, the technical solutions in the specific embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application.

[0028] As shown in Figure 1 , the present application provides an airborne hyperspectral imager radiation correction method, which comprises:

[0029] S101, obtaining the initial coordinates of the unmanned aerial vehicle, pre-storing ephemeris data, and obtaining atmospheric parameters;

[0030] S102, calculating the theoretical solar azimuth angle and the theoretical solar elevation angle according to the ephemeris data and the initial coordinates of the unmanned aerial vehicle;

[0031] S103, obtaining the IMU heading data and the wind speed and direction data of the unmanned aerial vehicle, performing IMU heading compensation on the theoretical solar azimuth angle based on the IMU heading data, performing wind speed differential compensation on the theoretical solar azimuth angle based on the wind speed and direction data, obtaining the corrected solar azimuth angle, and performing atmospheric refraction correction on the theoretical solar elevation angle based on the atmospheric parameters, obtaining the corrected solar elevation angle;

[0032] S104, as shown in Figure 2As shown, the hyperspectral image is collected by an airborne hyperspectral imager, and the irradiance data is synchronously collected by an irradiance sensor;

[0033] In S105, the hyperspectral images and the irradiance data collected at different times are spatially aligned based on the IMU heading data by using a quaternion interpolation compensation method.

[0034] In S106, based on the aligned hyperspectral image and the irradiance data and the corrected solar azimuth angle and the corrected solar elevation angle, a hemispherical space irradiance model is established to obtain a hemispherical space irradiance model.

[0035] In S107, the atmospheric optical thickness value and the variance of the atmospheric optical thickness value of a plurality of continuous frames are calculated based on the atmospheric parameters, the cloud layer mode is divided into a clear sky mode, a thin cloud mode and a thick cloud mode based on the atmospheric optical thickness value and the variance, if it is the clear sky mode, the hemispherical space irradiance model is corrected by using the proportion correction of the Lambert body assumption, if it is the thin cloud mode, the hemispherical space irradiance model is corrected by using the discrete ordinate method, and if it is the thick cloud mode, the hemispherical space irradiance model is corrected by using the deep convolution spectral network.

[0036] The airborne hyperspectral imager radiation correction method of the present application, in terms of illumination dynamic compensation, fuses multiple data sources such as ephemeris, IMU attitude, wind speed and direction and real-time atmospheric parameters, can jointly correct the solar azimuth angle and the solar elevation angle, significantly reduces the angle calculation error, and realizes the spatio-temporal accurate alignment of the hyperspectral image and the irradiance sensor by using the quaternion interpolation, solves the platform jitter and data misalignment problem. For cloud interference, the mode classification strategy based on the atmospheric optical thickness variance divides the cloud state into three categories of clear sky, thin cloud and thick cloud, and adaptively selects the correction model: the clear sky mode uses the Lambert body proportion correction, the thin cloud mode introduces the discrete ordinate method to calculate the multiple scattering effect, and the thick cloud mode uses the deep convolution spectral network to learn the nonlinear mapping relationship of the irradiance drop, effectively restores the shielded band characteristics. The present application systematically optimizes the irradiance modeling accuracy, and provides a high-fidelity data basis for agricultural vegetation monitoring, shallow water environment detection, mineral identification and other application scenarios.

[0037] The above-mentioned airborne hyperspectral imager radiation correction method is realized based on a UAV: a rigid substrate made of carbon fiber composite material is installed on the top of the UAV, and a gimbal capable of ±90° pitching motion and 360° rotation motion is installed on the rigid substrate; the airborne hyperspectral imager and the irradiance sensor are carried on the gimbal.

[0038] In the present application, the spectral range of the hyperspectral imager is 400-1000 nm, and the resolution is 6 nm; the irradiance sensor adopts a three-channel irradiance sensor of a silicon / indium gallium arsenide / tellurium cadmium mercury detector array, the silicon detector (400-700 nm; corresponding to visible light), the indium gallium arsenide detector (700-900 nm; corresponding to short-wave infrared light), and the tellurium cadmium mercury detector (900-1000 nm; corresponding to long-wave infrared light), and the spacing between the array devices is ≤50 mm, and power supply and communication synchronization are realized through a shared CANFD bus (5 Mbps).

[0039] In step S101, hardware self-checking is first performed:

[0040] The gimbal performs full-range scanning of the pitch motion of +90°→-90° and the rotation motion of 0°→360°.

[0041] During the full-range scanning performed by the gimbal, the dark current values of the hyperspectral imager and the irradiance sensor are checked.

[0042] The checking of the dark current values of the hyperspectral imager and the irradiance sensor includes: measuring the dark current in real time by using the light shielding pixels (OB pixels) outside the hyperspectral imager / irradiance sensor; and dynamically calculating the mean value of the dark current of each frame of OB region to adapt to temperature changes. .

[0043] In the present application, the dark current value of visible light needs to be less than 5 mV, the dark current value of short-wave infrared light needs to be less than 15 mV, and the dark current value of long-wave infrared light needs to be less than 25 mV. After the dark current value check is passed, the next step is performed.

[0044] In step S101, the initial position of the unmanned aerial vehicle is obtained by the GPS system of the airborne hyperspectral imager: longitude , latitude , and altitude . The ephemeris data is a positioning data set that records the orbital parameters of celestial bodies or artificial satellites in a list form, and the predetermined spatial position of the celestial bodies or artificial satellites at a specified time point is described through the coordinate data updated at a timing. In the present application, the validity period of the ephemeris data is 7 days, and the accuracy is ±0.01°. The atmospheric parameters include: temperature, humidity, and pressure.

[0045] In step S102, the initial sun vector is calculated according to the ephemeris data and the initial position of the unmanned aerial vehicle, including:

[0046] The sun declination angle and the hour angle are obtained through the ephemeris data.

[0047] The theoretical sun elevation angle is calculated:

[0048] ;

[0049] Theoretical solar azimuth angle:

[0050] ;

[0051] The initial sun vector is represented by the theoretical solar elevation angle and the theoretical solar azimuth angle.

[0052] In step S103, the IMU heading data of the UAV includes: three-axis angular velocity of the UAV , , ; three-axis acceleration of the UAV , , In this application, the IMU heading data is collected by using the IMU inertial sensor, and the sampling frequency is ≥200Hz. The wind speed is collected by using the wind speed sensor The wind speed differential term is calculated by using the first-order difference method , and the sampling frequency is ≥10Hz.

[0053] Step S103 further includes: IMU heading data and wind speed data preprocessing: aligning the IMU heading data and the wind speed data by using a hardware trigger signal, with a timestamp error <1ms; eliminating high-frequency noise of the IMU heading data by using Kalman filtering method; performing wind speed data filtering processing by using sliding window average filtering method, with a window size of 0.2s; performing zero offset correction of the gyroscope in the IMU inertial sensor; and performing gravity component stripping of the accelerometer in the IMU inertial sensor.

[0054] In step S103, the initial sun vector is compensated based on the IMU heading data, and the initial sun vector is compensated based on the wind speed differential based on the wind speed data, which includes:

[0055] Calculate the IMU heading compensation term of the sun vector: ;

[0056] Calculate the wind speed differential compensation term of the sun vector: ;

[0057] The wind speed differential compensation coefficient is , the air density is A, and the gimbal wind area is A.

[0058] The following method is used to obtain it: a dynamic model of wind speed disturbance is established, then the model parameters are calibrated by wind tunnel experiment, and finally adaptive adjustment is made in the real flight environment.

[0059] For example, in the wind tunnel experiment, the air density is 1.225 kg / m 3, the wind area A is 0.05m 2 , the wind speed v is 7.5m / s, the wind speed change rate is 5m / s 2 , the measured disturbance torque is 0.005N·m, according to the formula:

[0060] ;

[0061] can be obtained , then .

[0062]

[0063] ;

[0064] The directly calculated value is too small and has almost no effect in the compensation equation. This shows that the theoretical model is oversimplified. In fact, the gust effect is much more complex than the linear model.

[0065] Therefore, the core result of the wind tunnel experiment is to establish an empirical lookup table to fit a more practical empirical formula: , the experimenters found that, ignoring v in the formula, directly making the compensation angle proportional to , can achieve better compensation effect. Multiple linear regression is performed on the data of multiple wind tunnel experiments (under different wind speeds and different accelerations), and the optimal value of the coefficient is determined to be about 0.15.

[0066] Then the corrected solar azimuth angle is: , wherein is the prediction error compensation coefficient, which is obtained based on historical data regression analysis.

[0067] Specifically, continuous operation data of the unmanned aerial vehicle under typical working conditions of different seasons, regions and weather are collected, and the theoretical solar elevation angle, the RTK measured solar elevation angle and the turbulence intensity are recorded synchronously; an error sequence is generated to calculate the error; the solar elevation angle, the turbulence intensity, the time factor and the water vapor content characteristics are extracted; after solving the multicollinearity by Ridge Regression, a loss function is established to solve the error compensation coefficient.

[0068] For example, 12 months of historical operation data of the unmanned aerial vehicle under different working conditions are collected, and a total of 1500 groups of valid samples are obtained:

[0069] For example, the first set of data: theoretical solar elevation angle 32.5; measured solar elevation angle 32.1; turbulence intensity 0.12; water vapor content 8.2; time factor 0.35, which is collected from temperate zone;

[0070] The second set of data: theoretical solar elevation angle 48.7; measured solar elevation angle 37.9; turbulence intensity 0.18; water vapor content 6.5; time factor 0.60, which is collected from subtropical zone;

[0071] Calculate the theoretical-measured error of each set:

[0072] ; wherein is the measured solar elevation angle.

[0073] For example, the first set of data is -0.4; the second set of data is -0.8.

[0074] Then, construct the feature matrix:

[0075] For example, the first set of data feature matrix is [32.5, 0.12, 8.2, 0.35, 2]; the second set of data feature matrix is [48.7, 0.18, 6.5, 0.60, 3]. Wherein, the last element of the feature matrix represents the region, region 2 is temperate zone, and region 3 is subtropical zone.

[0076] Next, define the loss function as:

[0077] ;

[0078] Wherein, β is the coefficient vector to be solved, λ is the regularization strength, and the optimal λ = 0.1 is determined by cross-validation.

[0079] For example, the calculation is as follows:

[0080] The first set of data: intercept β0=0.0215, representing the inherent bias of the system; theoretical solar elevation angle β1=-0.0082, the higher the elevation angle, the smaller the error; turbulence intensity β2=0.1186, turbulence aggravates negative error; water vapor content β3=0.0029, water vapor increases positive error; time factor β4=-0.0357, the error decreases in the afternoon; region β5=-0.0123, representing the influence of the region.

[0081] Calculate the prediction error compensation coefficient:

[0082] ;

[0083] Where, taking typical working condition parameters: xtypical = [45.0, 0.15, 10.0, 0.50, 0, 1, 0, 0] (temperate noon), the empirical attenuation factor γ = 0.8;

[0084] Then .

[0085] Due to the refraction effect of light waves in the process of atmospheric transmission, the path is deviated in the process of light wave transmission, thereby affecting the target positioning accuracy. For high-altitude airborne photoelectric system long-distance reconnaissance, the influence of atmospheric refraction on target positioning is particularly serious. The light transmission deflection caused by atmospheric refraction consists of two parts: one part is the distance error caused by the increase of optical path caused by atmospheric refraction, and the other part is the elevation error caused by the deflection of light caused by atmospheric refraction. If the temperature, pressure and humidity information of the target point or the aircraft point are known, the distance error and elevation error caused by atmospheric refraction can be accurately calculated through the formula.

[0086] Specifically, in step S103, the atmospheric refraction correction of the initial sun vector based on the atmospheric parameters comprises:

[0087] The MODTRAN atmospheric database is called to obtain the atmospheric parameters. Exemplarily, the MODTRAN atmospheric database: MODTRAN4 (Moderate Resolution Atmospheric Transmission, synthetic moderate resolution atmospheric transmission) is a widely used software which uses a radiation transmission model to simulate the light radiation transmission process in the atmosphere. The atmospheric parameters include atmospheric pressure, air temperature, relative humidity, etc.

[0088] The atmospheric refraction correction term of the sun vector is calculated:

[0089] ;

[0090] ;

[0091] ;

[0092] ;

[0093] The atmospheric refraction correction term of the sun vector is corrected in combination with the atmospheric parameters:

[0094] ;

[0095] Wherein, is the atmospheric pressure, is the air temperature; is the relative humidity.

[0096] ​The corrected solar azimuth angle is calculated: ;

[0097] Based on the corrected solar azimuth angle and the corrected solar elevation angle , the real-time solar vector is obtained.

[0098] Further, in step 104, the unmanned aerial vehicle holder simultaneously performs a fuzzy PID control method, which includes: initializing PID parameters; reading the current azimuth angle and the current pitch angle of the unmanned aerial vehicle holder as feedback values, reading the corrected solar azimuth angle and the solar elevation angle as set points, and reading the current wind speed and wind direction data; calculating the control error and the error change rate of the azimuth angle and the pitch angle of the unmanned aerial vehicle holder; based on the wind speed and wind direction data, calculating the disturbance torque of the wind resistance on the unmanned aerial vehicle holder, and converting the disturbance torque into a compensation amount of the azimuth angle and the pitch angle offset of the unmanned aerial vehicle holder; inputting the control error and the error change rate into a fuzzy controller, and adjusting the PID parameters in real time according to a fuzzy rule base; and the fuzzy controller calculates a control amount based on the adjusted PID parameters and the compensation amount of the azimuth angle and the pitch angle offset of the unmanned aerial vehicle holder, and converts the control amount into a driving signal of the pitch movement and the rotation movement of the unmanned aerial vehicle holder. In this application, the control period of the fuzzy PID control method is 5 ms, the steady-state error is less than 0.05°, and the step response time is 0.25 seconds.

[0099] Still further, in step S104, the irradiance sensor adopts a three-channel irradiance sensor, which collects visible light channel irradiance data, near-infrared channel irradiance data and short-wave infrared channel irradiance data; in step S107, the weights of the visible light channel irradiance data, the near-infrared channel irradiance data and the short-wave infrared channel irradiance data are adjusted according to the atmospheric optical thickness value and the cloud layer mode, a three-channel weighted fusion function is established, the weighted and fused irradiance data is calculated, and the hemispherical space irradiance model is reconstructed according to the weighted and fused irradiance data.

[0100] For example, in the clear sky mode: the visible light weight starts to linearly decrease from the base value 0.40 with the increase of τ value, and the minimum is not less than 0.25, the near-infrared weight is fixed at 0.45, and the short-wave infrared weight is automatically allocated by the remaining weight so that the total weight is 1; in the thin cloud mode, the visible light weight has a base value of 0.35 and decays exponentially, the near-infrared weight starts to linearly increase from 0.40 with the increase of τ value, the short-wave infrared weight is fixed at 0.25, and finally the three weights are normalized; in the thick cloud mode, fixed weight allocation is adopted, the visible light weight is 0.15, the near-infrared weight is 0.25, and the short-wave infrared weight is 0.60.

[0101] In step S105, the hyperspectral images and irradiance data collected at different times are spatially aligned based on the IMU heading data by the quaternion interpolation compensation method, including:

[0102] The quaternion interpolation compensation method is used to solve the problem of IMU attitude jitter.

[0103] Step 1: IMU raw data preprocessing: accelerometer calibration and gyroscope denoising;

[0104] Step 2: Real-time quaternion solution:

[0105] wherein, is the IMU sampling interval, in this application, ; is the quaternion multiplication.

[0106] Step 3: Key point interpolation compensation:

[0107] wherein, is the interpolation factor, , is the midpoint time of the hyperspectral image or irradiance data at and time, is the quaternion angle, .

[0108] The hyperspectral images and irradiance data are spatially aligned by field of view projection transformation:

[0109] Step 1: Establish a dynamic projection plane:

[0110] ;

[0111] wherein, is the normal vector, ;

[0112] is the projection center (directly below the UAV);

[0113] Step 2: Pixel-level coordinate conversion:

[0114] Convert the hyperspectral image coordinates to the irradiance data coordinate system; convert the irradiance data coordinate system to the UAV coordinate system; and project the geographic coordinates;

[0115] Step 3: Realize irradiance field resampling based on bilinear interpolation method:

[0116] ,

[0117] wherein, is the weight, and is calculated based on the inverse distance of the projection plane coordinates.

[0118] In step S106, the hemispherical space irradiance modeling is performed to obtain a hemispherical space irradiance model, which includes:

[0119] Step 1, mesh division: divide the upper hemispherical space into 18 layers of zenith angles and 72 layers of azimuth angles; calculate the solid angle weight of each mesh unit;

[0120] Step 2, dynamic grouping: independently process three wave bands of visible light (400-700 nm), near-infrared (700-850 nm), and short-wave infrared (850-1000 nm).

[0121] Step 3, direct component separation: subtract the diffuse background radiation from the total irradiance based on the initial value of the Lambertian assumption, and lock the direct light core area through the solar direction vector;

[0122] Step 4, diffuse field modeling: decompose the spatial radiation distribution using a 4th order spherical harmonic function; inject the cloud optical thickness parameter to correct the scattering intensity;

[0123] Step 5, cloud scattering coupling: adaptively select the phase function model according to the cloud thickness; enable the gradient smoothing transition algorithm in the cloud edge area;

[0124] Step 6, slope correction: calculate the local incidence angle of the grid; scale the irradiance intensity according to the slope ratio;

[0125] Step 7, shadow processing: the radiation in the full shadow area (sun-terrain angle > 78°) is attenuated to 15% of the reference value; the linear attenuation model is used in the half shadow transition area (60°-78°).

[0126] Further, in step S107, the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption includes:

[0127] Step 1, when the optical thickness t < 0.1 (near cloudless state) and the irradiance data variance σ(t) < 0.05 (irradiance fluctuation is minimal), trigger the proportionally corrected hemispherical space irradiance model based on the Lambertian assumption;

[0128] Step 2, real-time synchronous acquisition of irradiance data;

[0129] Step 3, calculate the top-of-atmosphere irradiance based on the solar constant correction model, and perform atmospheric transmittance compensation;

[0130] Step 4, generate a proportionally corrected coefficient: ; thereby normalizing the measured irradiance of the irradiance sensor to the ideal value at the top of the atmosphere;

[0131] Proportionally corrected coefficient dynamic range constraint: ; and the over-limit triggers an abnormal alarm.

[0132] Step 5, calculate the reflectivity at the pixel level:

[0133] ;

[0134] wherein, is the original pixel gray value; is the dark current (measured in the initialization stage); is the system gain coefficient (pre-stored in the calibration file); is the quantization factor (4095 when 12-bit ADC);

[0135] Step 6, terrain effect compensation:

[0136] Slope / aspect correction (for non-flat farmland):

[0137] ;

[0138] wherein, is the solar incident angle (calculated based on the elevation data and the solar vector); is the local incident angle of the pixel (the terrain normal vector is generated in real time through the UAV LiDAR point cloud).

[0139] It should be noted that when <0.2 (slope > 78°), the correction of the pixel is frozen and marked as invalid; if 10 consecutive frames >1.15, automatically switch to the thin cloud mode (anti-miss judgment); three-channel independent correction to avoid band aliasing; real-time deduction of dark current (non-fixed value) to suppress thermal noise drift.

[0140] Further, in step S107, the correction of the hemispherical space irradiance model by the discrete ordinate method includes:

[0141] Step 1, parameterization of the irradiance field, output the solar incident angle , the spherical harmonic expansion order M, and the Legendre polynomial ; wherein, by default, M=4, which is sufficient under thin cloud conditions;

[0142] Step 2, dynamic discrete node optimization, first construct a node number decision function:

[0143] ;

[0144] For example, when t=0.5, N=10; when t=1.8, N=28;

[0145] Then the Gauss integral points are generated, specifically, Double-Gauss Quadrature is adopted to avoid the loss of extreme points.

[0146] Node weight factor: ; wherein .

[0147] Step 3, solve the radiation transfer equation:

[0148] ;

[0149] wherein, is the single scattering albedo (preset as 0.85-0.95); P is the scattering phase function (Henyey-Greenstein approximation is adopted);

[0150] Discretize the continuous radiation transfer equation into a 2N*2N linear equation group, and directly solve it by matrix exponential method:

[0151] ; wherein, is the coefficient matrix.

[0152] Step 4, synthesize the ground irradiance:

[0153] ;

[0154] wherein, is the direct component, is the diffuse component; is the terrain compensation term;

[0155] , wherein is the slope angle.

[0156] In step S107, the method comprises:

[0157] A deep convolutional spectral network is constructed, which comprises a multi-source input layer, a spectral-spatial double-flow encoder, a cross-modal attention fusion layer, a physical constraint decoder and a reflectance output layer.

[0158] The multi-source input layer is used to input visible light RGB image, near-infrared band data, irradiance time series data, solar position vector and cloud optical thickness. The visible light RGB image is preprocessed by adaptive histogram equalization; the near-infrared band data is preprocessed by cloud shadow enhancement filtering; the irradiance time series data is preprocessed by sliding window standardization; the solar position vector is converted from spherical coordinates to Cartesian coordinates; and the cloud optical thickness is exponentially scaled.

[0159] The spectral-spatial dual-flow encoder includes a spectral flow and a spatial flow. The spectral flow (1D convolution) adopts 5 layers of cavity convolution (expansion rate = 1, 2, 4, 8, 16) to capture long-range dependence in 400-1000 nm; the spatial flow (2D convolution) adopts a depthwise separable convolution (Depthwise Separable Conv) to extract cloud shadow spatial distribution features (convolution kernel 7x7).

[0160] The cross-modal attention fusion function of the cross-modal attention fusion layer is:

[0161] ;

[0162] wherein, is a normalized exponential function; is spatial data, from a spatial flow feature map; is a weighting coefficient of the visible light / near-infrared / short-wave infrared three channels, is the original observation value of each channel, , from a spectral flow; is a neural network activation function; the original irradiance feature is reserved for residual connection.

[0163] The physical constraint decoder includes a reflectivity physical boundary constraint layer and an inter-band continuity loss function.

[0164] wherein, the output of the reflectivity physical boundary constraint layer is calculated by using a torch.clamp function:

[0165] ;

[0166] wherein, is the output of the physical constraint decoder.

[0167] The inter-band continuity loss function is:

[0168] .

[0169] The reflectivity output layer outputs the calculated reflectivity.

[0170] The application also provides an airborne hyperspectral imager radiation correction system, which comprises a positioning and ephemeris module, an atmospheric parameter module, a sun position correction module, a data acquisition module, a spatial alignment module and a comprehensive correction module.

[0171] The positioning and ephemeris module is configured to obtain initial coordinates of the UAV and store ephemeris data. The atmospheric parameter module is configured to obtain atmospheric parameters and calculate an atmospheric optical thickness value and a variance of a plurality of continuous frames. The sun position correction module is connected to the positioning and ephemeris module and the atmospheric parameter module, and is configured to calculate a theoretical sun azimuth angle and a theoretical sun elevation angle based on the ephemeris data and the initial coordinates of the UAV, perform IMU heading compensation on the theoretical sun azimuth angle based on IMU heading data, perform wind speed differential compensation on the theoretical sun azimuth angle based on wind speed and direction data, and perform atmospheric refraction correction on the theoretical sun elevation angle based on the atmospheric parameters, and output a corrected sun azimuth angle and a corrected sun elevation angle. The data acquisition module is configured to acquire hyperspectral images and irradiance data. The spatial alignment module is connected to the data acquisition module, and is configured to align the hyperspectral images and the irradiance data based on the IMU heading data by using a quaternion interpolation compensation method. The comprehensive correction module is connected to the spatial alignment module, the sun position correction module, and the atmospheric parameter module, and is configured to construct a hemispherical space irradiance model based on the aligned images, the irradiance data, and the corrected sun azimuth angle and the corrected sun elevation angle, divide a cloud layer mode based on the atmospheric optical thickness value and the variance, and call a corresponding correction model according to the cloud layer mode.

[0172] It should be noted that the terms "first", "second", and similar terms used in the specification and claims of the present application do not denote any order, quantity, or importance, but are only used to distinguish different components. Similarly, the terms "one" or "a" and the like do not denote a quantity limitation, but mean that at least one exists. "Multiple" or "several" means at least two. Unless otherwise indicated, the terms "front", "back", "left", "right", "lower", and / or "upper" and the like are used for convenience only and are not limiting to a particular position or spatial orientation. The terms "include", "comprise", and the like are meant to encompass the elements listed thereafter as well as equivalents thereof, and do not exclude other elements. The terms "connected" and "coupled" and the like are not limited to physical or mechanical connections or couplings, and can include electrical connections or couplings, whether direct or indirect.

[0173] As used in the specification and the appended claims of the present application, the singular forms "a", "an" and "the" are intended to include plural forms as well, unless the context clearly indicates otherwise. It will also be understood that the term "and / or" as used herein refers to and encompasses any or all possible combinations of one or more of the associated listed items.

[0174] It should be understood that, for those of ordinary skill in the art, modifications or changes can be made according to the above description, and all such modifications and changes shall fall within the scope of protection of the appended claims of the present application.

Claims

1. A method for solar light correction in an airborne hyperspectral imager, characterized in that, The method includes: S101: Obtain the initial coordinates of the UAV, pre-store ephemeris data, and obtain atmospheric parameters; S102, calculate the theoretical solar azimuth and theoretical solar altitude angle based on ephemeris data and the initial coordinates of the UAV; S103: Acquire the IMU heading data and wind speed and direction data of the UAV; perform IMU heading compensation on the theoretical solar azimuth angle based on the IMU heading data; perform wind speed differential compensation on the theoretical solar azimuth angle based on the wind speed and direction data to obtain the corrected solar azimuth angle; and perform atmospheric refraction correction on the theoretical solar altitude angle based on atmospheric parameters to obtain the corrected solar altitude angle. S104 acquires hyperspectral images via an airborne hyperspectral imager and simultaneously acquires irradiance data via an irradiance sensor. S105 uses a quaternion interpolation compensation method to spatially align hyperspectral images and irradiance data collected at different times based on IMU heading data. S106, based on the aligned hyperspectral image and irradiance data, the corrected solar azimuth angle and the corrected solar elevation angle, hemispherical spatial irradiance modeling is performed to obtain the hemispherical spatial irradiance model; S107 calculates the atmospheric optical thickness value and the variance of the atmospheric optical thickness value over several consecutive frames based on atmospheric parameters. Based on the atmospheric optical thickness value and variance, the cloud pattern is divided into three types: clear sky, thin cloud, and thick cloud. If it is a clear sky mode, the proportional correction hemispherical spatial irradiance model based on the Lambertian body assumption is used. If it is a thin cloud mode, the discrete ordinate method is used to correct the hemispherical spatial irradiance model. If it is a thick cloud mode, the deep convolutional spectral network is used to correct the hemispherical spatial irradiance model.

2. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S103, the step of performing IMU heading compensation on the theoretical solar azimuth based on IMU heading data and wind speed differential compensation on the theoretical solar azimuth based on wind speed and wind direction data includes: calculating the solar vector IMU heading compensation term, calculating the solar vector wind speed differential compensation term, and calculating the corrected solar azimuth based on the theoretical solar azimuth, the solar vector IMU heading compensation term, and the solar vector wind speed differential compensation term; In step S103, the atmospheric refraction correction of the theoretical solar altitude angle based on atmospheric parameters to obtain the corrected solar altitude angle includes: calculating the solar vector atmospheric refraction correction term, correcting the solar vector atmospheric refraction correction term in combination with atmospheric parameters, and calculating the corrected solar altitude angle based on the theoretical solar altitude angle and the corrected solar vector atmospheric refraction correction term.

3. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S103, the airborne hyperspectral imager and irradiance sensor are mounted on the UAV gimbal, which can perform pitch and rotation movements. In step 104, the UAV gimbal simultaneously executes a fuzzy PID control method, which includes: initializing PID parameters; reading the current azimuth and pitch angles of the UAV gimbal as feedback values, reading the corrected solar azimuth and solar altitude angles as setpoints, and reading current wind speed and direction data; calculating the control error and error rate of change of the UAV gimbal's azimuth and pitch angles; calculating the disturbance torque of wind resistance on the UAV gimbal based on the wind speed and direction data, and converting the disturbance torque into compensation amounts for the offset of the UAV gimbal's azimuth and pitch angles; inputting the control error and error rate of change into the fuzzy controller, which adjusts the PID parameters in real time according to the fuzzy rule base; and calculating the control quantity based on the adjusted PID parameters and the compensation amounts for the offset of the UAV gimbal's azimuth and pitch angles, and converting the control quantity into drive signals for the pitch and rotation motion of the UAV gimbal.

4. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S104, the irradiance sensor is a three-channel irradiance sensor, which collects irradiance data from the visible light channel, the near-infrared channel, and the short-wave infrared channel. In step S107, the weights of the visible light channel irradiance data, the near-infrared channel irradiance data, and the short-wave infrared channel irradiance data are adjusted according to the atmospheric optical thickness value and cloud pattern to establish a three-channel weighted fusion function. The weighted fused irradiance data is calculated, and the hemispherical spatial irradiance model is reconstructed based on the weighted fused irradiance data.

5. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S105, the step of spatially aligning hyperspectral images and irradiance data acquired at different times based on IMU heading data using a quaternion interpolation compensation method includes: The problem of IMU attitude jitter is solved by using quaternion interpolation compensation method: IMU raw data preprocessing, including accelerometer calibration and gyroscope denoising; real-time quaternion calculation; interpolation compensation at key time points; By transforming the field of view projection, hyperspectral imagery and irradiance data are spatially aligned: a dynamic projection surface is established; pixel-level coordinate transformation of hyperspectral imagery and irradiance data is performed, as well as geographic coordinate projection; and irradiance field resampling is achieved based on bilinear interpolation.

6. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S106, the process of modeling hemispherical spatial irradiance to obtain a hemispherical spatial irradiance model includes: dividing the upper hemispherical space into several layers of zenith angle × several layers of azimuth angle; calculating the solid angle weight of each grid cell; independently processing the visible light, near-infrared, and short-wave infrared bands; subtracting diffuse background radiation from the total irradiance based on the initial value of the Lambertian body assumption, and locking the core region of direct light through the solar direction vector; decomposing the spatial radiation distribution using a fourth-order spherical harmonic function; injecting cloud optical thickness parameters to correct the scattering intensity; adaptively selecting the phase function model based on the cloud thickness; enabling a gradient smoothing transition algorithm in the cloud edge region; calculating the local incident angle of the grid; scaling the irradiance intensity according to the slope ratio; attenuating the radiation in the full shadow region to a certain proportion of the reference value; and using a linear attenuation model in the penumbra transition region.

7. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S107, the proportionally corrected hemispherical spatial irradiance model using the Lambertian assumption includes: triggering the proportionally corrected hemispherical spatial irradiance model using the Lambertian assumption when the optical thickness is less than a preset threshold and the variance of the irradiance data is less than a preset threshold; synchronously collecting irradiance data in real time; calculating the irradiance at the top of the atmosphere based on the solar constant correction model and performing atmospheric transmittance compensation; generating a proportional correction coefficient to normalize the measured irradiance of the irradiance sensor to the ideal value at the top of the atmosphere, and triggering an abnormal alarm if the dynamic range constraint of the proportional correction coefficient is exceeded; calculating the pixel-level reflectivity; and achieving terrain effect compensation through slope and aspect correction.

8. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S107, the step of correcting the hemispherical spatial irradiance model using the discrete ordinate method includes: parameterizing the irradiance field, outputting the solar incidence angle, the spherical harmonic expansion order, and the Legendre polynomial; dynamic discrete node optimization, first constructing a node number decision function, and then generating Gaussian integral points; discretizing the continuous radiative transfer equation into a 2N×2N linear equation system, and solving the radiative transfer equation using the matrix exponent method; and synthesizing the ground irradiance based on the direct component, diffuse component, and terrain compensation term.

9. The method for solar correction of an airborne hyperspectral imager according to claim 1, characterized in that, In step S107, the step of using a deep convolutional spectral network to correct the hemispherical spatial irradiance model includes: constructing a deep convolutional spectral network, which includes a multi-source input layer, a spectral-spatial dual-stream encoder, a cross-modal attention fusion layer, a physical constraint decoder, and a reflectance output layer; the multi-source input layer is used to input visible light RGB images, near-infrared band data, irradiance time-series data, solar position vector, and cloud optical thickness; the spectral-spatial dual-stream encoder includes a spectral stream and a spatial stream, the spectral stream outputs the weighted coefficients of the three channels of visible light / near-infrared / shortwave infrared and the original observation values ​​of each channel, and the spatial stream outputs spatial data; the cross-modal attention fusion function of the cross-modal attention fusion layer is used to normalize the spatial data and the weighted coefficients of the three channels of visible light / near-infrared / shortwave infrared and the original observation values ​​of each channel; a reflectance physical boundary constraint layer is constructed based on the physical constraint decoder, and an inter-band continuity loss function is established to calculate the reflectance; the reflectance output layer outputs the calculated reflectance.

10. A solar correction system for an airborne hyperspectral imager, characterized in that, The system includes: The positioning and ephemeris module is used to obtain the initial coordinates of the UAV and store ephemeris data; The atmospheric parameters module is used to acquire atmospheric parameters and calculate atmospheric optical thickness values ​​and their variance over multiple consecutive frames. The solar position correction module connects to the positioning and ephemeris module and the atmospheric parameters module. It is used to calculate the theoretical solar azimuth and altitude angles based on ephemeris data and the initial coordinates of the UAV, perform IMU heading compensation on the theoretical solar azimuth based on IMU heading data, perform wind speed differential compensation on the theoretical solar azimuth based on wind speed and direction data, and perform atmospheric refraction correction on the theoretical solar altitude angle based on atmospheric parameters. It outputs the corrected solar azimuth and altitude angles. The data acquisition module is used to acquire hyperspectral images and irradiance data; The spatial alignment module connects to the data acquisition module and is used to align hyperspectral imagery and irradiance data based on IMU heading data using a quaternion interpolation compensation method. The integrated correction module connects the spatial alignment module, the solar position correction module, and the atmospheric parameter module. It is used to construct a hemispherical spatial irradiance model based on the aligned image, irradiance data, and the corrected solar azimuth and elevation angles. It combines atmospheric optical thickness values ​​and variance to divide cloud patterns and calls the corresponding correction model according to the cloud pattern.

Images