Airborne hyperspectral imager radiation correction method and system

Abstract

The invention discloses a radiation correction method and system for an airborne hyperspectral imager. The radiation correction method and system are used for improving the radiation correction precision in complex weather. The method comprises the following steps: calculating a theoretical sun position based on unmanned aerial vehicle initial coordinates, ephemeris data and atmospheric parameters; performing course and differential compensation on a solar azimuth angle by combining IMU course and wind speed and direction data, performing refraction correction on an elevation angle, and outputting a high-precision sun position parameter; hyperspectral images and irradiance data are synchronously acquired, quaternion interpolation is combined with IMU data to realize space-time alignment, and attitude disturbance is eliminated; a hemispherical space irradiance model is constructed based on alignment data, cloud layer modes (clear sky, thin cloud and thick cloud) are divided according to atmospheric optical thickness variance, and Lambert body proportion correction, a discrete longitudinal standard method or a deep convolutional spectral network are respectively adopted for radiation correction. The method remarkably improves the correction precision and generalization ability, supports large-area image splicing, and is suitable for quantitative application of geological exploration, environment monitoring and the like.

Description

Technical Field

[0001] This application relates to the field of solar correction technology for spectral imagers, and in particular to a radiation correction method and system for an airborne hyperspectral imager. Background Technology

[0002] An airborne hyperspectral imager is a remote sensing device that performs continuous spectral sampling in the visible to short-wave infrared band (400-2500nm). Its pushbroom operating principle enables the acquisition of high-resolution surface spectral data. An 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 been successfully applied in fields such as geological exploration, agricultural monitoring, and environmental protection.

[0003] Airborne hyperspectral imagers are susceptible to interference from factors such as changes in solar altitude angle (fluctuating by up to 60° during the day), atmospheric turbulence (causing fluctuations in light intensity of ±35%), and cloud cover (leading to a sharp drop in irradiance of 40-80%), resulting in distorted spectral data. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this application aims to provide a radiometric correction method and system for an airborne hyperspectral imager. This method corrects the solar azimuth and elevation angles, performs atmospheric refraction correction based on atmospheric parameters, classifies cloud conditions into three categories—clear sky, thin clouds, and thick clouds—based on a pattern classification strategy using atmospheric optical thickness variance, and adaptively selects a correction model to avoid distortion of solar data spectral lines.

[0005] To achieve the above objectives, this application adopts the following technical solution: This application provides a radiometric correction method for an airborne hyperspectral imager, the method comprising: 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, as well as the corrected solar azimuth angle and the corrected solar azimuth 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.

[0006] As a preferred technical solution, in step S103, the IMU heading compensation for the theoretical solar azimuth based on IMU heading data and the wind speed differential compensation for the theoretical solar azimuth based on wind speed and direction data include: 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 for 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.

[0007] As a preferred technical solution, 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 change rate 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 change rate into a fuzzy controller, which adjusts the PID parameters in real time according to a 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 movements of the UAV gimbal.

[0008] As a preferred technical solution, in step S104, the irradiance sensor adopts a three-channel irradiance sensor, which collects irradiance data of 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, a three-channel weighted fusion function is established, the weighted fused irradiance data is calculated, and the hemispherical spatial irradiance model is reconstructed based on the weighted fused irradiance data.

[0009] As a preferred technical solution, step S105, which involves aligning hyperspectral images and irradiance data acquired at different times in space using quaternion interpolation compensation based on IMU heading data, includes: addressing IMU attitude jitter using quaternion interpolation compensation; preprocessing raw IMU data, including accelerometer calibration and gyroscope denoising; performing real-time quaternion calculation; performing interpolation compensation at key time points; achieving spatial alignment of hyperspectral images and irradiance data through field-of-view projection transformation; establishing a dynamic projection surface; performing pixel-level coordinate transformation of hyperspectral images and irradiance data and geographic coordinate projection; and resampling the irradiance field based on bilinear interpolation.

[0010] As a preferred technical solution, step S106, which involves 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.

[0011] As a preferred technical solution, 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 / aspect correction.

