Unmanned aerial vehicle-mounted multi-sensor automatic radiation calibration system and method

Abstract

The invention discloses an unmanned aerial vehicle-mounted multi-sensor automatic radiation calibration system and method, and the system comprises an unmanned aerial vehicle platform, a load integration module, a data collection module, a data processing module, and a data storage and output module, and can achieve the automatic recognition of PIF ground features, the relative normalization of multi-sensor data, and the real-time conversion of absolute reflectivity. And a high-precision and high-reliability data support is provided for quantitative analysis of the multi-source remote sensing data of the unmanned aerial vehicle in a complex field environment.

Description

Technical Field

[0001] This invention relates to the fields of UAV remote sensing technology and radiometric calibration technology, and is applicable to real-time radiometric calibration processing of UAVs equipped with various imaging sensors such as hyperspectral and multispectral sensors in dynamic field environments. Background Technology

[0002] Radiometric calibration is a crucial step in establishing the correspondence between sensor output digital quantization (DN) values ​​and the actual radiometric characteristics of ground features. It directly determines the comparability of multi-source, multi-temporal data and the reliability of quantitative analysis. Existing calibration methods mostly rely on ground calibration fields or calibration, which are not suitable for temporary field operations. Relative calibration methods based on pseudo-invariant features (PIFs) achieve radiometric normalization by identifying spectrally stable ground features. However, current technologies rely on manual selection of PIFs, resulting in low automation. Furthermore, they can only achieve relative radiometric normalization and cannot obtain absolute reflectance, thus failing to meet the requirements of high-precision analysis.

[0003] Online calibration technology using irradiance sensors on drones can only perform absolute calibration for a single sensor, failing to address the issue of data consistency across multiple sensors. Radiation bias still exists when fusing multi-source data. Therefore, a fully automated multi-sensor radiometric calibration scheme that does not require ground reference is urgently needed.

[0004] To address the aforementioned issues, this invention combines PIF relative calibration with online irradiance absolute calibration technology. Through hardware integration and algorithm innovation, it achieves automated normalization and absolute reflectance conversion of multi-sensor data, eliminating reliance on ground calibration and improving calibration performance in complex environments. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing UAV radiometric calibration methods, such as reliance on ground calibration targets, low automation, inability to adapt to dynamic lighting environments, and poor consistency of multi-sensor data. This invention provides an UAV-borne multi-sensor automatic radiometric calibration system and method. This system and method can achieve automatic identification of PIF ground features, relative normalization of multi-sensor data, and real-time conversion of absolute reflectance, providing high-precision and high-reliability data support for quantitative analysis of multi-source remote sensing data from UAVs in complex field environments.

[0006] To achieve the above objectives, this invention provides an unmanned aerial vehicle (UAV)-borne multi-sensor automatic radiometric calibration system, comprising an UAV platform, a payload integration module, a data acquisition module, a data processing module, and a data storage and output module.

[0007] Unmanned aerial vehicle (UAV) platform: adopts multi-rotor or fixed-wing models, has high-precision flight control capabilities, and can set flight paths, altitudes and speeds according to missions;

[0008] Payload integration module: includes a multispectral camera, a hyperspectral camera, and an upward-mounted high-precision spectral irradiance sensor, used to simultaneously acquire ground object images and downlink solar irradiance; it is fixed to the UAV by a shock-absorbing bracket to ensure imaging stability and time synchronization;

[0009] The data acquisition module consists of a field-programmable gate array (FPGA) synchronous control unit and a gigabit Ethernet / wireless transmission unit, ensuring that the timestamps of the image and irradiance data are consistent and transmitted to the data processing module in real time.

[0010] Data processing module: Based on an embedded processor or industrial computer, it runs PIF automatic identification, cross-sensor radiation normalization and absolute reflectance inversion algorithms to achieve real-time data processing;

[0011] Data storage and output module: Includes solid-state drive and USB 3.0, HDMI and network interfaces, used to store raw and calibration data, and supports local export and real-time display.

[0012] Based on the above system, this invention also proposes an automatic radiometric calibration method for UAV-borne multi-sensor systems, comprising the following steps:

[0013] Step S1: System Initialization and Parameter Configuration

[0014] Before the drone takes off, the entire calibration system is initialized and configured, specifically including:

[0015] Before the drone takes off, the calibration system is initialized, including the following configurations:

[0016] Sensor parameters: Input the intrinsic and extrinsic parameters and spectral response function of the multispectral and hyperspectral camera; input the sensitivity coefficient of the irradiance sensor. and dark current ;

[0017] Flight parameters: Set the flight altitude according to the operating area and resolution requirements. Forward overlap 80%, Lateral overlap 70%, Flight speed and image acquisition frequency ;

[0018] Algorithm parameters: Set the spectral stability threshold for PIF recognition. Spatial consistency threshold; Number of iterations for robust regression With convergence threshold and atmospheric correction parameters;

[0019] Step S2: Synchronous Data Acquisition

[0020] After the drone takes off according to the preset route, the synchronization control unit generates a cycle of... The synchronous trigger signals are sent to the main imaging sensor group and the irradiance measurement unit respectively, controlling both to synchronously start data acquisition:

[0021] The multispectral camera and hyperspectral camera in the main imaging sensor group simultaneously acquire images of ground features, obtaining multispectral image data. and hyperspectral image data ,in For pixel coordinates, The characteristic wavelength band of the multispectral camera, For the continuous band wavelengths of the hyperspectral camera, For the number of bands in a multispectral camera, This represents the number of bands in the hyperspectral camera.

[0022] The irradiance measurement unit simultaneously measures the downlink total solar radiation flux density to acquire irradiance data. ,in To collect timestamps, ensuring that each frame of image data corresponds to a unique irradiance data point,

[0023] After the data acquisition module adds metadata information such as timestamps and sensor numbers to the acquired image data and irradiance data, it transmits the data to the data processing module in real time through the data transmission unit.

[0024] Step S3: Automatic PIF Feature Identification and Algorithm Design

[0025] The PIF automatic identification module in the data processing module processes the synchronously acquired multispectral and hyperspectral images and automatically extracts spectrally stable PIF ground features;

[0026] Sub-step S3.1: Image preprocessing

[0027] The specific method is as follows:

[0028] S3.1.1 Radiation pretreatment: Eliminating system noise

[0029] (1) Dark current correction is performed using a "dark field calibration + band-by-band subtraction" strategy. The mathematical model is as follows:

[0030]

[0031] in:

[0032] Pixel positions in the original image In the band The numerical value;

[0033] The reference value for dark current in this band was obtained by capturing 50 dark field images with the lens blocked and taking the average value.

[0034] (2) Dead pixel repair

[0035] Repair is performed using a 3×3 neighborhood median filtering method: for each pixel, if its DN value deviates from the neighborhood mean by more than 3 times the standard deviation, it is judged as a bad pixel;

[0036] Replace the pixel value with the median value of the neighborhood center.

[0037] S3.1.2 Geometric preprocessing: Eliminating spatial distortion

[0038] Orthorectification is performed using a rational polynomial coefficient RPC model to convert image pixel coordinates into ground geographic coordinates, achieving a consistent representation of ground feature locations.

[0039] Corrected pixels Corresponding ground coordinates The calculation formula is as follows:

[0040]

[0041] in:

[0042] Pixel coordinates of the image before correction

[0043] : Corresponding ground projection coordinates

[0044] RPC model coefficients, obtained from camera calibration experiments;

[0045] Sub-step S3.2: Adaptive multi-scale image segmentation - generating a homogeneous region candidate set

[0046] A region-level recognition strategy based on "spectrally similar homogeneous regions" is proposed. A spectral-spatial collaborative adaptive mean-shift segmentation algorithm is designed. By pre-segmenting the image into homogeneous regions with consistent spectral features, subsequent recognition is carried out on a region-by-region basis, thereby improving the robustness and recognition accuracy of the system in complex environments.

[0047] S3.2.1 Data Input

[0048] Input data: spectral image

[0049] Image height, Image width, Number of bands;

[0050] Preset parameters: empirical coefficient

[0051] By band number Sure: Take 0.8 at that time. Time complexity is 1.2, spectral complexity coefficient Space kernel bandwidth coefficient Region constraint threshold, number of pixels Space compactness Spectral homogeneity ;

[0052] S3.2.2 Algorithm Flow

[0053] Step 1: Feature Vector Construction

[0054] For each pixel in the spectral image Construct a joint spectral-spatial feature vector:

[0055]

[0056] Pixels of 3D spectral vector;

[0057] Pixels The spatial coordinate vector;

[0058] Step 2: Calculation of local spectral complexity

[0059] For each pixel Take it Neighborhood Calculate the local spectral complexity:

[0060]

[0061] in For the neighboring region The variance of a band reflects the heterogeneity of the spectrum surrounding that pixel;

[0062] Step 3: Adaptive calculation of kernel bandwidth

[0063] Decompose the kernel bandwidth into spectral kernel bandwidth. With space kernel bandwidth Calculate separately:

[0064]

[0065] The value decreases as the local spectral complexity increases, adapting to the spectral characteristics of ground features in the spectral image.

[0066] It is positively correlated with the area of ​​the neighborhood, ensuring the spatial neighborhood continuity of the pixel;

[0067] Step 4: Iteration of mean drift in spectral-spatial dual-core weighted system

[0068] For each pixel , by Starting from the initial point, perform mean-shift iteration:

[0069] ① Determine the local region: Filter to meet the requirements and pixels , forming a local region ;

[0070] ② Calculate the dual-core weighted mean:

[0071]

[0072] Where: weighting coefficient ; spectral core Space Core ;

[0073] ③ Iterative update: Update the initial point to Repeat steps 1-2 until... ;

[0074] ④ Clustering and merging: Divide pixels that converge to the same mean into the same region. ;

[0075] Step 5: Region Constraint Filtering

[0076] For the segmented regions Regions that meet the following criteria are selected as the final homogeneous region candidate set:

[0077] ①Pixel count constraint:

[0078] ② Space compactness constraints: (Area Perimeter is the number of pixels contained in the region. (Number of pixels at the edge of the region)

[0079] ③ Spectral homogeneity constraint: (Std The standard deviation of the spectral vectors within the region, max (This is the global maximum value of the spectral vector of the spectral image);

[0080] S3.2.3 Algorithm Output

[0081] Output: Candidate set of homogeneous regions for this spectral image K represents the total number of regions, and each region satisfies the spectral-spatial homogeneity constraint.

[0082] Sub-step S3.3: Spectral stability analysis - dual-spectral index verification

[0083] For each homogeneous region The spectral feature vectors of the identified PIF ground features are extracted from multispectral and hyperspectral images. The spectral stability is evaluated by a two-dimensional index of "geometric similarity + statistical difference" to ensure that the spectral features of the identified PIF ground features remain stable under different time periods and different sensors.

[0084] S3.3.1 Calculation of Spectral Angle Matching Degree (SAM)

[0085] For the region At adjacent data collection times and average spectral vector and The SAM calculation formula is as follows:

[0086]

[0087] in:

[0088] Spectral angle

[0089] : region At any moment The average spectral reflectance or DN value, For wavelength,

[0090] The set of common bands for multispectral and hyperspectral imaging, and the range of common bands.

[0091] Spectral angle thresholds were set based on experimental experience data. ,

[0092] when At that time, it is assumed that the spectral characteristics of the region at the two moments are geometrically similar;

[0093] S3.3.2 Calculation of Spectral Information Divergence (SID)

[0094] set up , They are respectively regions The normalized spectral vectors at two time points, SID, are defined as follows:

[0095]

[0096] in The Kullback-Leibler divergence (KL divergence) is calculated using the following formula:

[0097]

[0098] in:

[0099] Extremely small positive numbers to prevent division by zero errors.

[0100] It needs to be pre-normalized to a probability distribution.

[0101] Set SID experience threshold ,

[0102] When SID At that time, it was assumed that the regional spectral distribution had statistical stability.

[0103] Regions that simultaneously meet the SAM and SID threshold conditions are marked as "spectrally stable candidate regions" and enter the subsequent spatial consistency verification stage.

[0104] Sub-step S3.4: Spatial consistency verification - removing spurious stable regions

[0105] S3.4.1 Spatial Heterogeneity Analysis

[0106] Computational area Spatial heterogeneity index This reflects the dispersion of pixel values ​​within a region, and the formula is as follows:

[0107]

[0108] in:

[0109] :area The standard deviation of the DN values ​​of all pixels within the range.

[0110] :area The mean of the DN values ​​of all pixels within the range.

[0111] It is a dimensionless index used to measure spectral uniformity within a region.

[0112] Specify spatial heterogeneity empirical threshold ,when When the pixel features within the region are uniform, the spatial consistency is good.

[0113] S3.4.2 Neighborhood Continuity Analysis

[0114] Statistical area The percentage of "spectrally stable candidate regions" in the 8-neighborhood The formula is as follows:

[0115]

[0116] in:

[0117] : 8. The number of spectrally stable candidate regions in the neighborhood, with a value range of : ,

[0118] The percentage of stable regions within the neighborhood is normalized to the [0,1] interval, and an empirical threshold for neighborhood consistency is set. ,when This indicates that the region is spatially connected to surrounding stable regions, thus excluding isolated, falsely stable regions.

[0119] Regions that simultaneously meet the following two conditions are ultimately identified as PIF (Picture-in-Flight) feature regions, and the DN (Domain Number) values ​​of all their pixels are extracted as a reference sample set for subsequent cross-sensor radiometric normalization:

[0120] ;

[0121] Step S4: Cross-sensor radiation normalization

[0122] Sub-step S4.1: Band matching

[0123] Band matching is performed using the convolution integral method of spectral response functions, and the fusion formula is as follows:

[0124]

[0125] in:

[0126] Multispectral bands after fusion The digital quantization value DN,

[0127] Hyperspectral images at wavelength The DN value at that location,

[0128] Multispectral camera Spectral response function SRM of the band,

[0129] Multispectral bands The wavelength range boundary;

[0130] Sub-step S4.2: Outlier Removal

[0131] S4.2.1 Initial Residual Calculation

[0132] Multispectral DN values ​​within the PIF region Matched hyperspectral DN values The difference is used as the initial residual:

[0133]

[0134] in:

[0135] : No. The residual of each pixel

[0136] Multispectral cameras in band The pixel value,

[0137] DN values ​​corresponding to hyperspectral data after band matching;

[0138] S4.2.2 Weight Calculation

[0139] The adaptive weight for each pixel is calculated based on the median absolute deviation (MAD) of the residuals, using the following weighting function:

[0140]

[0141] in:

[0142] The median absolute deviation (MAD) of the residuals.

[0143] 1.4826: The conversion factor between MAD and standard deviation under normal distribution (making MAD ≈ σ)

[0144] : No. The weight of each pixel is used to suppress the influence of outliers;

[0145] S4.2.3 Iterative Optimization

[0146] The model parameters are updated with weights, and new residuals are calculated based on the new parameters. This weight-parameter update process is repeated until the preset number of iterations is reached. Or the change in weights is less than the convergence threshold. ;

[0147] S4.2.4 Outlier Deletion

[0148] Remove weights The pixels, of which The weighted empirical threshold is used as the basis for determining the remaining pixels, which form the purified PIF sample set for subsequent radiometric normalization modeling.

[0149] Sub-step S4.3: Multispectral-Hyperspectral Adaptive Weighted Robust Regressive Radiometric Normalization Algorithm

[0150] S4.3.1 Algorithm Input

[0151] Input data:

[0152] ① Multispectral images in band DN value matrix ;

[0153] ② The DN value matrix of the hyperspectral image after band matching;

[0154] ③ Cleaned PIF sample set ( (sample size)

[0155] in, Multispectral cameras in spectral bands The DN value; : DN values ​​corresponding to hyperspectral data after band matching;

[0156] Preset parameters: Weight iteration convergence threshold Maximum number of iterations ;

[0157] S4.3.2 Algorithm Flow

[0158] Step 1: PIF sample weight initialization, for each sample in the PIF sample set Initialize weights ( );

[0159] Step 2: Adaptive weighted robust regression modeling for multispectral bands Iterative solution of the linear transformation model parameters :

[0160] ① Parameter solution for the t-th iteration: Construct the objective function with the goal of minimizing the weighted sum of squared residuals:

[0161]

[0162] right Taking the partial derivative and setting it to 0, we obtain the analytical solution:

[0163]

[0164] ,

[0165] ② Weight update: Calculate the residual of the t-th iteration. Update weights based on robust residual estimation:

[0166]

[0167]

[0168] in ,

[0169] ③ Convergence judgment: If or Stop iteration, take , Otherwise, repeat steps 1-2.

[0170] Step 3: Cross-band radiative normalization is performed on all pixels in the hyperspectral image. The result obtained by solving , Convert hyperspectral DN values ​​to multispectral reference DN values:

[0171]

[0172] Step 4: Model accuracy verification

[0173] Calculate the normalized residuals of the PIF sample set. Verify the following metrics:

[0174] residual mean (Unbiasedness);

[0175] residual standard deviation (Accuracy meets standards);

[0176] S4.3.3 Algorithm Output

[0177] Output result:

[0178] ① Multispectral band Corresponding gain coefficient Offset coefficient ;

[0179] ② Normalized hyperspectral DN value matrix;

[0180] This achieves radiometric consistency of multi-source remote sensing data, supporting subsequent quantitative analysis;

[0181] Step S5: Absolute reflectivity conversion - Real-time conversion algorithm based on coupled online irradiance

[0182] Even after cross-sensor normalization, the radiation data is still a relative value. Therefore, a "real-time absolute reflectance conversion algorithm coupled with online irradiance" is required. This algorithm uses the irradiance data collected in real time by the UAV's onboard sensor as a dynamic benchmark, and combines a simplified atmospheric correction model to eliminate environmental interference from light intensity fluctuations and atmospheric transmission attenuation. This algorithm converts the normalized DN value into a standardized surface absolute reflectance.

[0183] Sub-step S5.1 Algorithm preparation (input data and core assumptions)

[0184] S5.1.1 Input Data Structure

[0185] Radiation fundamental data: Normalized DN value dataset across sensors ,in For pixel coordinates, For the spectral bands, this data has eliminated differences in multi-sensor responses;

[0186] Dynamic irradiance data: Downlink total solar irradiance flux density synchronously acquired by the airborne irradiance sensor. , To collect timestamps, ensure a one-to-one correspondence with each frame of DN data;

[0187] Auxiliary parameter data: UAV GPS / IMU data, sensor laboratory calibration parameters, and real-time meteorological parameters;

[0188] S5.1.2 Conditional Assumptions

[0189] Based on the theory of remote sensing radiative transfer, the following assumptions are made:

[0190] Assumption 1: The sensor output DN value is linearly related to the input radiance, that is... This assumption forms the physical basis for radiance conversion;

[0191] Assumption 2: When the UAV flies at low altitude, the aerosol is vertically uniformly distributed and has a low concentration, and its scattering effect can be approximated by a simplified model;

[0192] Assumption 3: Downward solar radiation can be considered as uniformly distributed at the instant of observation. This assumption is guaranteed through sensor installation location optimization and data verification.

[0193] Sub-step S5.2 Dynamic correction of irradiance data

[0194] To obtain high-precision downlink solar irradiance flux density, multi-stage dynamic correction of the raw data is necessary, which is implemented as follows:

[0195] S5.2.1 Dark Current Dynamic Correction

[0196] The sensor still exhibits a weak output, known as dark current, even in the absence of light. This dark current needs to be subtracted in real time to eliminate baseline drift.

[0197] Correction strategy: The instantaneous dark current is calculated using a "sliding window mean" strategy, and the mathematical model is as follows:

[0198]

[0199] in:

[0200] : Moment The corrected dark current value;

[0201] The sliding window length is set to 5 frames to balance response speed and stability.

[0202] : No. Dark field data of the frame;

[0203] S5.2.2 Temperature Drift Compensation

[0204] The sensitivity of the sensor drifts with changes in operating temperature and needs to be compensated for using a temperature coefficient.

[0205] Compensation model: Linear temperature response model, temperature correction coefficient Described using a linear model:

[0206]

[0207] in:

[0208] Standard temperature Sensitivity coefficient at the following levels

[0209] : band Temperature coefficient (unit: The temperature was calibrated through high and low temperature environmental tests (-10℃ to 50℃).

[0210] : Moment The measured operating temperature of the sensor (acquired by the built-in temperature sensor);

[0211] S5.2.3 Random Noise Filtering

[0212] A 5-point moving average filter was used to remove Gaussian random noise. The filtered irradiance data The calculation is as follows:

[0213] S5.2.4 Final Irradiance Determination

[0214] Based on the above correction steps, the time is obtained. band True downlink solar irradiance The calculation formula is as follows:

[0215]

[0216] in:

[0217] : The original irradiance data after filtering;

[0218] Dark current after sliding window correction;

[0219] Temperature compensation coefficient;

[0220] Sub-step S5.3: Calculation of the radiance of the top of the atmosphere

[0221] Based on the normalized DN value and combined with the camera's laboratory calibration parameters, the radiance of ground features at the top of the atmosphere is calculated. The preliminary relationship of the sensor response characteristics is established, and the calculation formula is as follows:

[0222]

[0223] in:

[0224] Pixels At wavelength Atmospheric top radiance at that location.

[0225] DN is the digitally quantized value after cross-sensor normalization, i.e., the data after radiometric normalization.

[0226] Camera at wavelength The gain coefficient at a given point reflects the sensor's sensitivity.

[0227] Offset coefficient, reflecting the zero-point deviation of the system;

[0228] Sub-step S5.4: Atmospheric correction and absolute reflectance calculation

[0229] Based on simplification The model undergoes atmospheric correction. Considering the characteristics of low-altitude flight of UAVs, the influence of complex aerosol scattering is ignored, and only atmospheric molecular scattering and water vapor absorption are considered to calculate the absolute surface reflectance. The formula is as follows:

[0230]

[0231] in:

[0232] : Absolute surface reflectance

[0233] Atmospheric radiance

[0234] Atmospheric path radiation

[0235] Earth-Sun distance correction factor

[0236] Corrected solar irradiance

[0237] : Solar zenith angle

[0238] Atmospheric transmittance

[0239] S5.4.1 Atmospheric path radiation

[0240] Through dark-field experiments or by selecting deep water bodies and shadowy dark areas in images... Estimate using the following formula:

[0241]

[0242] S5.4.2 Earth-Sun Distance Correction Factor

[0243] Since Earth's orbit around the sun is elliptical, it is necessary to correct for the effect of changes in the Earth-Sun distance on solar radiation. The calculation formula is as follows:

[0244]

[0245] in:

[0246] DOY: Yearly Accumulation Days, Value Range Or 366,

[0247] The relative proportion of the Earth-Sun distance to the average distance.

[0248] S5.4.3 Solar Zenith Angle

[0249] Calculated based on the real-time geographic location and data collection time of the drone: First, calculate the solar altitude angle. :

[0250]

[0251] Then the solar zenith angle is:

[0252]

[0253] in:

[0254] Latitude of the observation point

[0255] The solar declination angle is calculated using the following formula:

[0256] ,

[0257] ω: Hour angle, calculated using the following formula:

[0258] ,

[0259] Collection time,

[0260] S5.4.4 Atmospheric transmittance

[0261] Divided into atmospheric molecular transmittance and water vapor transmission rate :

[0262]

[0263] in:

[0264] Atmospheric optical thickness, obtained from actual measurements at local meteorological stations or retrieved through hyperspectral data inversion.

[0265] Molecular scattering term, calculated using the Rayleigh scattering model:

[0266] ,

[0267] Rayleigh scattering coefficient, relative to wavelength Related,

[0268] Atmospheric water vapor content (unit: g / cm²) can be obtained through water vapor absorption spectrum inversion or by querying the US Standard Atmosphere database.

[0269] The water vapor absorption term is obtained from empirical models or by looking up tables.

[0270] Compared with the prior art, the present invention has the following significant advantages:

[0271] 1. No ground calibration target is required; fully autonomous field radiation calibration is achieved through a combination of PIF automatic identification and online irradiance measurement.

[0272] 2. Balancing relative consistency and absolute accuracy, it eliminates response differences between multiple sensors through PIF and utilizes airborne irradiance to achieve absolute reflectivity inversion.

[0273] 3. Strong environmental adaptability: Through multi-feature fusion PIF recognition and dynamic irradiance correction, it effectively suppresses interference such as sudden changes in illumination and attitude disturbances.

[0274] 4. The entire process is automated, from data acquisition to reflectivity output, which can be completed in real time on the embedded platform, meeting the timeliness requirements of emergency remote sensing. Detailed Implementation

[0275] The UAV-borne multi-sensor automatic radiometric calibration system includes a UAV platform, a payload integration module, a data acquisition module, a data processing module, and a data storage and output module.

[0276] Unmanned aerial vehicle (UAV) platform: adopts multi-rotor or fixed-wing models, has high-precision flight control capabilities, and can set flight paths, altitudes and speeds according to missions;

[0277] Payload integration module: includes a multispectral camera, a hyperspectral camera, and an upward-mounted high-precision spectral irradiance sensor, used to simultaneously acquire ground object images and downlink solar irradiance; it is fixed to the UAV by a shock-absorbing bracket to ensure imaging stability and time synchronization;

[0278] The data acquisition module consists of a field-programmable gate array (FPGA) synchronous control unit and a gigabit Ethernet / wireless transmission unit, ensuring that the timestamps of the image and irradiance data are consistent and transmitted to the data processing module in real time.

[0279] Data processing module: Based on an embedded processor or industrial computer, it runs PIF automatic identification, cross-sensor radiation normalization and absolute reflectance inversion algorithms to achieve real-time data processing;

[0280] Data storage and output module: Includes solid-state drive and USB 3.0, HDMI and network interfaces, used to store raw and calibration data, and supports local export and real-time display.

[0281] Based on the above system, this invention also proposes an automatic radiometric calibration method for UAV-borne multi-sensor systems, comprising the following steps:

[0282] Step S1: System Initialization and Parameter Configuration

[0283] Before the drone takes off, the entire calibration system is initialized and configured, specifically including:

[0284] Before the drone takes off, the calibration system is initialized, including the following configurations:

[0285] Sensor parameters: Intrinsic parameters (focal length) of the input multispectral and hyperspectral cameras Pixel size Main point ), External parameters (initial roll angle) Pitch angle Yaw angle ) and spectral response function; sensitivity coefficient of input irradiance sensor and dark current ;

[0286] Flight parameters: Set the flight altitude according to the operating area and resolution requirements. Forward overlap 80%, Lateral overlap 70%, Flight speed and image acquisition frequency ;

[0287] Algorithm parameters: Set the spectral stability threshold for PIF recognition. Spatial consistency threshold; Number of iterations for robust regression With convergence threshold ; and atmospheric correction parameters (atmospheric optical thickness) Water vapor content ;

[0288] Step S2: Synchronous Data Acquisition

[0289] After the drone takes off according to the preset route, the synchronization control unit generates a cycle of... The synchronous trigger signals are sent to the main imaging sensor group and the irradiance measurement unit respectively, controlling both to synchronously start data acquisition:

[0290] The multispectral camera and hyperspectral camera in the main imaging sensor group simultaneously acquire images of ground features, obtaining multispectral image data. and hyperspectral image data ,in For pixel coordinates, The characteristic wavelength band of the multispectral camera, For the continuous band wavelengths of the hyperspectral camera, For the number of bands in a multispectral camera, This represents the number of bands in the hyperspectral camera.

[0291] The irradiance measurement unit simultaneously measures the downlink total solar radiation flux density to acquire irradiance data. ,in To collect timestamps, ensuring that each frame of image data corresponds to a unique irradiance data point,

[0292] After the data acquisition module adds metadata information such as timestamps and sensor numbers to the acquired image data and irradiance data, it transmits the data to the data processing module in real time through the data transmission unit.

[0293] Step S3: Automatic PIF Feature Identification and Algorithm Design

[0294] The PIF automatic identification module in the data processing module processes the synchronously acquired multispectral and hyperspectral images and automatically extracts spectrally stable PIF features. This step designs a "multi-feature fusion PIF automatic identification algorithm" and constructs a three-level progressive process of "multi-scale segmentation - dual-spectral verification - spatial filtering" to achieve fully automatic and high-precision identification.

[0295] Sub-step S3.1: Image preprocessing

[0296] The accuracy of PIF recognition is highly dependent on the quality of the input image. To ensure the stability and accuracy of subsequent region segmentation and feature extraction, it is essential to effectively eliminate the two main interfering factors, "system noise" and "geometric distortion," during the preprocessing stage. Specific methods are as follows:

[0297] S3.1.1 Radiation pretreatment: Eliminating system noise

[0298] During the imaging process, the sensor introduces two typical types of system noise: dark current and dead pixels. Both of these affect the accuracy of the radiation value and need to be removed through correction.

[0299] (1) Dark current correction: Dark current is the weak output signal of the sensor under no-light conditions, which will cause the baseline of the DN value to shift. The correction is performed using the strategy of "dark field calibration + band-by-band subtraction", and its mathematical model is as follows:

[0300]

[0301] in:

[0302] Pixel positions in the original image In the band The numerical value;

[0303] The reference value for dark current in this band was obtained by capturing 50 dark field images with the lens blocked and taking the average value.

[0304] This model is built upon a physical mechanism: the true radiative response is composed of the superposition of the "ground object radiative response" and the "dark current". By subtracting the dark current, the true ground object radiative response can be reconstructed.

[0305] (2) Dead pixel repair

[0306] Defective pixels are isolated noise points caused by abnormal pixel response. A 3×3 neighborhood median filtering method is used for repair: for each pixel, if its DN value deviates from the neighborhood mean by more than 3 times the standard deviation, it is judged as a defective pixel.

[0307] Replace the pixel value with the median value of the neighborhood center.

[0308] S3.1.2 Geometric preprocessing: Eliminating spatial distortion

[0309] Attitude fluctuations during UAV flight (such as roll and pitch) can cause geometric distortion in images, resulting in shifts in the position of the same ground feature across different images. This severely impacts the accuracy of spectral consistency analysis and subsequent region segmentation. An orthorectification model with rational polynomial coefficients is employed to convert image pixel coordinates into ground geographic coordinates, achieving a consistent representation of ground feature locations.

[0310] Corrected pixels Corresponding ground coordinates The calculation formula is as follows:

[0311]

[0312] in:

[0313] Pixel coordinates of the image before correction

[0314] : Corresponding ground projection coordinates

[0315] RPC model coefficients, obtained from camera calibration experiments;

[0316] Sub-step S3.2: Adaptive multi-scale image segmentation - generating a homogeneous region candidate set

[0317] A region-level recognition strategy based on "spectrally similar homogeneous regions" is proposed. A spectral-spatial collaborative adaptive mean-shift segmentation algorithm is designed. By pre-segmenting the image into homogeneous regions with consistent spectral features, subsequent recognition is carried out on a region-by-region basis, thereby improving the robustness and recognition accuracy of the system in complex environments.

[0318] S3.2.1 Data Input

[0319] Input data: spectral image

[0320] Image height, Image width, Number of bands;

[0321] Preset parameters: empirical coefficient

[0322] By band number Sure: Take 0.8 at that time. Time complexity is 1.2, spectral complexity coefficient Space kernel bandwidth coefficient Region constraint threshold, number of pixels Space compactness Spectral homogeneity ;

[0323] S3.2.2 Algorithm Flow

[0324] Step 1: Feature Vector Construction

[0325] For each pixel in the spectral image Construct a joint spectral-spatial feature vector:

[0326]

[0327] Pixels of 3D spectral vector;

[0328] Pixels The spatial coordinate vector;

[0329] Step 2: Calculation of local spectral complexity

[0330] For each pixel Take it Neighborhood Calculate the local spectral complexity:

[0331]

[0332] in For the neighboring region The variance of a band reflects the heterogeneity of the spectrum surrounding that pixel;

[0333] Step 3: Adaptive calculation of kernel bandwidth

[0334] Decompose the kernel bandwidth into spectral kernel bandwidth. With space kernel bandwidth Calculate separately:

[0335]

[0336] The value decreases as the local spectral complexity increases, adapting to the spectral characteristics of ground features in the spectral image.

[0337] It is positively correlated with the area of ​​the neighborhood, ensuring the spatial neighborhood continuity of the pixel;

[0338] Step 4: Iteration of mean drift in spectral-spatial dual-core weighted system

[0339] For each pixel , by Starting from the initial point, perform mean-shift iteration:

[0340] ① Determine the local region: Filter to meet the requirements and pixels , forming a local region ;

[0341] ② Calculate the dual-core weighted mean:

[0342]

[0343] Where: weighting coefficient (Prioritize spectral similarity); spectral kernel (Gaussian kernel, weights balancing spectral similarity); Spatial kernel (Epanechnikov kernel, enhancing the locality of spatial neighborhood);

[0344] ③ Iterative update: Update the initial point to Repeat steps 1-2 until... (Convergence threshold);

[0345] ④ Clustering and merging: Divide pixels that converge to the same mean into the same region. ;

[0346] Step 5: Region Constraint Filtering

[0347] For the segmented regions Regions that meet the following criteria are selected as the final homogeneous region candidate set:

[0348] ①Pixel count constraint:

[0349] ② Space compactness constraints: (Area Perimeter is the number of pixels contained in the region. (Number of pixels at the edge of the region)

[0350] ③ Spectral homogeneity constraint: (Std The standard deviation of the spectral vectors within the region, max (This is the global maximum value of the spectral vector of the spectral image);

[0351] S3.2.3 Algorithm Output

[0352] Output: Candidate set of homogeneous regions for this spectral image K represents the total number of regions, and each region satisfies the spectral-spatial homogeneity constraint.

[0353] Sub-step S3.3: Spectral stability analysis - dual-spectral index verification

[0354] For each homogeneous region The spectral feature vectors of the identified PIF ground features are extracted from multispectral and hyperspectral images. The spectral stability is evaluated by a two-dimensional index of "geometric similarity + statistical difference" to ensure that the spectral features of the identified PIF ground features remain stable under different time periods and different sensors.

[0355] S3.3.1 Calculation of Spectral Angle Matching Degree (SAM)

[0356] Spectral angle This reflects the directional similarity of two spectral vectors; a smaller angle indicates more stable spectral characteristics and less susceptibility to changes in light intensity. For a region... At adjacent data collection times and average spectral vector and The SAM calculation formula is as follows:

[0357]

[0358] in:

[0359] Spectral angle (unit: radians rad)

[0360] : region At any moment The average spectral reflectance or DN value, For wavelength,

[0361] The set of common bands for multispectral and hyperspectral imaging, and the range of common bands.

[0362] Spectral angle thresholds were set based on experimental experience data. ,

[0363] when At that time, it is assumed that the spectral characteristics of the region at the two moments are geometrically similar;

[0364] S3.3.2 Calculation of Spectral Information Divergence (SID)

[0365] SID measures the difference in spectral distribution from an information theory perspective, compensating for the deficiency of SAM, which only considers direction and ignores amplitude distribution. Let... , They are respectively regions The normalized spectral vectors at two time points, SID, are defined as follows:

[0366]

[0367] in The Kullback-Leibler divergence (KL divergence) is calculated using the following formula:

[0368]

[0369] in:

[0370] Extremely small positive numbers to prevent division by zero errors.

[0371] It needs to be pre-normalized to a probability distribution.

[0372] Set SID experience threshold ,

[0373] When SID At that time, it was assumed that the regional spectral distribution had statistical stability.

[0374] Regions that simultaneously meet the SAM and SID threshold conditions are marked as "spectrally stable candidate regions" and enter the subsequent spatial consistency verification stage.

[0375] Sub-step S3.4: Spatial consistency verification - removing spurious stable regions

[0376] False stable regions may exist in the spectrally stable candidate regions due to image noise, shadow edges, etc., and further screening is required through spatial feature analysis to ensure the spatial continuity and consistency of PIF features.

[0377] S3.4.1 Spatial Heterogeneity Analysis

[0378] Computational area Spatial heterogeneity index This reflects the dispersion of pixel values ​​within a region, and the formula is as follows:

[0379]

[0380] in:

[0381] :area The standard deviation of the DN values ​​of all pixels within the range.

[0382] :area The mean of the DN values ​​of all pixels within the range.

[0383] It is a dimensionless index used to measure spectral uniformity within a region.

[0384] Specify spatial heterogeneity empirical threshold ,when When the pixel features within the region are uniform, the spatial consistency is good.

[0385] S3.4.2 Neighborhood Continuity Analysis

[0386] Statistical area The percentage of "spectrally stable candidate regions" in the 8-neighborhood The formula is as follows:

[0387]

[0388] in:

[0389] : 8. The number of spectrally stable candidate regions in the neighborhood, with a value range of : ,

[0390] The percentage of stable regions within the neighborhood is normalized to the [0,1] interval, and an empirical threshold for neighborhood consistency is set. ,when This indicates that the region is spatially connected to surrounding stable regions, thus excluding isolated, falsely stable regions.

[0391] Regions that simultaneously meet the following two conditions are ultimately identified as PIF (Picture-in-Flight) feature regions, and the DN (Domain Number) values ​​of all their pixels are extracted as a reference sample set for subsequent cross-sensor radiometric normalization:

[0392] ;

[0393] Step S4: Cross-sensor radiation normalization

[0394] The core of cross-sensor radiometric normalization is to establish the radiometric response conversion relationship between multispectral and hyperspectral sensors, normalizing the output data of different sensors to a unified radiometric reference and eliminating radiometric bias caused by differences in sensor responses. This step uses the DN value of the PIF land cover area as a reference and is achieved through three sub-steps: band matching, outlier removal, and robust regression modeling.

[0395] Sub-step S4.1: Band matching

[0396] Multispectral cameras use discrete feature bands, while hyperspectral cameras use continuous spectral bands. The hyperspectral data must first be fused into feature bands corresponding to the multispectral data to ensure comparability in the spectral dimension. Band matching is performed using the spectral response function convolution integral method, and the fusion formula is as follows:

[0397]

[0398] in:

[0399] Multispectral bands after fusion The digital quantization value DN,

[0400] Hyperspectral images at wavelength The DN value at that location,

[0401] Multispectral camera The spectral response function (SRM) for a given band.

[0402] Multispectral bands The wavelength range boundary;

[0403] Integration can be achieved using the response function obtained through laboratory calibration, or it can be approximated by numerical discretization and summation.

[0404] Sub-step S4.2: Outlier Removal

[0405] Within a PIF (Peripherally Identifiable Information Foundation) feature area, there may be a small number of unstable feature pixels. These outliers can affect the accuracy of the radiative response model. An iterative weighted least squares algorithm is designed to remove them.

[0406] S4.2.1 Initial Residual Calculation

[0407] Multispectral DN values ​​within the PIF region Matched hyperspectral DN values The difference is used as the initial residual:

[0408]

[0409] in:

[0410] : No. The residual of each pixel

[0411] Multispectral cameras in band The pixel value,

[0412] DN values ​​corresponding to hyperspectral data after band matching;

[0413] S4.2.2 Weight Calculation

[0414] The adaptive weight for each pixel is calculated based on the median absolute deviation (MAD) of the residuals, using the following weighting function:

[0415]

[0416] in:

[0417] The median absolute deviation (MAD) of the residuals.

[0418] 1.4826: The conversion factor between MAD and standard deviation under normal distribution (making MAD ≈ σ)

[0419] : No. The weight of each pixel is used to suppress the influence of outliers;

[0420] S4.2.3 Iterative Optimization

[0421] The model parameters are updated with weights, and new residuals are calculated based on the new parameters. This weight-parameter update process is repeated until the preset number of iterations is reached. Or the change in weights is less than the convergence threshold. ;

[0422] S4.2.4 Outlier Deletion

[0423] Remove weights The pixels, of which The weighted empirical threshold is used as the basis for determining the remaining pixels, which form the purified PIF sample set for subsequent radiometric normalization modeling.

[0424] Sub-step S4.3: Multispectral-Hyperspectral Adaptive Weighted Robust Regressive Radiometric Normalization Algorithm

[0425] S4.3.1 Algorithm Input

[0426] Input data:

[0427] ① Multispectral images in band DN value matrix ;

[0428] ②DN value matrix of hyperspectral image after band matching (with multispectral bands) correspond);

[0429] ③ Cleaned PIF sample set ( (Number of samples);

[0430] in, Multispectral cameras in spectral bands The DN value; : DN values ​​corresponding to hyperspectral data after band matching;

[0431] Preset parameters: Weight iteration convergence threshold Maximum number of iterations ;

[0432] S4.3.2 Algorithm Flow

[0433] Step 1: PIF sample weight initialization, for each sample in the PIF sample set Initialize weights ( );

[0434] Step 2: Adaptive weighted robust regression modeling for multispectral bands Iterative solution of the linear transformation model parameters :

[0435] ① Parameter solution for the t-th iteration: Construct the objective function with the goal of minimizing the weighted sum of squared residuals:

[0436]

[0437] right Taking the partial derivative and setting it to 0, we obtain the analytical solution:

[0438]

[0439] ,

[0440] ② Weight update: Calculate the residual of the t-th iteration. Weights are updated based on robust estimation of residuals (MAD: median absolute deviation):

[0441]

[0442]

[0443] in (Huber weighting function, balancing noise robustness and estimation efficiency).

[0444] ③ Convergence judgment: If or Stop iteration, take , Otherwise, repeat steps 1-2.

[0445] Step 3: Cross-band radiative normalization is performed on all pixels in the hyperspectral image. The result obtained by solving , Convert hyperspectral DN values ​​to multispectral reference DN values:

[0446]

[0447] Step 4: Model accuracy verification

[0448] Calculate the normalized residuals of the PIF sample set. Verify the following metrics:

[0449] residual mean (Unbiasedness);

[0450] residual standard deviation (Accuracy meets standards);

[0451] S4.3.3 Algorithm Output

[0452] Output result:

[0453] ① Multispectral band Corresponding gain coefficient Offset coefficient ;

[0454] ② Normalized hyperspectral DN value matrix (consistent with multispectral radiation reference).

[0455] This achieves radiometric consistency of multi-source remote sensing data, supporting subsequent quantitative analysis;

[0456] Step S5: Absolute reflectivity conversion - Real-time conversion algorithm based on coupled online irradiance

[0457] Even after cross-sensor normalization, the radiation data is still a relative value. Therefore, a "real-time absolute reflectance conversion algorithm coupled with online irradiance" is required. This algorithm uses the irradiance data collected in real time by the UAV's onboard sensor as a dynamic benchmark, and combines a simplified atmospheric correction model to eliminate environmental interference from light intensity fluctuations and atmospheric transmission attenuation. This algorithm converts the normalized DN value into a standardized surface absolute reflectance.

[0458] Sub-step S5.1 Algorithm preparation (input data and core assumptions)

[0459] S5.1.1 Input Data Structure

[0460] Radiation fundamental data: Normalized DN value dataset across sensors ,in For pixel coordinates, For the spectral bands, this data has eliminated differences in multi-sensor responses;

[0461] Dynamic irradiance data: Downlink total solar irradiance flux density synchronously acquired by the airborne irradiance sensor. , To collect timestamps, ensure a one-to-one correspondence with each frame of DN data;

[0462] Auxiliary parameter data: UAV GPS / IMU data (longitude Lon, latitude Lat, flight attitude), sensor laboratory calibration parameters (gain) Offset ) and real-time meteorological parameters (atmospheric optical thickness) Water vapor content );

[0463] S5.1.2 Conditional Assumptions

[0464] Based on the theory of remote sensing radiative transfer, the following assumptions are made:

[0465] Assumption 1: The sensor output DN value is linearly related to the input radiance (this has been experimentally verified), that is... This assumption forms the physical basis for radiance conversion;

[0466] Assumption 2: When the UAV flies at low altitude (≤1000m), the aerosol is vertically uniform and has a low concentration. Its scattering effect can be approximated by a simplified model (compared to high-altitude remote sensing, the contribution of low-altitude aerosols to radiative transfer is reduced by more than 60%).

[0467] Assumption 3: Downward solar radiation can be considered as uniformly distributed at the moment of observation (the single-point data collected by the irradiance sensor can represent the illumination conditions of the entire imaging area). This assumption is guaranteed by optimizing the sensor installation position (no obstruction on the top of the fuselage) and data verification (the PIF ground object spectral variation coefficient within the same frame image is ≤5%).

[0468] Sub-step S5.2 Dynamic correction of irradiance data

[0469] In actual operation, the raw output of an irradiance sensor is affected by a variety of factors, including:

[0470] Dark current: a weak output under conditions of no light.

[0471] Temperature drift: Sensitivity shift caused by changes in operating temperature;

[0472] Random noise: fluctuations introduced by electronic circuits and environmental interference.

[0473] To obtain high-precision downlink solar irradiance flux density, multi-stage dynamic correction of the raw data is necessary, which is implemented as follows:

[0474] S5.2.1 Dark Current Dynamic Correction

[0475] The sensor still exhibits a weak output, known as dark current, even in the absence of light. This dark current needs to be subtracted in real time to eliminate baseline drift.

[0476] Correction strategy: The instantaneous dark current is calculated using a "sliding window mean" strategy, and the mathematical model is as follows:

[0477]

[0478] in:

[0479] : Moment The corrected dark current value;

[0480] The sliding window length is set to 5 frames to balance response speed and stability.

[0481] : No. Dark field data of the frame (collected periodically through the sensor's built-in occlusion mechanism);

[0482] S5.2.2 Temperature Drift Compensation

[0483] The sensitivity of the sensor drifts with changes in operating temperature and needs to be compensated for using a temperature coefficient.

[0484] Compensation model: Linear temperature response model, temperature correction coefficient Described using a linear model:

[0485]

[0486] in:

[0487] Standard temperature Sensitivity coefficient at the following levels

[0488] : band Temperature coefficient (unit: The temperature was calibrated through high and low temperature environmental tests (-10℃ to 50℃).

[0489] : Moment The measured operating temperature of the sensor (acquired by the built-in temperature sensor);

[0490] S5.2.3 Random Noise Filtering

[0491] The raw data contains random noise, which affects the accuracy of irradiance estimation.

[0492] A 5-point moving average filter was used to remove Gaussian random noise. The filtered irradiance data The calculation is as follows:

[0493] ,

[0494] S5.2.4 Final Irradiance Determination

[0495] Based on the above correction steps, the time is obtained. band True downlink solar irradiance The calculation formula is as follows:

[0496]

[0497] in:

[0498] : The original irradiance data after filtering;

[0499] Dark current after sliding window correction;

[0500] Temperature compensation coefficient;

[0501] Sub-step S5.3: Calculation of the radiance of the top of the atmosphere

[0502] Based on the normalized DN value and combined with the camera's laboratory calibration parameters, the radiance of ground features at the top of the atmosphere is calculated. The preliminary relationship of the sensor response characteristics is established, and the calculation formula is as follows:

[0503]

[0504] in:

[0505] Pixels At wavelength Atmospheric radiance at the top of the atmosphere (unit: ),

[0506] DN is the digitally quantized value after cross-sensor normalization, i.e., the data after radiometric normalization.

[0507] Camera at wavelength The gain coefficient at a given point reflects the sensor's sensitivity.

[0508] Offset coefficient, reflecting the zero-point deviation of the system;

[0509] Sub-step S5.4: Atmospheric correction and absolute reflectance calculation

[0510] Based on simplification The model performs atmospheric correction for low-altitude flight of drones (typically below). Based on the characteristics of atmospheric molecular scattering and neglecting the complex scattering effects of aerosols (assuming low and uniform aerosol concentration), and considering only atmospheric molecular scattering and water vapor absorption, the absolute surface reflectance is calculated. The formula is as follows:

[0511]

[0512] in:

[0513] Surface absolute reflectance (unit: dimensionless, range) ),

[0514] Atmospheric top radiance (from sub-step S5.2).

[0515] Atmospheric path radiation (by) Model calculations reflect atmospheric scattering contributions.

[0516] Earth-Sun distance correction factor (unit: AU).

[0517] Corrected solar irradiance (from sub-step S5.1).

[0518] Solar zenith angle (unit: radians rad).

[0519] Atmospheric transmittance (including molecular scattering and water vapor absorption).

[0520] S5.4.1 Atmospheric path radiation

[0521] This refers to the amount of solar radiation that reaches the ground and directly enters the sensor after being scattered by the atmosphere. This is achieved through dark-field experiments or by selecting deep water bodies or shadowy areas in the image. (Its surface reflectance is approximately 0), estimated using the following formula:

[0522]

[0523] S5.4.2 Earth-Sun Distance Correction Factor

[0524] Since Earth's orbit around the sun is elliptical, it is necessary to correct for the effect of changes in the Earth-Sun distance on solar radiation. The calculation formula is as follows:

[0525]

[0526] in:

[0527] DOY: Yearly Accumulation Days, Value Range Or 366,

[0528] The relative proportion of the Earth-Sun distance to the average distance (unit: AU).

[0529] S5.4.3 Solar Zenith Angle

[0530] The solar altitude angle is the angle between sunlight and the local ground plane. It is calculated using the real-time geographic location of the drone (longitude Lon, latitude Lat) and the data collection time (UTC time): First, the solar altitude angle is calculated. :

[0531]

[0532] Then the solar zenith angle is:

[0533]

[0534] in:

[0535] Latitude of the observation point (unit: degrees).

[0536] The solar declination angle (unit: degrees) is calculated using the following formula:

[0537] ,

[0538] ω: Hour angle (unit: degrees), calculated using the following formula:

[0539] ,

[0540] Data acquisition time (UTC time, unit: hours).

[0541] S5.4.4 Atmospheric transmittance

[0542] This refers to the transmittance of solar radiation as it travels through the atmosphere from the top of the atmosphere to the ground and then back to the sensor. It is divided into atmospheric molecular transmittance. and water vapor transmission rate :

[0543]

[0544] in:

[0545] Atmospheric optical thickness, obtained from actual measurements at local meteorological stations or retrieved through hyperspectral data inversion.

[0546] Molecular scattering term, calculated using the Rayleigh scattering model:

[0547] ,

[0548] Rayleigh scattering coefficient, relative to wavelength Related,

[0549] Atmospheric water vapor content (unit: g / cm²) can be obtained through water vapor absorption spectrum inversion or by querying the US Standard Atmosphere database.

[0550] The water vapor absorption term is obtained from empirical models or by looking up tables.

Claims

1. An unmanned aerial vehicle (UAV) multi-sensor automatic radiometric calibration system, characterized in that: This includes an unmanned aerial vehicle (UAV) platform, payload integration module, data acquisition module, data processing module, and data storage and output module. Unmanned aerial vehicle (UAV) platform: adopts multi-rotor or fixed-wing models, has high-precision flight control capabilities, and can set flight paths, altitudes and speeds according to missions; Payload integration module: includes a multispectral camera, a hyperspectral camera, and an upward-mounted high-precision spectral irradiance sensor, used to simultaneously acquire ground object images and downlink solar irradiance; it is fixed to the UAV by a shock-absorbing bracket to ensure imaging stability and time synchronization; The data acquisition module consists of a field-programmable gate array (FPGA) synchronous control unit and a gigabit Ethernet / wireless transmission unit, ensuring that the timestamps of the image and irradiance data are consistent and transmitted to the data processing module in real time. Data processing module: Based on an embedded processor or industrial computer, it runs algorithms for automatic identification of pseudo-invariant feature PIF, cross-sensor radiation normalization, and absolute reflectance inversion to achieve real-time data processing; Data storage and output module: Includes solid-state drive and USB 3.0, HDMI and network interfaces, used to store raw and calibration data, and supports local export and real-time display.

2. Based on the above system, this invention proposes an automatic radiometric calibration method for unmanned aerial vehicles (UAVs) with multiple sensors, characterized in that: Includes the following steps: Step S1: System Initialization and Parameter Configuration Before the drone takes off, the entire calibration system is initialized and configured, specifically including: Before the drone takes off, the calibration system is initialized, including the following configurations: Sensor parameters: Input the intrinsic and extrinsic parameters and spectral response function of the multispectral and hyperspectral camera; input the sensitivity coefficient of the irradiance sensor. and dark current ; Flight parameters: Set the flight altitude according to the operating area and resolution requirements. Forward overlap 80%, Lateral overlap 70%, Flight speed and image acquisition frequency ; Algorithm parameters: Set the spectral stability threshold for PIF recognition. Spatial consistency threshold; Number of iterations for robust regression With convergence threshold and atmospheric correction parameters; Step S2: Synchronous Data Acquisition After the drone takes off according to the preset route, the synchronization control unit generates a cycle of... The synchronous trigger signals are sent to the main imaging sensor group and the irradiance measurement unit respectively, controlling both to synchronously start data acquisition: The multispectral camera and hyperspectral camera in the main imaging sensor group simultaneously acquire images of ground features, obtaining multispectral image data. and hyperspectral image data ,in For pixel coordinates, The characteristic wavelength band of the multispectral camera, For the continuous band wavelengths of the hyperspectral camera, For the number of bands in a multispectral camera, This represents the number of bands in the hyperspectral camera. The irradiance measurement unit simultaneously measures the downlink total solar radiation flux density to acquire irradiance data. ,in To collect timestamps, ensuring that each frame of image data corresponds to a unique irradiance data point, After the data acquisition module adds metadata information such as timestamps and sensor numbers to the acquired image data and irradiance data, it transmits the data to the data processing module in real time through the data transmission unit. Step S3: Automatic PIF Feature Identification and Algorithm Design The PIF automatic identification module in the data processing module processes the synchronously acquired multispectral and hyperspectral images and automatically extracts spectrally stable PIF ground features; Step S4: Cross-sensor radiation normalization The core of cross-sensor radiometric normalization is to establish the radiometric response conversion relationship between multispectral and hyperspectral sensors, normalize the output data of different sensors to a unified radiometric reference, and eliminate the radiometric bias caused by the differences in the response of the sensors themselves. This step takes the DN value of the PIF ground cover area as a reference and is achieved through three sub-steps: band matching, outlier removal, and robust regression modeling. Step S5: Absolute reflectivity conversion - Real-time conversion algorithm based on coupled online irradiance Even after cross-sensor normalization, the radiation data is still a relative value. Therefore, a "real-time absolute reflectance conversion algorithm coupled with online irradiance" is required. This algorithm uses the irradiance data collected in real time by the UAV's onboard sensors as a dynamic benchmark, and combines a simplified atmospheric correction model to eliminate environmental interference from light intensity fluctuations and atmospheric transmission attenuation. This process converts the normalized DN value into a standardized surface absolute reflectance.

3. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 2, characterized in that, In step S3, the specific details are as follows: Step S3.1: Image preprocessing S3.1.1 Radiation pretreatment: Eliminating system noise (1) Dark current correction is performed using a "dark field calibration + band-by-band subtraction" strategy. The mathematical model is as follows: , in: Pixel positions in the original image In the band The numerical value; The reference value for dark current in this band was obtained by capturing 50 dark field images with the lens blocked and taking the average value. (2) Dead pixel repair Repair is performed using a 3×3 neighborhood median filtering method: for each pixel, if its DN value deviates from the neighborhood mean by more than 3 times the standard deviation, it is judged as a bad pixel; Replace the pixel value with the median value of the neighborhood center. S3.1.2 Geometric preprocessing: Eliminating spatial distortion Orthorectification is performed using a rational polynomial coefficient RPC model to convert image pixel coordinates into ground geographic coordinates, achieving a consistent representation of ground feature locations. Corrected pixels Corresponding ground coordinates The calculation formula is as follows: in: : Pixel coordinates of the image before correction : Corresponding ground projection coordinates : RPC model coefficients, obtained from camera calibration experiments.

4. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 3, characterized in that, In step S3, the specific details are as follows: Step S3.2: Adaptive multi-scale image segmentation - generating a homogeneous region candidate set This paper proposes a region-level recognition strategy based on "spectrally similar homogeneous regions". A spectral-spatial collaborative adaptive mean-shift segmentation algorithm is designed. By pre-segmenting the image into homogeneous regions with consistent spectral features, subsequent recognition is performed on a region-by-region basis, improving the system's robustness and recognition accuracy in complex environments. S3.2.1 Data Input Input data: spectral image Image height, Image width, Number of bands; Preset parameters: empirical coefficient By band number Sure: Take 0.8 at that time. Time complexity is 1.2, spectral complexity coefficient Space kernel bandwidth coefficient Region constraint threshold, number of pixels Space compactness Spectral homogeneity ; S3.2.2 Algorithm Flow Step 1 Feature vector construction For each pixel in the spectral image Construct a joint spectral-spatial feature vector: Pixels of 3D spectral vector; Pixels The spatial coordinate vector; Step 2: Calculation of local spectral complexity For each pixel Take it Neighborhood Calculate the local spectral complexity: in For the neighboring region The variance of a band reflects the heterogeneity of the spectrum surrounding that pixel; Step 3: Adaptive calculation of kernel bandwidth Decompose the kernel bandwidth into spectral kernel bandwidth. With space kernel bandwidth Calculate separately: The value decreases as the local spectral complexity increases, adapting to the spectral characteristics of ground features in the spectral image. It is positively correlated with the area of ​​the neighborhood, ensuring the spatial neighborhood continuity of the pixel; Step 4: Iteration of mean drift in spectral-spatial dual-core weighted system For each pixel , by Starting from the initial point, perform mean-shift iteration: ① Determine the local region: Filter to meet the requirements and pixels , forming a local region ; ② Calculate the dual-core weighted mean: Where: weighting coefficient ; spectral nucleus Space Core ; ③ Iterative update: Update the initial point to Repeat steps 1-2 until... ; ④ Clustering and merging: Divide pixels that converge to the same mean into the same region. ; Step 5: Region Constraint Filtering For the segmented regions Regions that meet the following criteria are selected as the final homogeneous region candidate set: ②Pixel count constraint: ② Space compactness constraints: ( The number of pixels contained in the region. (Number of pixels at the edge of the region) ③ Spectral homogeneity constraint: (Std The standard deviation of the spectral vectors within the region, max (This is the global maximum value of the spectral vector of the spectral image); S3.2.3 Algorithm Output Output: Candidate set of homogeneous regions for this spectral image K represents the total number of regions, and each region satisfies the spectral-spatial homogeneity constraint.

5. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 4, characterized in that, In step S3, the specific details are as follows: Step S3.3: Spectral stability analysis - dual-spectral index verification For each homogeneous region The spectral feature vectors of the identified PIF ground features in multispectral and hyperspectral images are extracted, and their spectral stability is evaluated using a two-dimensional index of "geometric similarity + statistical dissimilarity" to ensure that the spectral features of the identified PIF ground features remain stable under different time periods and different sensors. S3.3.1 Calculation of Spectral Angle Matching Degree (SAM) For the region At adjacent data collection times average spectral vector and The SAM calculation formula is as follows: in: Spectral angle : region At any moment The average spectral reflectance or DN value, For wavelength, The set of common bands for multispectral and hyperspectral imaging, and the range of common bands. Spectral angle thresholds were set based on experimental experience data. , when At that time, it is assumed that the spectral characteristics of the region at the two moments are geometrically similar; S3.3.2 Calculation of Spectral Information Divergence (SID) set up , They are respectively regions The normalized spectral vectors at two time points, SID, are defined as follows: in The Kullback-Leibler divergence (KL divergence) is calculated using the following formula: in: Extremely small positive numbers to prevent division by zero errors. It needs to be pre-normalized to a probability distribution. Set SID experience threshold , When SID At that time, it was assumed that the regional spectral distribution had statistical stability. Regions that simultaneously meet the SAM and SID threshold conditions are marked as "spectrally stable candidate regions" and enter the subsequent spatial consistency verification stage.

6. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 5, characterized in that, In step S3, the specific details are as follows: Step S3.4: Spatial Consistency Verification - Removing False Stable Regions S3.4.1 Spatial Heterogeneity Analysis Computational area Spatial heterogeneity index This reflects the dispersion of pixel values ​​within a region, and the formula is as follows: in: :area The standard deviation of the DN values ​​of all pixels within the range, :area The mean of the DN values ​​of all pixels within the range. It is a dimensionless index used to measure spectral uniformity within a region. Specify spatial heterogeneity empirical threshold ,when When the pixel features within the region are uniform, the spatial consistency is good. S3.4.2 Neighborhood Continuity Analysis Statistical area The percentage of "spectrally stable candidate regions" in the 8-neighborhood The formula is as follows: in: :

8. The number of spectrally stable candidate regions in the neighborhood, with a value range of : , The percentage of stable regions within the neighborhood is normalized to the [0,1] interval, and an empirical threshold for neighborhood consistency is set. ,when This indicates that the region is spatially connected to surrounding stable regions, thus eliminating isolated, spurious stable regions. Regions that simultaneously meet the following two conditions are ultimately identified as PIF (Picture-in-Flight) feature regions, and the DN (Domain Number) values ​​of all their pixels are extracted as a reference sample set for subsequent cross-sensor radiometric normalization: 。 7. The automatic radiometric calibration method for multiple sensors on a UAV according to claim 2, characterized in that, In step S4, the specific details are as follows: Step S4.1: Band Matching Band matching is performed using the convolution integral method of spectral response functions, and the fusion formula is as follows: in: Multispectral bands after fusion The digital quantization value DN, Hyperspectral images at wavelength The DN value at that location, Multispectral camera Spectral response function SRM of the band, Multispectral bands The wavelength range boundary; Step S4.2: Outlier Removal S4.2.1 Initial Residual Calculation Multispectral DN values ​​within the PIF region Matched hyperspectral DN values The difference is used as the initial residual: in: : No. The residual of each pixel Multispectral cameras in band The pixel value, DN values ​​corresponding to hyperspectral data after band matching; S4.2.2 Weight Calculation The adaptive weight for each pixel is calculated based on the median absolute deviation (MAD) of the residuals, and the weighting function is as follows: in: The median absolute deviation (MAD) of the residuals. 1.4826: The conversion factor between MAD and standard deviation under normal distribution (making MAD ≈ σ) : No. The weight of each pixel is used to suppress the influence of outliers; S4.2.3 Iterative Optimization The model parameters are updated with weights, and new residuals are calculated based on the new parameters. This weight-parameter update process is repeated until the preset number of iterations is reached. Or the change in weights is less than the convergence threshold. ; S4.2.4 Outlier Deletion Remove weights The pixels, of which The weighted empirical threshold is used as the basis for determining the weights. The remaining pixels form the purified PIF sample set, which is used for subsequent radiometric normalization modeling.

8. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 7, characterized in that, In step S4, the specific details are as follows: Step S4.3: Multispectral-Hyperspectral Adaptive Weighted Robust Regressive Radiometric Normalization Algorithm S4.3.1 Algorithm Input Input data: ① Multispectral images in band DN value matrix ; ②DN value matrix of hyperspectral image after band matching ; ③ Cleaned PIF sample set ( (sample size) in, Multispectral cameras in spectral bands The DN value; : DN values ​​corresponding to hyperspectral data after band matching; Preset parameters: Weight iteration convergence threshold Maximum number of iterations ; S4.3.2 Algorithm Flow Step 1: PIF sample weight initialization, for each sample in the PIF sample set Initialize weights ; Step 2: Adaptive weighted robust regression modeling for multispectral bands Iterative solution of the linear transformation model parameters : ① Parameter solution for the t-th iteration: Construct the objective function with the goal of minimizing the weighted sum of squared residuals: right Taking the partial derivative and setting it to 0, we obtain the analytical solution: ， ② Weight update: Calculate the residual of the t-th iteration. Update weights based on robust residual estimation: in , ③ Convergence judgment: If or Stop iteration, take , Otherwise, repeat steps 1-2. Step 3: Cross-band radiative normalization is performed on all pixels in the hyperspectral image. The result obtained by solving , Convert hyperspectral DN values ​​to multispectral reference DN values: Step 4: Model Accuracy Verification Calculate the normalized residuals of the PIF sample set. Verify the following metrics: residual mean (Unbiasedness); residual standard deviation (Accuracy meets standards); S4.3.3 Algorithm Output Output result: ① Multispectral band Corresponding gain coefficient Offset coefficient ; ② Normalized hyperspectral DN value matrix; This enables radiometric consistency of multi-source remote sensing data, supporting subsequent quantitative analysis.

9. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 2, characterized in that, In step S5, the specific details are as follows: Step S5.1 Algorithm Preparatory Work S5.1.1 Input Data Structure Radiation fundamental data: Normalized DN value dataset across sensors ,in For pixel coordinates, For the spectral bands, this data has eliminated differences in multi-sensor responses; Dynamic irradiance data: Downlink total solar irradiance flux density synchronously acquired by the airborne irradiance sensor. , To collect timestamps, ensure a one-to-one correspondence with each frame of DN data; Auxiliary parameter data: UAV GPS / IMU data, sensor laboratory calibration parameters, and real-time meteorological parameters; S5.1.2 Conditional Assumptions Based on the theory of remote sensing radiative transfer, the following assumptions are made: Assumption 1: The sensor output DN value is linearly related to the input radiance, that is... This assumption forms the physical basis for radiance conversion; Assumption 2: When the UAV flies at low altitude, the aerosol is vertically uniformly distributed and has a low concentration, and its scattering effect can be approximated by a simplified model; Assumption 3: Downward solar radiation can be considered as uniformly distributed at the moment of observation. This assumption is guaranteed through sensor installation location optimization and data verification.

10. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 9, characterized in that, In step S5, the specific details are as follows: Step S5.2 Dynamic correction of irradiance data To obtain high-precision downlink solar irradiance flux density, multi-stage dynamic correction of the raw data is necessary, which is implemented as follows: S5.2.1 Dark Current Dynamic Correction The sensor still exhibits a weak output, known as dark current, even in the absence of light. This dark current needs to be subtracted in real time to eliminate baseline drift. Correction strategy: The instantaneous dark current is calculated using a "sliding window mean" strategy. The mathematical model is as follows: in: : Moment The corrected dark current value; The sliding window length is set to 5 frames to balance response speed and stability. : No. Dark field data of the frame; S5.2.2 Temperature Drift Compensation The sensitivity of the sensor drifts with changes in operating temperature and needs to be compensated for using a temperature coefficient. Compensation model: Linear temperature response model, temperature correction coefficient Described using a linear model: in: Standard temperature Sensitivity coefficient at the following levels : band Temperature coefficient (unit: The temperature was calibrated by high and low temperature environmental testing (-10℃ to 50℃). : Moment The measured operating temperature of the sensor (acquired by the built-in temperature sensor); S5.2.3 Random Noise Filtering A 5-point moving average filter was used to remove Gaussian random noise. The filtered irradiance data The calculation is as follows: ， S5.2.4 Final Irradiance Determination Based on the above correction steps, the time is obtained. band True downlink solar irradiance The calculation formula is as follows: in: : The original irradiance data after filtering; Dark current after sliding window correction; Temperature compensation coefficient.

11. The UAV-borne multi-sensor automatic radiometric calibration method according to claim 10, characterized in that, In step S5, the specific details are as follows: Step S5.3: Calculation of Top Atmosphere Radiance Based on the normalized DN value and combined with the camera's laboratory calibration parameters, the radiance of ground features at the top of the atmosphere is calculated. The preliminary relationship of the sensor response characteristics is established, and the calculation formula is as follows: in: Pixels At wavelength Atmospheric top radiance at that location. DN is the digitally quantized value after cross-sensor normalization, i.e., the data after radiometric normalization. Camera at wavelength The gain coefficient at a given point reflects the sensor's sensitivity. Offset coefficient, reflecting the zero-point deviation of the system; Sub-step S5.4: Atmospheric correction and absolute reflectance calculation Based on simplification The model undergoes atmospheric correction. Considering the characteristics of low-altitude flight of UAVs, the influence of complex aerosol scattering is ignored, and only atmospheric molecular scattering and water vapor absorption are considered to calculate the absolute surface reflectance. The formula is as follows: in: : Absolute surface reflectance Atmospheric radiance Atmospheric path radiation Earth-Sun distance correction factor Corrected solar irradiance : Solar zenith angle Atmospheric transmittance S5.4.1 Atmospheric path radiation Through dark-field experiments or by selecting deep water bodies and shadowy dark areas in images... Estimate using the following formula: ， S5.4.2 Earth-Sun Distance Correction Factor Since Earth's orbit around the sun is elliptical, it is necessary to correct for the effect of changes in the Earth-Sun distance on solar radiation. The calculation formula is as follows: in: DOY: Yearly Accumulation Days, Value Range Or 366, The relative proportion of the Earth-Sun distance to the average distance. S5.4.3 Solar Zenith Angle Calculated based on the real-time geographic location and data collection time of the drone: First, calculate the solar altitude angle. : Then the solar zenith angle is: in: Latitude of the observation point The solar declination angle is calculated using the following formula: ， ω: Hour angle, calculated using the following formula: ， Collection time, S5.4.4 Atmospheric transmittance Divided into atmospheric molecular transmittance and water vapor transmission rate : , in: Atmospheric optical thickness, obtained from actual measurements at local meteorological stations or retrieved through hyperspectral data inversion. Molecular scattering term, calculated using the Rayleigh scattering model: ， Rayleigh scattering coefficient, relative to wavelength Related, Atmospheric water vapor content (unit: g / cm²) can be obtained through water vapor absorption spectrum inversion or by querying the US Standard Atmosphere database. The water vapor absorption term is obtained from empirical models or by looking up tables.

Images