[0012] As a preferred technical solution, step S107, which involves using the discrete ordinate method to correct the hemispherical spatial irradiance model, includes: irradiance field parameterization, outputting the solar incidence angle, spherical harmonic expansion order, and 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 exponential method; and synthesizing ground irradiance based on the direct component, diffuse component, and terrain compensation term.

[0013] As a preferred technical solution, step S107, the step of correcting the hemispherical spatial irradiance model using a deep convolutional spectral network, includes: constructing a deep convolutional spectral network, which comprises 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 comprises a spectral stream and a spatial stream, the spectral stream outputs the weighted coefficients of the three channels (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 (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.

[0014] This application also provides an airborne hyperspectral imager radiometric correction system. Based on the aforementioned airborne hyperspectral imager radiometric correction method, the system includes: a positioning and ephemeris module for acquiring the initial coordinates of the UAV and storing ephemeris data; an atmospheric parameter module for acquiring atmospheric parameters and calculating atmospheric optical thickness values ​​and their variance across multiple frames; and a solar position correction module, connected to the positioning and ephemeris module and the atmospheric parameter module, for calculating the theoretical solar azimuth and altitude angles based on ephemeris data and the UAV's initial coordinates, performing IMU heading compensation on the theoretical solar azimuth based on IMU heading data, performing wind speed differential compensation on the theoretical solar azimuth based on wind speed and direction data, and performing other functions based on the large... The system employs atmospheric parameters to correct the theoretical solar altitude angle using atmospheric refraction, outputting the corrected solar azimuth and altitude angles. A data acquisition module collects hyperspectral imagery and irradiance data. A spatial alignment module, connected to the data acquisition module, aligns hyperspectral imagery and irradiance data using quaternion interpolation compensation based on IMU heading data. A comprehensive correction module, connecting the spatial alignment module, solar position correction module, and atmospheric parameters module, constructs a hemispherical spatial irradiance model based on the aligned imagery, irradiance data, and corrected solar azimuth and altitude angles. It also combines atmospheric optical thickness and variance to define cloud patterns and calls the corresponding correction model based on the cloud pattern.

[0015] Compared with the prior art, the beneficial effects of this application are as follows: The radiometric correction method for airborne hyperspectral imagers in this application, in terms of dynamic illumination compensation, integrates multi-source data (such as ephemeris, IMU attitude, wind speed and direction, and real-time atmospheric parameters) to jointly correct the solar azimuth and elevation angles, significantly reducing angle calculation errors. Furthermore, it achieves precise spatiotemporal alignment between hyperspectral images and irradiance sensors through quaternion interpolation, resolving platform jitter and data misalignment issues. Regarding cloud interference, a pattern classification strategy based on atmospheric optical thickness variance categorizes cloud states into three types: clear sky, thin clouds, and thick clouds, and adaptively selects the correction model: the clear sky model uses Lambertian scale correction, the thin cloud model introduces the discrete ordinate method to calculate multiple scattering effects, and the thick cloud model utilizes a deep convolutional spectral network to learn the nonlinear mapping relationship of sudden irradiance drops, effectively recovering the characteristics of obscured bands. This application systematically optimizes the accuracy of irradiance modeling, providing a high-fidelity data foundation for applications such as agricultural vegetation monitoring, shallow water environment detection, and mineral identification. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the steps of a radiometric correction method for an airborne hyperspectral imager according to this application. Figure 2 This is an example diagram illustrating how the UAV in this application acquires hyperspectral images via an airborne hyperspectral imager and simultaneously acquires irradiance data via an irradiance sensor. Detailed Implementation

[0017] To enable those skilled in the art to better understand the present application, the technical solutions in specific embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings.

[0018] like Figure 1 As shown, this application provides a radiometric correction method for an airborne hyperspectral imager, the method comprising: 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, such as Figure 2 As shown, hyperspectral images are acquired using an airborne hyperspectral imager, and irradiance data is simultaneously acquired using 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, as well as the corrected solar azimuth angle and the corrected solar azimuth 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.

[0019] This application's airborne hyperspectral imager radiometric correction method, in terms of dynamic illumination compensation, integrates multi-source data, such as ephemeris, IMU attitude, wind speed and direction, and real-time atmospheric parameters, to jointly correct the solar azimuth and elevation angles, significantly reducing angle calculation errors. Furthermore, it achieves precise spatiotemporal alignment between the hyperspectral image and the irradiance sensor through quaternion interpolation, resolving platform jitter and data misalignment issues. Regarding cloud interference, a pattern classification strategy based on atmospheric optical thickness variance categorizes cloud states into three types: clear sky, thin clouds, and thick clouds, and adaptively selects the correction model: the clear sky model uses Lambertian scale correction, the thin cloud model introduces the discrete ordinate method to calculate multiple scattering effects, and the thick cloud model utilizes a deep convolutional spectral network to learn the nonlinear mapping relationship of sudden irradiance drops, effectively recovering the characteristics of obscured bands. This application systematically optimizes the accuracy of irradiance modeling, providing a high-fidelity data foundation for applications such as agricultural vegetation monitoring, shallow water environment detection, and mineral identification.

[0020] The above-mentioned airborne hyperspectral imager radiation correction method is based on a UAV: ​​a rigid substrate made of carbon fiber composite material is installed on the top of the UAV, and a gimbal that can move ±90° pitch and 360° rotation is installed on the rigid substrate; an airborne hyperspectral imager and an irradiance sensor are mounted on the gimbal.

[0021] In this application, the hyperspectral imager has a spectral range of 400-1000nm and a resolution of 6nm; the irradiance sensor adopts a three-channel irradiance sensor with a silicon / indium gallium arsenide / mercury cadmium telluride detector array, with silicon detector (400-700nm; corresponding to visible light), indium gallium arsenide detector (700-900nm; corresponding to short-wave infrared light), and mercury cadmium telluride detector (900-1000nm; corresponding to long-wave infrared light). The spacing between each device in the array is ≤50mm, and power supply and communication are synchronized through a shared CANFD bus (5Mbps).

[0022] In step S101, a hardware self-test is performed first: The gimbal performs a full-range scan of pitch motion from +90° to -90° and rotation motion from 0° to 360°. During the full-stroke scanning process of the gimbal, the dark current values ​​of the hyperspectral imager and the irradiance sensor are verified.

[0023] The calibration of the dark current values ​​of the hyperspectral imager and irradiance sensor includes: real-time measurement of dark current using the light-shielding pixels (OB pixels) around the hyperspectral imager / irradiance sensor; and dynamic calculation of the average dark current value of the OB area for each frame, adapting to temperature changes. .

[0024] In this application, the dark current value for visible light must be less than 5mV, the dark current value for short-wave infrared light must be less than 15mV, and the dark current value for long-wave infrared light must be less than 25mV. After the dark current value verification is passed, proceed to the next step.

[0025] In step S101, the initial position of the UAV is obtained through the GPS system of the airborne hyperspectral imager: longitude. ,latitude ,altitude Ephemeris data is a positioning dataset that records the orbital parameters of celestial bodies or artificial satellites in list form. It describes the predetermined spatial position of a celestial body or artificial satellite at a specified point in time through periodically updated coordinate data. In this application, the ephemeris data is valid for 7 days and has an accuracy of ±0.01°. Atmospheric parameters include temperature, humidity, and pressure.

[0026] In step S102, the initial solar vector is calculated based on the ephemeris data and the initial position of the UAV, including: Obtain the solar declination angle from ephemeris data Hour angle ; Calculate the theoretical solar altitude angle: ; Calculate the theoretical solar azimuth angle: ; The initial solar vector is represented by the theoretical solar altitude angle and the theoretical solar azimuth angle.

[0027] In step S103, the UAV's IMU heading data includes: the UAV's three-axis angular velocities. , , UAV three-axis acceleration , , In this application, an IMU inertial sensor is used to collect IMU heading data, with a sampling frequency ≥200Hz. A wind speed sensor is used to collect the actual wind speed. The wind speed differential term is calculated using the first-order difference method. The sampling frequency is ≥10Hz.

[0028] Step S103 also includes: IMU heading data and wind speed data preprocessing: aligning IMU heading data and wind speed data using hardware trigger signals with a timestamp error of <1ms; eliminating high-frequency noise in IMU heading data using Kalman filtering; filtering wind speed data using a sliding window averaging filter with a window size of 0.2s; performing zero-bias correction on the gyroscope in the IMU inertial sensor; and stripping the gravity component from the accelerometer in the IMU inertial sensor.

[0029] In step S103, IMU heading compensation is performed on the initial solar vector based on IMU heading data, and wind speed differential compensation is performed on the initial solar vector based on wind speed data, including: Calculate the heading compensation term for the solar vector IMU: ; Calculate the differential compensation term for solar vector wind speed: ; The wind speed differential compensation coefficient is... Where A is the air density and A is the windward area of ​​the gimbal.

[0030] The following method was used to obtain the model: a dynamic model of wind speed disturbance was established, the model parameters were calibrated through wind tunnel experiments, and finally adaptive adjustments were made in a real flight environment.

[0031] For example, in a wind tunnel experiment, the following is obtained: air density It is 1.225 kg / m 3 The windward area A is 0.05m². 2 The wind speed v is 7.5 m / s, and the rate of change of wind speed is... 5m / s 2 The disturbance torque was measured. It is 0.005 N·m, according to the formula: ;

[0032] achievable ,but .

[0033]

[0034] ; Directly calculated The value is too small, and it has almost no effect in the compensation equation. This indicates that the theoretical model has been oversimplified. In reality, the gust effect is far more complex than a linear model.

[0035] Therefore, the core result of the wind tunnel experiment is the establishment of an empirical lookup table to fit a more practical empirical formula: The researchers found that ignoring v in the formula and directly letting the compensation angle be the same as... A direct proportionality leads to better compensation results. Multiple linear regression was performed on multiple wind tunnel test data (under different wind speeds and accelerations) to ultimately determine the coefficients. The optimal value is approximately 0.15.

[0036] The corrected solar azimuth is as follows: ,in, The error compensation coefficient is derived from regression analysis of historical data.

[0037] Specifically, continuous operation data of UAVs under typical working conditions in different seasons, regions, and weather are collected, and theoretical solar altitude angle, RTK measured solar altitude angle, and turbulence intensity are recorded simultaneously. An error sequence is generated to calculate the error. Features of solar altitude angle, turbulence intensity, time factor, and water vapor content are extracted. After solving multicollinearity using ridge regression, a loss function is established to solve for the error compensation coefficient.

[0038] For example, 12 months of historical operational data from drones under different working conditions were collected, resulting in 1500 valid samples: For example, the first set of data: theoretical solar altitude angle 32.5°; measured solar altitude angle 32.1°; turbulence intensity 0.12°; water vapor content 8.2°; time factor 0.35°. This set of data was collected from a temperate region. Group 2 data: theoretical solar altitude angle 48.7°; measured solar altitude angle 37.9°; turbulence intensity 0.18; water vapor content 6.5%; time factor 0.60. This group of data was collected from the subtropical region. Calculate the theoretical-to-measured error for each group: ;in This is the measured solar altitude angle.

[0039] For example, the first set of data -0.4; Group 2 data It is -0.8.

[0040] Then, construct the feature matrix: For example, the first set of data feature matrices is [32.5, 0.12, 8.2, 0.35, 2]; the second set of data feature matrices is [48.7, 0.18, 6.5, 0.60, 3]. The last element of each feature matrix represents a region; region 2 is a temperate region, and region 3 is a subtropical region.

[0041] Next, the loss function is defined as: ; Where β is the coefficient vector to be determined, and λ is the regularization strength. The optimal λ = 0.1 is determined through cross-validation.

[0042] For example, the calculation yields: Group 1 data: intercept β0 = 0.0215, representing the inherent bias of the system; theoretical solar altitude angle β1 = -0.0082, the higher the altitude angle, the smaller the error; turbulence intensity β2 = 0.1186, turbulence exacerbates the negative error; water vapor content β3 = 0.0029, water vapor increases the positive error; time factor β4 = -0.0357, the error decreases in the afternoon; region β5 = -0.0123, representing regional influence.

[0043] Calculate the prediction error compensation coefficient: ; Among them, the typical operating condition parameters are: xtypical = [45.0, 0.15, 10.0, 0.50, 0, 1, 0, 0] (temperate noon), and the empirical attenuation factor γ = 0.8; but .

[0044] Because light waves undergo refraction during atmospheric transmission, their path deviates, affecting target positioning accuracy. For long-range ground reconnaissance by high-altitude airborne optoelectronic systems, the impact of atmospheric refraction on target positioning is particularly severe. The light transmission deflection caused by atmospheric refraction consists of two parts: one is the distance error due to the increased optical path length caused by atmospheric refraction, and the other is the elevation angle error caused by the deflection. If the temperature, pressure, and humidity information of the target point or the aircraft's location are known, the distance and elevation angle errors caused by atmospheric refraction can be accurately calculated using formulas.

[0045] Specifically, in step S103, the atmospheric refraction correction of the initial solar vector based on atmospheric parameters includes: Atmospheric parameters are obtained by calling the MODTRAN atmospheric database. For example, the MODTRAN atmospheric database, MODTRAN4 (Moderate Resolution Atmospheric Transmission), is a widely used software that uses a radiative transfer model to simulate the light radiation transmission process in the atmosphere. Atmospheric parameters include atmospheric pressure, air temperature, and relative humidity.

[0046] Calculate the solar vector atmospheric refraction correction term: ; ; ; ; Incorporating atmospheric parameters, the solar vector atmospheric refraction correction term is adjusted. Make corrections: ; in, Atmospheric pressure For temperature; This refers to relative humidity.

[0047] Calculate the corrected solar altitude angle: ; Based on the corrected solar azimuth angle and the corrected solar altitude angle This allows us to obtain the real-time solar vector.

[0048] Furthermore, 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 change rate 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 azimuth and pitch angle offsets of the UAV gimbal; inputting the control error and error change rate into a fuzzy controller, which adjusts the PID parameters in real time according to a fuzzy rule base; and calculating the control quantity based on the adjusted PID parameters and the compensation amounts for the azimuth and pitch angle offsets of the UAV gimbal, and converting the control quantity into drive signals for the pitch and rotation motions of the UAV gimbal. In this application, the control cycle of the fuzzy PID control method is 5ms, the steady-state error is <0.05°, and the step response time is 0.25 seconds.

[0049] Furthermore, 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, a three-channel weighted fusion function is established, the weighted fused irradiance data is calculated, and the hemispherical spatial irradiance model is reconstructed based on the weighted fused irradiance data.

[0050] For example, in the clear sky mode: the visible light weight starts from a base value of 0.40 and decreases linearly with the τ value, with a minimum of 0.25; the near-infrared weight is fixed at 0.45; and the short-wave infrared weight is automatically allocated from the remaining weights so that the total weights are 1. In the thin cloud mode, the visible light weight starts from a base value of 0.35 and decays exponentially; the near-infrared weight starts from 0.40 and increases linearly with the τ value; the short-wave infrared weight is fixed at 0.25; and finally, the three weights are summed and normalized. In the thick cloud mode, a fixed weight allocation is used, with the visible light weight at 0.15, the near-infrared weight at 0.25, and the short-wave infrared weight at 0.60.

[0051] 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: Step 1, IMU raw data preprocessing: accelerometer calibration and gyroscope noise reduction; Step 2, Real-time Quaternion Solving: ;in, For the IMU sampling interval, in this application, ; This is quaternion multiplication.

[0052] Step 3, Interpolation compensation at key time points: ;in, The interpolation factor is... , For hyperspectral images or irradiance data in and Midpoint of time, The angle between quaternions, .

[0053] By transforming the field of view projection, hyperspectral imagery and irradiance data are spatially aligned: Step 1, Create a dynamic projection surface: ; in, It is the normal vector. ; The projection center (directly below the drone); Step 2, Pixel-level coordinate transformation: Transform the hyperspectral image coordinates to the irradiance data coordinate system; transform the irradiance data coordinate system to the UAV coordinate system; perform geographic coordinate projection; Step 3: Resampling of the irradiance field based on bilinear interpolation: , in, The weight is calculated based on the inverse ratio of the distance to the projection plane coordinates.

[0054] In step S106, the process of performing hemispherical spatial irradiance modeling to obtain a hemispherical spatial irradiance model includes: Step 1, Grid generation: Divide the upper hemisphere space into 18 zenith angle layers × 72 azimuth angle layers; calculate the solid angle weight of each grid cell; Step 2, Dynamic Grouping: Independently process the three bands of visible light (400-700nm), near-infrared (700-850nm), and short-wave infrared (850-1000nm).

[0055] Step 3, Direct Component Separation: Based on the initial value of the Lambertian assumption, diffuse background radiation is subtracted from the total irradiance, and the core region of direct light is locked by the solar direction vector; Step 4, Diffuse field modeling: Decompose the spatial radiation distribution using a fourth-order spherical harmonic function; inject cloud optical thickness parameters to correct the scattering intensity; Step 5, Cloud Scattering Coupling: Adaptively select the phase function model based on cloud thickness; enable gradient smoothing transition algorithm in cloud edge regions; Step 6, Slope Correction: Calculate the local incident angle of the grid; scale the irradiance intensity according to the slope ratio; Step 7, Shadowing: Radiation in the full shadow area (solar-terrain angle > 78°) is attenuated to 15% of the baseline value; a linear attenuation model is used in the penumbra transition area (60°-78°).

[0056] Furthermore, in step S107, the proportionally corrected hemispherical spatial irradiance model using the Lambertian body assumption includes: Step 1: When the optical thickness t < 0.1 (nearly cloudless state) and the variance of the irradiance data σ(t) < 0.05 (irradiance fluctuation is very small), the proportionally corrected hemispherical spatial irradiance model using the Lambertian body assumption is triggered. Step 2: Real-time synchronous acquisition of irradiance data; Step 3: Calculate the top-of-atmosphere irradiance based on the solar constant correction model and perform atmospheric transmittance compensation; Step 4, generate the proportional correction coefficient: This normalizes the measured irradiance from the irradiance sensor to the ideal value at the top of the atmosphere. Dynamic range constraint of proportional correction factor: If the limit is exceeded, an abnormal alarm will be triggered.

[0057] Step 5, calculate pixel-level reflectivity: ; in, This represents the original pixel grayscale value. This is the dark current (measured during the initialization phase). This is the system gain coefficient (pre-saved in the calibration file); This is the quantization factor (4095 for a 12-bit ADC). Step 6, Terrain Effect Compensation: Slope / aspect correction (for non-flat farmland): ; in, The solar incidence angle (calculated based on elevation data and solar vector); The local incident angle of a pixel (generated in real time using terrain normal vectors from UAV LiDAR point clouds).

[0058] It should be noted that when When the slope is <0.2 (slope >78°), freeze the pixel correction and mark it as invalid; if it is 10 consecutive frames... >1.15, automatically switches to thin cloud mode (to prevent leakage detection); three-channel independent calibration avoids band aliasing; dark current real-time subtraction (non-fixed value) suppresses thermal noise drift.

[0059] Furthermore, in step S107, the step of correcting the hemispherical spatial irradiance model using the discrete ordinate method includes: Step 1: Parameterize the irradiance field and output the solar incidence angle. spherical harmonic expansion order M and Legendre polynomial The default value is M=4, which is sufficient under thin cloud conditions. Step 2, Dynamic Discrete Node Optimization: First, construct the node number decision function: ; For example, when t=0.5, N=10; when t=1.8, N=28; Then Gaussian integration points are generated; specifically, double-Gauss quadrature is used to avoid losing extreme points. Node weight factor: ;in .

[0060] Step 3, Solve the radiative transfer equation: ; in, is the single-scattering albedo (preset to 0.85-0.95); P is the scattering phase function (using the Henyey-Greenstein approximation); The continuous radiative transfer equations are discretized into a system of 2N×2N linear equations, which can be solved directly using the matrix exponentiation method: ;in, It is a coefficient matrix.

[0061] Step 4, synthesize ground irradiance: ; in, For direct component, For diffuse components; This is a terrain compensation item; ,in The slope angle is denoted by .

[0062] In step S107, the step of correcting the hemispherical spatial irradiance model using a deep convolutional spectral network includes: A deep convolutional spectral network is constructed, which includes a multi-source input layer, a spectral-spatial dual-stream encoder, a cross-modal attention fusion layer, a physically constrained decoder, and a reflectance output layer.

[0063] The multi-source input layer is used to input visible light RGB imagery, near-infrared band data, irradiance time-series data, solar position vector, and cloud optical thickness. Visible light RGB imagery is preprocessed using adaptive histogram equalization; near-infrared band data is preprocessed using cloud shadow enhancement filtering; irradiance time-series data is preprocessed using sliding window normalization; the solar position vector undergoes spherical coordinate to Cartesian coordinate transformation; and cloud optical thickness is exponentially scaled.

[0064] The spectral-spatial dual-stream encoder comprises a spectral stream and a spatial stream. The spectral stream (1D convolution) employs 5 layers of dilated convolutions (dilation rates = 1, 2, 4, 8, 16) to capture long-range dependencies of 400–1000 nm; the spatial stream (2D convolution) employs depthwise separable convolutions (7×7 kernels) to extract spatial distribution features of cloud shadows.

[0065] The cross-modal attention fusion function of the cross-modal attention fusion layer is: ; in, It is a normalized exponential function; The data is spatial, derived from spatial flow feature maps; These are the weighting coefficients for the visible light, near-infrared, and short-wave infrared channels. These are the raw observations from each channel. , From spectral stream; For neural network activation functions; The original irradiance characteristics are preserved for residual connections.

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

[0067] The reflectivity physical boundary constraint layer output is calculated using the torch.clamp function: ; in, This is the output of the physical constraint decoder.

[0068] The inter-band continuity loss function is: .

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

[0070] This application also provides a radiometric correction system for an airborne hyperspectral imager, which includes: a positioning and ephemeris module, an atmospheric parameter module, a solar position correction module, a data acquisition module, a spatial alignment module, and a comprehensive correction module.

[0071] The system comprises several modules: a positioning and ephemeris module for acquiring the UAV's initial coordinates and storing ephemeris data; an atmospheric parameters module for acquiring atmospheric parameters and calculating atmospheric optical thickness and its variance across multiple frames; a solar position correction module connecting the positioning and ephemeris module and the atmospheric parameters module for calculating the theoretical solar azimuth and altitude angles based on ephemeris data and the UAV's initial coordinates, performing IMU heading compensation on the theoretical solar azimuth based on IMU heading data, wind speed differential compensation on the theoretical solar azimuth based on wind speed and direction data, and atmospheric refraction correction on the theoretical solar altitude based on atmospheric parameters, outputting the corrected solar azimuth and altitude angles; a data acquisition module for acquiring hyperspectral imagery and irradiance data; and a spatial alignment module connected to the data acquisition module for aligning the hyperspectral imagery and irradiance data using quaternion interpolation compensation based on IMU heading data. 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.

[0072] It should be noted that the terms "first," "second," and similar terms used in this application specification and claims do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Similarly, "a" or "one," and similar terms do not indicate a quantity limitation, but rather indicate the presence of at least one. "A plurality" or "several" indicates at least two. Unless otherwise stated, terms such as "front," "back," "left," "right," "lower," and / or "upper" are for illustrative purposes only and are not limited to a location or spatial orientation. Terms such as "comprising" or "including" indicate that the elements or objects preceding "comprising" encompass the elements or objects listed following "comprising" or "including" and their equivalents, and do not exclude other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect.

[0073] The singular forms “a,” “the,” and “the” used in this application specification and appended claims are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used herein refers to and includes any or all possible combinations of one or more of the associated listed items.

[0074] It should be understood that those skilled in the art can make improvements or modifications based on the above description, and all such improvements and modifications should fall within the protection scope of the appended claims.

Claims

1. A radiometric correction method for 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, as well as the corrected solar azimuth angle and the corrected solar azimuth 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 radiometric correction method for 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 performing 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 radiometric correction method for 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 radiometric correction method for 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 radiometric correction method for 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 radiometric correction method for 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 radiometric correction method for 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 / aspect correction.

8. The radiometric correction method for 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 radiometric correction method for 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 radiation correction system for an airborne hyperspectral imager, characterized in that, The system is based on the airborne hyperspectral imager radiometric correction method according to any one of claims 1-9, and the system comprises: 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