A method and system for karst landform identification based on UAV hyperspectral imaging

Abstract

This invention discloses a method and system for karst landform identification based on UAV hyperspectral imaging, relating to the field of UAV remote sensing technology. The method includes acquiring range data of the area to be explored and conducting surveying and deployment; determining the UAV take-off and landing points and flight path map based on the surveying and deployment results; collecting hyperspectral data, ground-measured benchmark data, and topographic spatial reference data according to the flight path map based on the UAV take-off and landing points; performing radiometric geometric joint correction on the hyperspectral data based on the ground-measured benchmark data and topographic spatial reference data to obtain corrected hyperspectral reflectance data; and using the corrected hyperspectral reflectance data, topographic spatial reference data, and ground-measured benchmark data to identify karst landforms, thus obtaining the karst landform identification results for the area to be explored. This method achieves low-cost, high-efficiency, and full-area coverage of karst landform exploration through UAVs, and ensures high accuracy in identifying karst geological features through multi-source data collaborative correction and fusion identification.

Description

Technical Field

[0001] This invention relates to the field of UAV remote sensing technology, and in particular to a method and system for identifying karst landforms based on UAV hyperspectral imaging. Background Technology

[0002] In recent years, environmental emergencies related to petroleum pollutants have occurred frequently. When these products leak due to production accidents, traffic accidents, or illegal discharges, they can easily seep and spread rapidly along underground channels and surface in areas with complex geological conditions and well-developed fissures, causing cross-regional pollution and posing a serious threat to the ecological environment and drinking water safety. Responding to such emergencies requires, first and foremost, accurately identifying the basic structure of the regional geological system to lay a theoretical foundation for effectively responding to pollution incidents.

[0003] Response speed is crucial for handling such emergencies. Existing geological structure identification methods rely on a single detection data source, which is not effective in identifying geological structures in complex geological areas. In emergency response, accurate identification mainly relies on manual on-site detection. However, due to the complexity of geological structures and limited operating conditions, manual detection is costly, inefficient, and difficult, making it hard to meet the requirements of emergency response.

[0004] Single-source detection methods are limited by their own characteristics and often fail to capture the complex and ever-changing geological structural features comprehensively and accurately. Existing data processing and analysis technologies generally suffer from data redundancy and low information extraction efficiency when processing multiple types of detection data, which restricts the speed and accuracy of detection results. Therefore, there is an urgent need for a method to detect and identify karst landforms in order to quickly and accurately ascertain the basic situation of the geological structure and provide data support for effectively responding to sudden environmental events. Summary of the Invention

[0005] This invention provides a method and system for karst landform identification based on UAV hyperspectral imaging, which solves the technical problem that existing karst landform detection methods cannot simultaneously achieve low-cost, high-efficiency data acquisition and high-precision geological feature identification.

[0006] The first aspect of this invention provides a method for identifying karst landforms based on UAV hyperspectral imaging, comprising:

[0007] Acquire the range data of the area to be explored and conduct survey and deployment. Based on the survey and deployment results, determine the take-off and landing points and flight route map of the UAV.

[0008] Based on the take-off and landing point of the UAV, hyperspectral data and ground-measured reference data of the area to be detected are collected according to the flight route map, and terrain spatial reference data are acquired simultaneously.

[0009] Based on the ground-measured reference data and the topographic spatial reference data, the hyperspectral data is subjected to radiometric geometric joint correction to obtain corrected hyperspectral reflectance data.

[0010] Karst landform identification is performed using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground measured benchmark data to obtain the karst landform identification results for the area to be detected.

[0011] Optionally, the step of acquiring the range data of the area to be detected and conducting surveying and deployment, and determining the UAV take-off and landing points and flight route map based on the surveying and deployment results, includes:

[0012] Obtain the range data of the area to be detected;

[0013] Based on the range data, obtain the edge survey information of the area to be detected;

[0014] Based on the on-site edge survey information, the location layout information of ground measurement base stations and image control points is determined within the area to be detected;

[0015] The take-off and landing points of the UAV are determined by using the on-site edge survey information and the point layout information as constraints.

[0016] A flight path map is generated based on the range data and the on-site survey information of the edge.

[0017] Optionally, the step of performing radiometric geometric joint correction on the hyperspectral data based on the ground-measured reference data and the topographic spatial reference data to obtain corrected hyperspectral reflectance data includes:

[0018] The hyperspectral data were radiometrically calibrated to obtain hyperspectral radiance data;

[0019] Based on the ground-measured benchmark data, the hyperspectral radiance data is converted to spectral reflectance to obtain hyperspectral reflectance data.

[0020] Using the ground-measured benchmark data and the topographic spatial reference data, the hyperspectral reflectance data is geometrically corrected to obtain corrected hyperspectral reflectance data.

[0021] Optionally, the step of radiometrically calibrating the hyperspectral data to obtain hyperspectral radiance data includes:

[0022] Extract the hyperspectral raw data cube from the hyperspectral data;

[0023] Extract the original DN value matrix corresponding to each spectral band from the hyperspectral raw data cube;

[0024] Based on the original DN value matrix and the radiation intensity of the preset standard plate, calculate the radiation calibration coefficient corresponding to each spectral band;

[0025] A radiometric calibration transformation function is constructed based on the radiometric calibration coefficients, and the original DN value matrices are substituted into the function to obtain the two-dimensional radiance matrix corresponding to each spectral band.

[0026] The two-dimensional radiance matrices are reorganized according to the spectral band order to obtain the hyperspectral radiance data of the region to be detected.

[0027] Optionally, the ground-measured reference data includes ground object spectral data, and the step of performing spectral reflectance conversion on the hyperspectral radiance data based on the ground-measured reference data to obtain hyperspectral reflectance data includes:

[0028] Extract measured data of standard targets from the spectral data of the ground features;

[0029] Using the measured data of the standard target and the preset standard target reference reflectance, the illumination correction coefficient corresponding to each of the spectral bands is calculated.

[0030] The target reflectivity conversion function is obtained by correcting the preset basic reflectivity conversion model using the illumination correction coefficient.

[0031] Substitute each two-dimensional radiance matrix in the hyperspectral radiance data into the target reflectance conversion function to obtain the two-dimensional reflectance matrix corresponding to each spectral band.

[0032] The two-dimensional reflectance matrices are reorganized according to the spectral band order to obtain the hyperspectral reflectance data of the region to be detected.

[0033] Optionally, the terrain spatial reference data includes DEM data and orthophoto data, the ground-measured reference data includes image control point data, and the step of using the ground-measured reference data and the terrain spatial reference data to perform geometric correction on the hyperspectral reflectance data to obtain corrected hyperspectral reflectance data includes:

[0034] Using the orthophoto data as a spatial reference, the initial pixel coordinates corresponding to each two-dimensional reflectance matrix within the hyperspectral reflectance data are determined;

[0035] Calculate the actual ground elevation corresponding to each pixel in the area to be detected based on the DEM data;

[0036] Based on the reference plane elevation corresponding to the orthophoto data and the actual ground elevation, the positional offset of each pixel is calculated as an elevation correction term.

[0037] Using the measured planar coordinates of the control point data as the truth benchmark, and combining the initial pixel coordinates with the elevation correction term, a geometric correction mathematical model is constructed.

[0038] The model coefficients of the geometric correction mathematical model are solved using the least squares method to obtain the target geometric correction model;

[0039] The hyperspectral reflectance data is substituted into the target geometric correction model for spatial transformation to obtain the corrected hyperspectral reflectance data.

[0040] Optionally, the step of using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground-measured benchmark data to identify karst landforms and obtain karst landform identification results for the area to be detected includes:

[0041] The corrected hyperspectral reflectance data is preprocessed, including interference masking, dimensionality reduction, and noise suppression.

[0042] Extract the measured limestone spectrum from the aforementioned ground feature spectral data;

[0043] Based on the spectral angle mapping method, using the measured limestone spectrum as the end element and setting an angle threshold, the angle between the spectral vector of each pixel in the preprocessed corrected hyperspectral reflectance data and the limestone reference spectral vector is calculated.

[0044] Pixels with included angles below the specified angle threshold are marked as carbonate rock target pixels to obtain carbonate rock extraction results.

[0045] Linear spectral unmixing was performed on the carbonate rock extraction results to determine the carbonate rock region;

[0046] Calculate the terrain factor corresponding to each pixel in the area to be detected based on the DEM data;

[0047] Based on the aforementioned topographic factors, topographic features are screened to determine candidate regions for peak forests and depressions, respectively.

[0048] The candidate areas of peak forests and depressions are spatially superimposed with the carbonate rock areas to obtain the karst landform identification results of the area to be detected.

[0049] Optionally, it also includes:

[0050] Morphological screening was performed on the candidate region set of peak forests and depressions to obtain a second-screened candidate region set of peak forests and depressions.

[0051] The candidate areas of peak forests and depressions selected in the secondary screening are spatially superimposed with the carbonate rock area to obtain new karst landform identification results for the area to be explored.

[0052] Optionally, it also includes:

[0053] Based on the corrected hyperspectral reflectance data, the relative surface water content data of the area to be detected is retrieved.

[0054] Extract the relative surface water content values ​​of each candidate region within the secondary screening set of peak forest and depression candidate regions from the relative surface water content data.

[0055] Determine whether the relative surface water content values ​​corresponding to each candidate region meet the associated preset hydrological determination conditions;

[0056] Candidate areas that meet the preset hydrological judgment conditions are retained to obtain a set of candidate areas for peak forests and depressions after three rounds of screening;

[0057] The candidate areas of peak forests and depressions selected through the three screenings are spatially superimposed with the carbonate rock area to obtain new karst landform identification results for the area to be explored.

[0058] A second aspect of the present invention provides a karst landform identification system based on UAV hyperspectral imaging, comprising:

[0059] The acquisition module is used to acquire range data of the area to be detected and to conduct survey and deployment. Based on the survey and deployment results, the take-off and landing points and flight route map of the UAV are determined.

[0060] The data acquisition module is used to acquire hyperspectral data and ground-measured reference data of the area to be detected based on the take-off and landing point of the UAV and according to the flight route map, and simultaneously acquire terrain spatial reference data.

[0061] The correction module is used to perform radiometric geometric joint correction on the hyperspectral data based on the ground measured reference data and the topographic spatial reference data to obtain corrected hyperspectral reflectance data.

[0062] The identification module is used to identify karst landforms using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground measured benchmark data, and to obtain the karst landform identification results of the area to be detected.

[0063] As can be seen from the above technical solutions, the present invention has the following advantages:

[0064] This invention provides a method and system for karst landform identification based on UAV hyperspectral imaging. Relying on precise planning of UAV take-off and landing points and flight routes through prior surveying and deployment, the UAV simultaneously acquires hyperspectral data, ground-measured benchmark data, and topographic spatial reference data. This enables the acquisition of comprehensive karst area data in a low-cost and high-efficiency manner, avoiding the limitations of high cost, low efficiency, high operational risks, and limited coverage of manual on-site karst landform exploration. This ensures the economy and timeliness of the karst landform exploration process. Furthermore, the hyperspectral data undergoes radiometric and geometric joint correction using ground-measured benchmark data and topographic spatial reference data, simultaneously achieving spectral information fidelity and spatial location accuracy optimization. This effectively compensates for the disconnect between spectral and geometric information processing in traditional correction methods, avoiding the information loss caused by single-dimensional correction. Information deviations are eliminated to ensure the spectral reliability and spatial accuracy of the corrected data, providing a high-quality data foundation for karst landform lithology identification and spatial morphology characterization. Subsequent fusion and correction of hyperspectral reflectance data, topographic spatial reference data, and ground-measured benchmark data for karst landform identification can accurately capture lithological differences, fissure development characteristics, and landform details in karst areas, effectively improving the accuracy and reliability of karst geological feature identification. This invention achieves low-cost, high-efficiency, and comprehensive coverage of karst landform exploration through UAVs, and ensures high accuracy in karst geological feature identification through multi-source data collaborative correction and fusion identification. It effectively solves the problem in existing karst landform exploration technologies that cannot simultaneously achieve low-cost, high-efficiency exploration and high-precision geological feature identification, providing reliable technical support for rapid and accurate exploration of karst areas. Attached Figure Description

[0065] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0066] Figure 1 This is a flowchart illustrating the steps of a karst landform identification method based on UAV hyperspectral imaging, provided in Embodiment 1 of the present invention.

[0067] Figure 2 This is a flowchart illustrating the steps of a karst landform identification method based on UAV hyperspectral imaging, provided in Embodiment 2 of the present invention.

[0068] Figure 3 This is a schematic diagram of illumination changes provided in Embodiment 2 of the present invention;

[0069] Figure 4This is a structural block diagram of a karst landform identification system based on UAV hyperspectral imaging provided in Embodiment 3 of the present invention. Detailed Implementation

[0070] This invention provides a method and system for identifying karst landforms based on UAV hyperspectral imaging, which addresses the technical problem that existing karst landform detection methods struggle to balance low-cost, high-efficiency data acquisition with high-precision geological feature identification.

[0071] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0072] Please see Figure 1 , Figure 1 This is a flowchart illustrating the steps of a karst landform identification method based on UAV hyperspectral imaging, as provided in Embodiment 1 of the present invention.

[0073] This invention provides a method for identifying karst landforms based on UAV hyperspectral imaging, comprising:

[0074] Step 101: Obtain the range data of the area to be detected and conduct survey and layout. Based on the survey and layout results, determine the take-off and landing points and flight route map of the UAV.

[0075] In this embodiment of the invention, the first step is to acquire the range data of the karst area to be explored. This data can be retrieved through a basic geographic information platform or directly input by the user as the region's boundary vector information. It also includes basic information such as the region's approximate elevation and the distribution of major features. Based on the acquired range data, a survey and deployment are conducted, focusing on investigating the safe take-off and landing conditions of the area to be explored and its surroundings. Open areas free from tall obstacles and high-voltage line interference are marked as candidate locations for UAV take-off and landing points. Simultaneously, key areas affecting flight safety, such as those with mountain obstructions and complex airflow, are identified. Combining the survey and deployment results, with the goal of covering the entire area to be explored without blind spots, a main flight path is planned based on the shape of the region's boundaries. This avoids identified safety risk areas, and transitional routes are set between the take-off and landing points and the main flight path. Finally, a flight route map is generated, including take-off and landing point coordinates, route node coordinates, flight altitude, and safe operating boundaries.

[0076] Step 102: Based on the UAV take-off and landing point, collect hyperspectral data and ground-measured benchmark data of the area to be detected according to the flight route map, and simultaneously acquire terrain spatial reference data.

[0077] In this embodiment of the invention, the UAV takes off from a determined take-off and landing point and carries out flight operations according to the flight route, flight altitude, and safe operating boundaries planned in the flight route map. During the operation, the hyperspectral imager carried by the UAV collects hyperspectral image data of the area to be detected according to preset operating parameters, and simultaneously records position and attitude auxiliary data. Ground personnel at pre-deployed control points simultaneously collect the three-dimensional coordinates of the control points and the spectral data of ground objects as ground measurement reference data. At the same time, the UAV simultaneously collects lidar point cloud data of the area to be detected during flight, or retrieves pre-stored regional digital elevation model data as terrain spatial reference data, completing the synchronous collection and storage of multiple types of detection data.

[0078] Step 103: Based on ground-measured benchmark data and topographic spatial reference data, perform radiometric geometric joint correction on the hyperspectral data to obtain corrected hyperspectral reflectance data.

[0079] In this embodiment of the invention, the hyperspectral data is first radiometrically calibrated using standard plate radiance data from ground-measured reference data as the calibration basis, and then uniformly converted into radiance data corresponding to each band. Next, the radiance data is converted to spectral reflectance based on ground-measured spectral reflectance data to obtain hyperspectral reflectance data. Simultaneously, using elevation and planar coordinate information from topographic spatial reference data, combined with control point coordinates from ground-measured reference data, a geometric correction model is constructed to perform spatial position correction on the hyperspectral reflectance data, ultimately obtaining corrected hyperspectral reflectance data that simultaneously possesses both faithful spectral information and accurate spatial position information.

[0080] Step 104: Karst landform identification is performed using calibrated hyperspectral reflectance data, topographic spatial reference data, and ground-measured benchmark data to obtain karst landform identification results for the area to be explored.

[0081] In this embodiment of the invention, firstly, based on corrected hyperspectral reflectance data and combined with karst landform spectral samples from ground-measured benchmark data, key spectral features such as lithology and karst fissures are extracted to complete the initial spectral classification of karst-related landforms. Then, topographic feature information such as elevation, slope, and curvature from topographic spatial reference data is superimposed to supplement the morphological features of karst landforms. Subsequently, spectral and topographic feature data are fused, and comprehensive classification and identification are carried out through a preset karst landform identification model to initially obtain the karst landform distribution results of the area to be explored. Finally, using verification control points from ground-measured benchmark data, the identification results are checked and corrected for accuracy, and the final karst landform identification result is output. This result includes the distribution range and boundary information of different karst landform types within the area to be explored.

[0082] UAV hyperspectral imaging refers to the technology of acquiring continuous spectral and spatial information of ground objects by using UAVs equipped with hyperspectral imaging equipment. Karst landforms are landform types such as peak forests and depressions formed by the dissolution and erosion of carbonate rocks by water. The area to be detected is the geographical range for which karst landform identification needs to be carried out. The range data is the geographical range information that characterizes the boundary and area of ​​the area to be detected. Survey and layout is the planning and layout of stations and points for UAV data collection and ground measurement. The UAV take-off and landing point is the designated site for UAV take-off and landing. The flight route map is the flight path planning map that guides the UAV to collect data according to the preset trajectory. Hyperspectral data is the raw data containing the spectral and spatial information of ground objects collected by the UAV hyperspectral equipment. Ground measured reference data is the reference data obtained by field measurement for calibration and auxiliary identification. Topographic spatial reference data is the data that provides geographic spatial reference and topographic features. Radiometric and geometric joint correction is a comprehensive processing process that simultaneously corrects the radiometric accuracy and geometric accuracy of hyperspectral data. Corrected hyperspectral reflectance data is the data that accurately characterizes the spectral reflectance characteristics of ground objects after eliminating radiometric and geometric distortions. The karst landform identification result is the final identification information such as the distribution and type of karst landforms in the area to be detected.

[0083] This embodiment establishes a complete and collaborative technical path for karst landform exploration through an integrated process of survey and deployment planning, simultaneous acquisition of multi-source data, joint radiation geometric correction, and multi-source fusion identification. Early survey and deployment, along with flight path planning, proactively identify and mitigate flight safety risks associated with the complex terrain of karst areas, ensuring the continuity and full coverage of subsequent data acquisition. Furthermore, by relying on UAVs to simultaneously acquire multiple types of exploration data, the operational costs and safety risks of manual field exploration in karst areas are significantly reduced, achieving efficient acquisition of full-area exploration data and avoiding the high cost and low efficiency of traditional manual exploration. Limitations; however, by simultaneously performing radiation-geometric joint correction to ensure spectral fidelity and optimize spatial location accuracy of hyperspectral data, this invention effectively overcomes the shortcomings of traditional step-by-step correction methods where spectral and geometric information processing are disconnected, providing high-quality foundational data for karst lithology identification and spatial morphology characterization. Finally, by fusing corrected hyperspectral reflectance data, topographic spatial reference data, and ground-measured benchmark data, karst landform identification is conducted, fully leveraging the complementary advantages of spectral features, topographic features, and ground benchmark data to accurately capture lithological differences and landform details in karst areas, effectively improving the accuracy and reliability of karst landform identification results. This invention achieves low-cost and high-efficiency karst landform detection through UAVs, while ensuring the accuracy of karst geological feature identification through high-precision data processing and multi-source fusion identification. It effectively solves the core problem in existing karst landform detection technologies that cannot simultaneously achieve low-cost, high-efficiency detection and high-precision geological feature identification, providing practical technical support for rapid and accurate detection in karst areas.

[0084] Please see Figure 2 , Figure 2 This is a flowchart illustrating the steps of a karst landform identification method based on UAV hyperspectral imaging, as provided in Embodiment 2 of the present invention.

[0085] This invention provides a method for identifying karst landforms based on UAV hyperspectral imaging, comprising:

[0086] Step 201: Obtain the range data of the area to be explored and conduct survey and layout. Based on the survey and layout results, determine the take-off and landing points and flight route map of the UAV.

[0087] Further, step 201 may include the following sub-steps:

[0088] S11. Obtain the range data of the area to be detected.

[0089] In this embodiment of the invention, the range of the karst area to be explored is first determined, and the corresponding geographical range data is obtained, including basic data such as the geographical coordinates, boundary coordinates, and regional description information of the area to be explored. At the same time, the types and characteristics of land features in the exploration area are initially identified based on the basic data, such as typical karst landform types such as water bodies, soil, rock bodies, peak forests, and depressions.

[0090] S12. Based on the range data, obtain the edge survey information of the area to be detected.

[0091] In this embodiment of the invention, based on the acquired data of the area to be detected, an on-site survey is conducted on the edge of the survey area to obtain information on the topography, landforms, distribution of land features, traffic conditions and surrounding environment of the edge area. Key information such as the terrain undulations, obstacle distribution and electromagnetic interference source location of the edge area is investigated to form a complete edge on-site survey report.

[0092] S13. Based on the on-site survey information of the edge, determine the location layout information of the ground measurement base station and the image control point within the area to be detected.

[0093] In this embodiment of the invention, based on the edge field survey information, the layout of ground measurement base stations and image control points is determined within the area to be detected. The ground measurement base stations need to be deployed in locations with open fields of view and unobstructed signals within the survey area for real-time calibration of GPS data during UAV flight. The image control points need to be evenly distributed within the survey area, with priority given to locations with obvious ground features and easy identification, for subsequent geometric correction and accuracy verification of hyperspectral images, ensuring that the distribution of points can cover the entire area to be detected.

[0094] S14. Determine the take-off and landing points of the UAV based on the edge field survey information and point layout information.

[0095] In this embodiment of the invention, by combining edge field survey information and point layout information, a UAV take-off and landing point is selected in a flat and open location at the edge of the area to be detected. The take-off and landing point must avoid no-fly zones and restricted areas to ensure no conflict with other aircraft. At the same time, it should be far away from obstacles such as tall buildings and high-voltage power lines, and avoid strong electromagnetic interference sources such as radar stations and base stations to ensure the safety and stability of UAV take-off, landing and flight. Priority is given to locations with good communication conditions with ground measurement base stations to facilitate signal transmission and data interaction during flight.

[0096] S15. Generate a flight path map based on the range data and edge field survey information.

[0097] In this embodiment of the invention, based on the range data of the area to be detected and the on-site survey information of the edge, the ground resolution is first determined according to the detection accuracy requirements of karst landforms, and then the key flight parameters of the UAV, such as flight altitude, flight speed, heading overlap rate, and lateral overlap rate, are determined. Then, combined with the boundary range of the area to be detected and the minimum turning radius of the UAV, the main flight route is planned to ensure that the route can cover the entire area to be detected without blind spots. At the same time, a transition route is planned between the take-off and landing point and the main route to avoid safety risk areas such as mountain obstruction and complex airflow identified during the survey. Finally, a complete flight route map containing the coordinates of the take-off and landing point, the coordinates of the route node, the flight altitude, the flight speed, and the safe operation boundary is generated to ensure the continuity and full coverage of subsequent data collection operations.

[0098] Edge field survey information refers to the terrain, obstacles, and other related information obtained from the edge survey of the area to be explored. Ground measurement base stations are ground stations deployed for precise positioning and data calibration. Image control points are control points with precise field coordinates used for geometric correction. Point layout information refers to the location and distribution density of ground measurement base stations and image control points.

[0099] Step 202: Based on the UAV take-off and landing point, collect hyperspectral data and ground-measured benchmark data of the area to be detected according to the flight route map, and simultaneously acquire terrain spatial reference data.

[0100] In this embodiment of the invention, a suitable period of operation with suitable light intensity, low wind speed and stable weather is first selected. The UAV takes off from the determined UAV take-off and landing point and carries out flight operations along the planned route according to the flight parameters preset in the flight route map (including flight altitude, heading overlap rate, lateral overlap rate, etc.). During the operation, the hyperspectral imaging equipment carried by the UAV collects image data of the area to be detected at a preset sampling frequency, and simultaneously records the position, attitude information, flight number and flight path number of the UAV, and generates hyperspectral data of the area to be detected.

[0101] While collecting hyperspectral data, ground personnel simultaneously carried out the collection of ground-based benchmark data: they went to the designated locations, used high-precision positioning equipment to collect the three-dimensional coordinate data of the locations, and took on-site images of the locations; for typical ground features in the area to be detected, they used ground feature spectral acquisition equipment to collect their ground spectral data, and recorded the location of each collection point, the collection time, and the corresponding UAV flight sorties and flight paths, forming complete ground-based benchmark data.

[0102] In addition, during flight, the drone can simultaneously collect relevant data that can be used to build a digital elevation model, or directly retrieve existing high-precision terrain data of the area to be explored as terrain spatial reference data. This enables the simultaneous collection and storage of hyperspectral data, ground-measured benchmark data, and terrain spatial reference data, providing a complete multi-source data foundation for subsequent data processing and karst landform identification.

[0103] Step 203: Radiometrically calibrate the hyperspectral data to obtain hyperspectral radiance data.

[0104] Furthermore, step 203 may include the following sub-steps:

[0105] S21. Extract the hyperspectral raw data cube from the hyperspectral data.

[0106] In this embodiment of the invention, after the hyperspectral data acquisition is completed, the stored raw hyperspectral image data is first parsed and the dimensions are extracted to construct a hyperspectral raw data cube containing spectral and spatial dimensions. The data cube uses the row and column numbers of spatial pixels as spatial indexes and the spectral band numbers as spectral indexes to completely record the raw imaging data of each spatial pixel in each spectral band within the area to be detected.

[0107] S22. Extract the original DN value matrix corresponding to each spectral band from the original hyperspectral data cube.

[0108] In this embodiment of the invention, based on the extracted hyperspectral raw data cube, the data is split band by band according to the spectral band number, and the two-dimensional pixel matrix corresponding to each spectral band is extracted, thus obtaining the raw DN value matrix corresponding to each spectral band. The row and column dimensions of each raw DN value matrix correspond to the spatial pixel rows and columns of the hyperspectral image, respectively. The element value in the matrix is ​​the image output DN value of the corresponding pixel in that band, which represents the quantized value of the raw radiation signal intensity received by the sensor. Since the DN value output by the imaging spectrometer itself does not have physical meaning, it needs to be converted into a quantifiable radiance value through subsequent radiometric calibration.

[0109] S23. Calculate the radiation calibration coefficients corresponding to each spectral band based on the original DN value matrix and the radiation intensity of the preset standard plate.

[0110] In this embodiment of the invention, a visible-shortwave infrared band reflected radiation calibration method is used, and standard plate radiance data synchronously acquired during hyperspectral imaging is selected as the calibration benchmark; the radiance of the standard plate is a spectral radiance value with known physical meaning, denoted as . The DN value corresponding to the standard plate acquired by the imaging spectrometer is denoted as The mapping relationship between spectral radiance and DN value is established based on a linear calibration model, and its expression is as follows:

[0111]

[0112] In the formula, The spectral radiance obtained within the standard field of view of the imaging spectrometer, in units of ; The DN value for the image output; The slope in the radiation calibration coefficient; This is the intercept in the radiation calibration coefficients. The radiation calibration coefficients for each spectral band can be calculated by using the ratio of the known radiation intensity on a standard plate to the instrument's measured DN value. and This enables the quantitative conversion of remote sensing information.

[0113] S24. Construct a radiometric calibration transformation function based on the radiometric calibration coefficients, and substitute it into each original DN value matrix to solve for the two-dimensional radiance matrix corresponding to each spectral band.

[0114] In this embodiment of the invention, the radiometric calibration of the imaging spectrometer is performed band by band. Based on the dynamic range of the imaging spectrometer, Equation 1 is extended to band-level calculation to obtain the relationship between the radiance input value and the remote sensor output DN value, expressed as:

[0115]

[0116] In the formula, For the first Group 1 Band radiance input value, in units of ; For the first Group 1 Band image grayscale output value; , For the first Group 1 The band radiation calibration coefficients correspond to the slope in Equation 1. With intercept The specific values ​​for each band.

[0117] Radiometric calibration coefficients calculated based on each spectral band , A radiometric calibration transformation function for the band is constructed. The original DN value matrix corresponding to the band is substituted into the transformation function, and the radiance value corresponding to each spatial pixel is calculated through linear transformation. Finally, a two-dimensional radiance matrix corresponding to each spectral band is generated.

[0118] S25. Reassemble each two-dimensional radiance matrix according to the spectral band order to obtain the hyperspectral radiance data of the region to be detected.

[0119] In this embodiment of the invention, the two-dimensional radiance matrices corresponding to all spectral bands are reorganized and stitched together according to the original spectral band order to restore the dimensional structure consistent with the original hyperspectral data cube, thereby constructing hyperspectral radiance data of the area to be detected. This hyperspectral radiance data is a spectral data cube with clear physical meaning, and the pixel value corresponding to each spectral band is a calibrated radiance value, which can be directly used for subsequent processing such as spectral reflectance inversion, radiometric correction, and karst landform feature identification, realizing the accurate conversion of hyperspectral data from original quantized values ​​to physical radiance.

[0120] Radiometric calibration is the process of converting raw hyperspectral DN values ​​into radiance data. Hyperspectral radiance data is data that characterizes the radiation intensity of ground objects in each band after radiometric calibration. Spectral reflectance conversion is the process of converting radiance data into reflectance data to eliminate the influence of illumination and atmosphere. Hyperspectral reflectance data is the core identification data that characterizes the spectral reflectance capability of ground objects in each band. Geometric correction is the process of correcting the spatial position distortion of hyperspectral data so that pixels correspond to accurate ground coordinates.

[0121] The hyperspectral raw data cube is a three-dimensional raw hyperspectral data matrix containing rows, columns, and spectral bands. The raw DN value matrix is ​​a two-dimensional matrix composed of the raw gray values ​​of single-spectral band pixels. The preset standard plate radiance intensity is the radiance intensity data of the standard reflector plate used as the radiometric calibration reference. The radiometric calibration coefficients are the coefficients that establish the conversion relationship between DN values ​​and radiance. The radiometric calibration conversion function is a mathematical relationship between DN values ​​and radiance constructed based on the calibration coefficients. The two-dimensional radiance matrix is ​​a two-dimensional matrix composed of the radiance values ​​of single-band pixels.

[0122] Step 204: Based on the ground-based measured benchmark data, perform spectral reflectance conversion on the hyperspectral radiance data to obtain hyperspectral reflectance data.

[0123] Furthermore, the ground-based measured reference data includes spectral data of ground features, and step 204 may include the following sub-steps:

[0124] S31. Extract measured data of standard targets from ground object spectral data.

[0125] In this embodiment of the invention, the ground object spectral data includes ground-measured spectral data of standard targets (such as standard target reflective cloth) deployed in the area to be detected, as well as target image data collected by the UAV during flight. First, the image grayscale values ​​(DN values) of the standard targets under different flight sorties, different flight paths, and different spectral bands are extracted from the ground object spectral data, including target image data collected by the UAV before and during its return from passing over the standard targets. At the same time, the reference reflectance data of the standard targets in each band monitored by the ground spectrometer and the illumination change monitoring data at different times are extracted.

[0126] S32. Using the measured data of the standard target and the preset standard target reference reflectance, calculate the illumination correction coefficient corresponding to each spectral band.

[0127] In this embodiment of the invention, reflectance calculation employs a real-time reference target method to eliminate the influence of atmospheric disturbances and changes in illumination radiation. First, the UAV flies over the standard target cloth for the first time before the survey line (e.g., at 10:00 AM) to collect the initial image grayscale value of the target cloth, denoted as... During subsequent flights, as the drone passes over the standard target in different sorties and flight paths, it will collect images of the target again, obtaining image grayscale values ​​at different times. Combined with a diagram illustrating changes in illumination, such as... Figure 3 As shown in the figure, this graph visually illustrates the changing trend of light intensity over time in different spectral bands (such as 1000nm and 1600nm). The horizontal axis represents time, the vertical axis corresponds to the image grayscale value, and the peak of the curve corresponds to the time of strongest light at noon (14:30). The baseline grayscale value is at 10:00 AM. The ratio of the two values ​​represents the grayscale value for the subsequent flight period, and is the illumination correction factor for that period. The calculation formula is:

[0128]

[0129] In the formula, For drones Gray values ​​of images captured each time the target cloth is passed through; For the first Lighting correction factor for different time periods; The initial image grayscale value (baseline grayscale value) is the value collected when the UAV first passes over the target cloth.

[0130] Based on the preset baseline reflectance of the standard target, the initial reflectance calculation relationship is constructed using the real-time reference method. The basic calculation formula is as follows:

[0131]

[0132] In the formula, For the first The reflectivity of the line image; For the first Grayscale values ​​of the line image; For the first The grayscale value of the target cloth corresponding to the line image (the grayscale value has been corrected by the illumination correction factor).

[0133] This method enables matching and converting reflectance changes every second, effectively addressing errors caused by atmospheric disturbances and variations in light radiation, and improving the accuracy of reflectance calculations.

[0134] S33. The preset basic reflectivity conversion model is corrected using the illumination correction coefficient to obtain the target reflectivity conversion function.

[0135] In this embodiment of the invention, the reflectance spectrum data of the standard target reflective cloth is first measured to obtain its reference reflectance in each spectral band. Based on this reference reflectance, combined with the illumination correction coefficient calculated in step S32, the basic reflectance model in Equation 3 is corrected to eliminate the systematic deviation caused by illumination changes and atmospheric disturbances at different times. During the correction process, the radiation spectrum value of the corresponding pixel point of the standard target on the aerial survey image is used as a reference to calibrate the model parameters, and finally a target reflectance conversion function applicable to each spectral band is obtained. This function integrates the standard target reference reflectance, real-time illumination correction and spectral band characteristics, and can adapt to the illumination conditions of different flight times in the area to be detected, ensuring the stability and accuracy of reflectance conversion.

[0136] S34. Substitute each two-dimensional radiance matrix in the hyperspectral radiance data into the target reflectance conversion function to obtain the two-dimensional reflectance matrix corresponding to each spectral band.

[0137] In this embodiment of the invention, for each spectral band in the hyperspectral radiance data, the corresponding two-dimensional radiance matrix is ​​extracted. The gray value (DN value) of each pixel in the matrix is ​​substituted into the target reflectance conversion function constructed in step S33. Combined with the illumination correction coefficient corresponding to the band and the standard target reference reflectance, the reflectance value of each pixel is calculated. Through band-by-band and pixel-by-pixel conversion calculation, the influence of atmosphere and illumination is eliminated, and the two-dimensional reflectance matrix corresponding to each spectral band is obtained. The row and column dimensions of each matrix are consistent with the original image, and the matrix elements are the spectral reflectance values ​​of the corresponding pixels in that band.

[0138] S35. Reassemble the two-dimensional reflectance matrices according to the spectral band order to obtain the hyperspectral reflectance data of the area to be detected.

[0139] In this embodiment of the invention, the two-dimensional reflectance matrices corresponding to all spectral bands are reorganized according to the original spectral band order to restore the three-dimensional data cube structure consistent with the hyperspectral radiance data, thereby constructing the hyperspectral reflectance data of the area to be detected. In this data cube, the value of each spatial pixel in each spectral band is the reflectance value after illumination correction and atmospheric disturbance correction, which has a clear physical meaning and can be directly used for spectral feature extraction and identification analysis of karst landforms, providing a reliable spectral data foundation for the accurate identification of karst landforms.

[0140] Ground feature spectral data is ground reference data containing measured spectral reflectance of ground features such as limestone and vegetation. Standard target measured data is measured spectral and image grayscale data of the standard target. Preset standard target reference reflectance is the inherent reference reflectance value of the standard target. Illumination correction coefficient is a coefficient that corrects the influence of illumination changes during flight. Preset basic reflectance conversion model is an initially constructed radiance-to-reflectance basic model. Target reflectance conversion function is a precise reflectance conversion mathematical function after illumination correction. Two-dimensional reflectance matrix is ​​a two-dimensional matrix composed of single-band pixel reflectance values.

[0141] Step 205: Using ground-measured benchmark data and topographic spatial reference data, perform geometric correction on the hyperspectral reflectance data to obtain corrected hyperspectral reflectance data.

[0142] Furthermore, the topographic spatial reference data includes DEM data and orthophoto data, and the ground-measured reference data includes control point data. Step 205 may include the following sub-steps:

[0143] S41. Using orthophoto data as a spatial reference, determine the initial pixel coordinates corresponding to each two-dimensional reflectance matrix within the hyperspectral reflectance data.

[0144] In this embodiment of the invention, the POS (Position and Orientation System) data acquired synchronously during the hyperspectral data acquisition process is first processed. The original POS data is smoothed, and then the WGS84 (World Geodetic System 1984) frame coordinates and attitude parameters of each line of hyperspectral image data at the time of exposure are precisely calculated using Precise Point Positioning (PPP). The calculated POS data is then interpolated to obtain the exterior orientation information corresponding to each scan line of the hyperspectral image, including position and attitude parameters. Subsequently, a rotation matrix is ​​established to perform coordinate transformation, converting the navigation coordinate system to the projected coordinate system. The coordinate transformation process is as follows: imaging coordinate system (m) → navigation coordinate system (g) → IMU (Inertial Measurement Unit) coordinate system (b) → sensor coordinate system (c) → image space coordinate system (i). The expression for the rotation matrix transformation from the imaging coordinate system to the image space coordinate system is:

[0145]

[0146] In the formula, This is the rotation matrix from the imaging coordinate system to the image space coordinate system; This is the transformation matrix from the navigation coordinate system to the imaging coordinate system; Let be the transformation matrix from the IMU coordinate system to the navigation coordinate system, where IMU attitude angle; This is the transformation matrix from the sensor coordinate system to the IMU coordinate system; This is the transformation matrix from the image space coordinate system to the sensor coordinate system; The three exterior orientation elements are attitude parameters, which can be obtained through matrix transformation.

[0147] Based on this, combined with the coordinate position offset in the IMU coordinate system ( , , The position parameters are corrected, and the coordinate transformation formula is as follows:

[0148]

[0149] In the formula, , , The corrected position parameters for the exterior orientation elements; , , These are the initial position parameters in the navigation coordinate system; This is the transformation matrix from the navigation coordinate system to the imaging coordinate system; This is the transformation matrix from the IMU coordinate system to the navigation coordinate system; , , This represents the coordinate position offset in the IMU coordinate system. Based on the corrected exterior orientation information and rotation matrix, and using orthophoto data as a spatial reference, the pixel positions of each two-dimensional reflectance matrix in the hyperspectral reflectance data are matched to determine the initial pixel coordinates corresponding to each pixel.

[0150] S42. Calculate the actual ground elevation corresponding to each pixel in the area to be detected based on the DEM data.

[0151] In this embodiment of the invention, the initial pixel coordinates of each pixel obtained in S41 are... Spatial raster matching and bilinear interpolation are performed on the DEM data to obtain the actual ground elevation corresponding to the pixel. The calculation formula is:

[0152]

[0153] In the formula, For pixels The corresponding actual ground elevation, in meters; For DEM data and pixels Elevation values ​​of four adjacent grid points; The bilinear interpolation weights are the values ​​corresponding to each adjacent grid point. and These are the row and column index parameters that mark the locations of adjacent DEM raster points of the target cell; both values ​​are 0 or 1. This indicates the row direction (vertical) index of the target cell (x, y) within the DEM raster cell. The column direction (horizontal) index is represented by the four combinations of the two ((0,0), (0,1), (1,0), (1,1)) which exactly traverse the four adjacent grid points of the target cell (corresponding to the top left, top right, bottom left, and bottom right corners of the grid cell respectively), and are used to locate the adjacent grid points participating in the elevation interpolation calculation.

[0154] S43. Based on the reference plane elevation and the actual ground elevation corresponding to the orthophoto data, calculate the position offset of each pixel as an elevation correction term.

[0155] In this embodiment of the invention, the reference plane elevation corresponding to the orthophoto data is used. For reference, first calculate the elevation difference between the actual ground elevation of the pixel and the elevation of the reference plane. Then, by combining the imaging geometry, the elevation difference is converted into the pixel plane position offset. , This is used as an elevation correction term. The formula for calculating the elevation difference is:

[0156]

[0157] In the formula, This represents the pixel elevation difference, in meters (m). This represents the actual ground elevation corresponding to a pixel, in meters. The elevation of the reference plane for the orthophoto is in meters.

[0158] Formula for calculating pixel position offset:

[0159]

[0160] In the formula, , These are the positional offsets of the pixels in the row and column directions (elevation correction items), respectively. This is the sensor focal length, in mm; The relative flight altitude of the drone is expressed in meters (m). , These are the initial pixel coordinates.

[0161] S44. Using the measured plane coordinates of the control point data as the true value benchmark, and combining the initial pixel coordinates with the elevation correction term, construct a geometric correction mathematical model.

[0162] In this embodiment of the invention, a geometric correction mathematical model is constructed using collinearity condition equations. This model associates image point coordinates, ground point coordinates, and exterior orientation elements, as shown in the following formula:

[0163]

[0164] In the formula, , These are the planar coordinates of the ground point corresponding to the pixel; , , For the position parameters of the exterior orientation element; This represents the actual ground elevation corresponding to the pixel. , The initial pixel coordinates; , This is an elevation correction item; , , , , , , , , The direction cosine is related to the attitude parameters; The focal length is the sensor's focal length. During model construction, the measured planar coordinates of the control point data are used as the true values. Combined with the initial pixel coordinates, elevation correction terms, and exterior orientation parameters of each pixel, these are substituted into the collinearity condition equation to form a geometric correction mathematical model containing model coefficients.

[0165] S45. Solve the model coefficients of the geometric correction mathematical model using the least squares method to obtain the target geometric correction model.

[0166] In this embodiment of the invention, the model coefficients of the geometric correction mathematical model (collinearity condition equation) specifically include:

[0167] Image exterior orientation parameters: , , (The three-dimensional spatial coordinates of the sensor projection center);

[0168] Image exterior orientation and attitude parameters: (Roll angle, pitch angle, yaw angle);

[0169] Direction cosine coefficients derived from attitude parameters: , , , , , , , , .

[0170] Measured plane coordinates of control point data Using the true value as a benchmark, a residual equation is constructed between the model's calculated coordinates and the measured coordinates. The sum of squared residuals is used as the optimization objective, and the least squares method is used to iteratively solve for the undetermined model coefficients. Through multiple iterations, the coordinate residuals are eliminated, and the optimal combination of model coefficients is obtained by convergence. The optimal coefficients are then substituted into the collinearity condition equation to construct a target geometric correction model adapted to the area to be detected. This model can accurately correct the geometric distortion caused by terrain undulations and attitude errors.

[0171] S46. Substitute the hyperspectral reflectance data into the target geometric correction model and perform spatial transformation to obtain the corrected hyperspectral reflectance data.

[0172] In this embodiment of the invention, the hyperspectral reflectance data is a three-dimensional reflectance data cube sorted by spectral bands, containing two-dimensional reflectance matrices corresponding to multiple spectral bands. Each matrix element is the original spectral reflectance value of a pixel, with only spatial coordinate distortion. The two-dimensional reflectance matrices of each spectral band are successively substituted into the target geometric correction model. The initial pixel coordinates of each pixel within the matrix are combined with an elevation correction term to perform spatial coordinate transformation, obtaining the accurate ground geographic coordinates corresponding to the pixel while preserving the original spectral reflectance value of the pixel. After completing the spatial coordinate correction of all pixels, the corrected two-dimensional reflectance matrices of each band are reassembled according to the original spectral band order to restore the three-dimensional data cube structure, ultimately obtaining the corrected hyperspectral reflectance data. This corrected data possesses both accurate geographic spatial coordinates and faithful spectral reflectance information. The spatial position of each pixel is fully registered with DEM (Digital Elevation Model) data and orthophoto data, and the spectral reflectance value is not distorted, allowing it to be directly used for subsequent accurate identification of karst landforms.

[0173] DEM data, or Digital Elevation Model data, is three-dimensional topographic data characterizing the elevation of a region's ground. Orthophoto data is precise planar image data that eliminates topographic and projection distortions and serves as a spatial reference. Control point data is reference data containing the measured coordinates of control points and their corresponding coordinates in the image. Initial pixel coordinates are the initial spatial coordinates of hyperspectral pixels determined based on the orthophoto. Actual ground elevation is the true elevation of the ground point corresponding to the pixel. Reference plane elevation is the elevation of the reference plane corresponding to the orthophoto. Position offset is the pixel coordinate deviation caused by topographic undulations. Elevation correction terms are supplementary data. The coordinate correction parameters compensate for the influence of terrain undulations. The measured plane coordinates are the precise plane coordinates measured in the field by the image control points. The true reference is the measured coordinates of the image control points used as the correction standard. The geometric correction mathematical model is a mathematical model that describes the mapping relationship between the pixel image coordinates and the real ground coordinates. The least squares method is an optimization method that minimizes the sum of squared residuals to solve for the optimal parameters of the model. The model coefficients are the exterior orientation parameters and direction cosine coefficients in the geometric correction model. The target geometric correction model is the precise coordinate correction model obtained after solving for the coefficients. Spatial transformation is the process of converting the initial coordinates of the pixels to the real ground coordinates through the correction model.

[0174] Step 206: Karst landform identification is performed using calibrated hyperspectral reflectance data, topographic spatial reference data, and ground-measured benchmark data to obtain karst landform identification results for the area to be explored.

[0175] Furthermore, step 206 may include the following sub-steps:

[0176] S51. Preprocess the corrected hyperspectral reflectance data, including interference masking, dimensionality reduction, and noise suppression.

[0177] In this embodiment of the invention, the corrected hyperspectral reflectance data is first preprocessed. The interference masking process mainly removes vegetation, water bodies, and shadowed areas, with a focus on exposed bedrock areas to avoid non-target features interfering with subsequent karst lithology identification. The dimensionality reduction and noise suppression processes retain effective spectral information and improve the signal-to-noise ratio of the data while compressing the data dimension, reducing the data volume, and improving the subsequent processing speed, providing clean and efficient preprocessed data for karst feature identification.

[0178] S52. Extract the measured limestone spectrum from the spectral data of ground features.

[0179] In this embodiment of the invention, the characteristic absorption bands of carbonate rocks (taking limestone as an example) are mainly composed of Induced by the group, located in the 2300-2350nm band, the center wavelength of its strong absorption band is usually around 2340nm; the spectral curve of the measured limestone in the field is extracted from the ground-measured reference data as a reference endmember spectrum. This spectrum contains the typical reflectance characteristics of limestone in the visible-shortwave infrared band, providing a reliable benchmark for subsequent spectral angle matching.

[0180] S53. Based on the spectral angle mapping method, using the measured limestone spectrum as the end element and setting an angle threshold, calculate the angle between the spectral vector of each pixel in the preprocessed corrected hyperspectral reflectance data and the limestone reference spectral vector.

[0181] In this embodiment of the invention, the Spectral Angle Mapper (SAM) method is used, with the measured limestone spectrum extracted by S52 as the endmember, and the maximum angle threshold is set according to the accuracy requirements for karst landform identification. (Typically 3°~10°, can be adjusted according to the actual situation in the area), calculate the angle between the spectral vector of each pixel in the preprocessed corrected hyperspectral reflectance data and the reference spectral vector of the limestone. The calculation formula is:

[0182]

[0183] In the formula, The angle between the pixel spectral vector and the limestone reference spectral vector is expressed in degrees. To correct the first hyperspectral reflectance data after preprocessing Pixel reflectance values ​​for the band; The measured reference spectrum of limestone was obtained in the [number]th [year]. The reflectivity value of the band; The total number of spectral bands involved in the calculation; the smaller the angle, the higher the similarity between the pixel spectrum and the limestone reference spectrum.

[0184] S54. Mark the pixels with included angles below the angle threshold as carbonate rock target pixels to obtain carbonate rock extraction results.

[0185] In this embodiment of the invention, the calculated spectral angle of each pixel is compared with a preset angle threshold. Comparison, if the pixel spectral angle If the spectral characteristics of a pixel are highly similar to those of a limestone reference spectrum, it is marked as a carbonate rock target pixel; otherwise, it is marked as a non-carbonate rock pixel, thus generating the initial extraction results of the carbonate rock distribution in the area to be detected.

[0186] S55. Perform linear spectral unmixing on the carbonate rock extraction results to determine the carbonate rock region.

[0187] In this embodiment of the invention, since karst areas often contain mixed spectra of limestone and other associated minerals (such as argillaceous and siliceous minerals), based on the initial extraction results of carbonate rocks obtained by the spectral angle mapping method, a linear spectral demixing model is used to decompose the mixed pixels, separate the limestone endmembers and associated mineral endmembers, correct the misjudgments and omissions caused by the mixed spectra in the initial extraction results, and obtain a more accurate range of carbonate rock areas.

[0188] S56. Calculate the terrain factor corresponding to each pixel in the area to be detected based on the DEM data.

[0189] In this embodiment of the invention, based on high-precision DEM data, multiple terrain factors corresponding to each pixel within the area to be detected are calculated, specifically including:

[0190] Elevation features: Extract the absolute elevation value of each pixel and use the stratified characteristics of the elevation distribution of karst peak forest depressions (higher elevation in the peak cluster area and lower elevation at the bottom of the depression) as the basis for karst landform identification.

[0191] Slope characteristics: Calculate the ground slope of each pixel to distinguish the steepness of the terrain between peak forests and depressions;

[0192] Slope aspect characteristics: Calculate the slope aspect of each pixel to provide a reference for depression identification;

[0193] Curvature characteristics: Calculate the planar curvature and the profile curvature. Peak forests (or ridges) usually exhibit a convex shape with a positive planar curvature (dispersion of water flow) and a negative profile curvature (accelerated erosion of water flow). Depressions usually exhibit a concave shape with a negative planar curvature (convergence of water flow) and a positive profile curvature (decelerated deposition of water flow).

[0194] Surface relative humidity index: As a hydrological feature derived from topography, this index reflects the impact of topography on water accumulation. The value is relatively high at the bottom of depressions and relatively low in peak cluster areas. It can help distinguish between peak forests and depression areas and be included in the category of topographic factors.

[0195] S57. Based on topographic factors, select topographic features and determine candidate regions for peak forests and depressions respectively.

[0196] In this embodiment of the invention, based on various calculated topographic factors, targeted screening conditions are set to extract candidate karst landform areas:

[0197] Peak forest candidate area screening: Areas with a slope >30° and a planar curvature >0.2 were screened out, and the low surface relative humidity index was also considered to determine the peak forest candidate area set;

[0198] Depression candidate area screening: Closed areas with a slope of <10° and a plane curvature of <-0.15 were screened out. At the same time, the high surface relative humidity index was combined to determine the depression candidate area set.

[0199] The above screening criteria match the typical topographic features of karst peak forest depressions, which can effectively eliminate interference from other topographic types and narrow the identification range.

[0200] S58. Spatially overlay the candidate areas of peak forests and depressions with the carbonate rock areas to obtain the karst landform identification results of the area to be explored.

[0201] In this embodiment of the invention, the candidate areas for peak forests and depressions are spatially overlaid with the carbonate rock areas obtained in S55. Only areas that simultaneously meet the criteria of being carbonate rocks and having topography that conforms to the characteristics of peak forests or depressions are retained, while pseudo-topographic candidate areas that are not carbonate rocks are eliminated. After overlay, the karst landform types such as peak forests and depressions in the area to be detected are classified and labeled to generate karst landform identification results containing the distribution range and boundary information of each karst landform type. Finally, the identification results can be verified on-site by combining aerial photography data to verify the distribution location and morphological characteristics of karst landforms, ensuring the accuracy and reliability of the identification results.

[0202] Preprocessing is the initial data processing to eliminate interference and improve data quality. Interference masking is the process of removing non-target features such as vegetation and water bodies. Dimensionality reduction and noise suppression are the processes of compressing data dimensions, suppressing random noise, and improving the signal-to-noise ratio. Measured limestone spectra are the spectral curves of limestone measured in the field. Spectral angle mapping is an identification method based on the angle between spectral vectors to determine the spectral similarity of ground features. Endmembers are the pure ground feature spectra used as the identification benchmark. Angle thresholds are the critical angle values ​​for determining spectral similarity. Pixel spectral vectors are feature vectors composed of the reflectance of each band of a single pixel. Limestone reference spectral vectors are the benchmark vectors composed of the measured reflectance of limestone. It is the spatial angle between two types of spectral vectors. The target pixel of carbonate rock is the pixel to be judged by spectral matching of limestone. The carbonate rock extraction result is a preliminary set of carbonate rock pixels. Linear spectral unmixing is a method to decompose the mixed pixel spectrum and correct the identification error. The carbonate rock region is the true distribution range of carbonate rock after unmixing correction. The topographic factor is the topographic feature parameters such as elevation, slope, curvature, and surface humidity calculated based on DEM. Topographic feature screening is the process of screening candidate regions based on karst topographic features. The peak forest and depression candidate region set is a preliminary set of peak forest and depression candidate regions. Spatial overlay is an analysis method that matches spatial data and retains overlapping areas.

[0203] Furthermore, step 206 may also include the following sub-steps:

[0204] S59. Morphological screening of the candidate regions of peak forests and depressions is performed to obtain a second-stage candidate region set of peak forests and depressions.

[0205] In this embodiment of the invention, morphological analysis is performed on the candidate region set of karst peaks and depressions obtained in S57. Screening conditions are set in combination with the typical morphological characteristics of karst peaks and depressions to further eliminate false candidate regions:

[0206] Morphological screening of peak forests: The area of ​​the candidate peak forest area must be greater than 500m². 2 Furthermore, the elevation difference within the region is greater than 30m, ensuring that the candidate area possesses the scale and elevation difference characteristics of karst peak forests, and eliminating pseudo-peak forests caused by small protrusions or artificial features.

[0207] Morphological screening of depressions: The closure index of candidate depression areas is required to be greater than 0.8. The formula for calculating the closure index is as follows:

[0208]

[0209] In the formula, The closure index of the candidate depression area; The perimeter of the candidate depression area is expressed in meters (m). The area of ​​the candidate depression is expressed in m². 2This index reflects the degree of morphological enclosure of a depression. A higher value indicates that the depression is closer to a circle in shape and has better water catchment enclosure. It can effectively identify the typical depression characteristics of karst depressions and eliminate pseudo-depressions that do not meet the conditions for water catchment. Through the above screening criteria, a set of candidate areas for peak forests and depressions was obtained after secondary screening.

[0210] S510. The candidate areas of peak forests and depressions selected in the second screening are spatially superimposed with the carbonate rock areas to obtain new karst landform identification results for the area to be explored.

[0211] In this embodiment of the invention, hyperspectral reflectance data and DEM data generated by lidar are spatially registered to ensure accurate correspondence at the pixel level. First, the candidate areas of peak forests and depressions selected in S59 are spatially overlaid with the carbonate rock areas determined in S55. Only areas that simultaneously meet the requirements of carbonate rock lithology and topographic and morphological features that conform to the peak forest or depression requirements are retained, while false topographic candidate areas of non-carbonate rock areas are eliminated. Subsequently, the relative surface water content is retrieved from the hyperspectral data to verify the overlaid candidate areas: the relative surface water content at the bottom of the depression is higher, while the relative surface water content in the peak forest area is lower, thereby further verifying the karst landform attributes of the area and eliminating false identification results that do not match the hydrological characteristics. Finally, the DEM and orthophoto data are comprehensively interpreted to generate a karst landform identification result map of the area to be explored, which includes the distribution range and boundary information of peak forests and depressions. This result significantly reduces the misjudgment rate of traditional single-factor identification methods and achieves accurate identification of karst landforms.

[0212] Morphological screening is a screening process that eliminates false candidate areas based on the morphological characteristics of peak forests and depressions. The set of candidate areas for peak forests and depressions after secondary screening is a set of high-precision candidate areas after morphological screening. The new karst landform identification result is a high-precision identification result obtained by superimposing carbonate rock areas after secondary screening.

[0213] Furthermore, step 206 may also include the following sub-steps:

[0214] S511. Based on the corrected hyperspectral reflectance data, the relative content of surface water in the area to be detected is retrieved.

[0215] In this embodiment of the invention, the sensitive spectral response of hyperspectral data to water bodies and moist surfaces is utilized. Based on corrected hyperspectral reflectance data, spectral bands sensitive to water absorption characteristics (such as the 1300-1500 nm and 1800-2000 nm bands) are selected, and the relative surface water content data of the area to be detected is obtained by inversion using the spectral index method. Specifically, a normalized water index (NDWI) type spectral inversion model is constructed, and the calculation formula is as follows:

[0216]

[0217] In the formula, This represents the relative surface water content value corresponding to each pixel. The corrected reflectance is for the near-infrared band (e.g., 800-900nm). It is the corrected reflectance for short-wave infrared bands (such as 1300-1500nm or 1800-2000nm); the larger the value, the higher the surface moisture content, and vice versa, which can effectively distinguish between dry peak forest areas and humid depression areas.

[0218] S512. Extract the relative surface water content values ​​of each candidate region in the secondary screening of peak forest and depression candidate region set from the relative surface water content data.

[0219] In this embodiment of the invention, based on the spatial boundary vector data of the candidate peak forest and depression areas obtained by secondary screening in S59, spatial pixel-level extraction and statistics are performed on the relative surface water content data; the relative surface water content values ​​of all pixels in each candidate peak forest area and candidate depression area are extracted respectively, and the average value in the area is calculated as the overall relative surface water content feature value of the candidate area.

[0220] S513. Determine whether the relative surface water content values ​​corresponding to each candidate area meet the preset hydrological judgment conditions for association.

[0221] In this embodiment of the invention, based on the hydrological characteristics of karst landforms, corresponding hydrological judgment conditions are preset:

[0222] Criteria for determining candidate areas for karst peak forests: regional average relative surface water content. Less than the preset threshold (Usually a value of 0.2~0.35), indicating dry surface, consistent with the characteristics of peak forest topography;

[0223] Criteria for determining candidate depression areas: regional average relative surface water content. Greater than the preset threshold (Usually a value of 0.5~0.7) indicates that the surface is moist and water is easy to accumulate, which is consistent with the characteristics of depression terrain; the average relative content of surface water in each candidate area is compared with the corresponding threshold to determine whether the hydrological judgment conditions are met.

[0224] S514. Retain candidate areas that meet the preset hydrological judgment conditions to obtain a set of candidate areas for peak forests and depressions after three rounds of screening.

[0225] In this embodiment of the invention, the judgment results of S513 are screened: all candidate peak forest areas whose relative surface water content meets the hydrological judgment conditions of peak forests are retained, and pseudo peak forest areas with excessively high water content and not meeting the dry characteristics of peak forests are removed; all candidate depression areas whose relative surface water content meets the hydrological judgment conditions of depressions are retained, and pseudo depression areas with excessively low water content and not meeting the water catchment characteristics of depressions are removed; through three screenings of hydrological conditions, a set of candidate peak forest and depression areas with higher accuracy is obtained, providing a reliable spatial data foundation for the final karst landform identification.

[0226] S515. The candidate areas of peak forests and depressions selected in three screenings are spatially superimposed with the carbonate rock areas to obtain new karst landform identification results for the area to be explored.

[0227] In this embodiment of the invention, the candidate area set of peak forests and depressions obtained from the three screenings in S514 is spatially overlaid with the carbonate rock lithology area determined in S55. Only areas that simultaneously meet the requirements of carbonate rock lithology, topographic features, and hydrological conditions are retained, and areas that are misjudged due to a single factor such as lithology, topography, or hydrology are completely eliminated. After the overlay is completed, the peak forests and depressions in the area to be detected are delineated and classified. The data are then combined with DEM and orthophoto data for comprehensive verification. Finally, the karst landform identification result of the area to be detected is generated. This result has multiple constraints of lithology, topography, and hydrology, which greatly improves the accuracy and reliability of karst landform identification.

[0228] The relative surface water content data is the data representing the surface moisture level of the region obtained by inversion. The relative surface water content value is the average surface water content of the pixels in the candidate area. The preset hydrological judgment conditions are the verification standards set based on the dry peak forest and moist depression characteristics. The candidate area set of peak forest and depression after three screenings is the highest precision candidate area set after hydrological screening.

[0229] It is worth mentioning that this invention employs a progressive three-stage screening mechanism—topographic feature screening, morphological screening, and hydrological feature screening—introducing relative surface water content as a third hydrological constraint on top of the first two topographic and morphological constraints. Starting from the essential hydrological evolution law of karst landforms, this invention solves the problem that traditional single-factor identification methods relying solely on topography or spectrum easily misidentify non-karst topography as the target landform. This significantly reduces the false target misjudgment rate, substantially improves the accuracy of karst landform identification and the geological reliability of the results, and adapts to the identification needs of complex and fractured karst environments. It avoids the problem of insufficient identification accuracy of traditional methods under fractured terrain, reduces the workload of subsequent field verification, improves the overall engineering application efficiency of karst landform surveys, and reduces implementation costs.

[0230] Please see Figure 4 , Figure 4This is a structural block diagram of a karst landform identification system based on UAV hyperspectral imaging provided in Embodiment 3 of the present invention.

[0231] This invention provides a karst landform identification system based on UAV hyperspectral imaging, comprising:

[0232] The acquisition module 401 is used to acquire range data of the area to be detected and to conduct survey and layout, and to determine the take-off and landing points and flight route map of the UAV based on the survey and layout results;

[0233] The acquisition module 402 is used to acquire hyperspectral data and ground-measured reference data of the area to be detected based on the take-off and landing point of the UAV and the flight route map, and simultaneously acquire terrain spatial reference data.

[0234] The correction module 403 is used to perform radiometric geometric joint correction on hyperspectral data based on ground-measured benchmark data and topographic spatial reference data to obtain corrected hyperspectral reflectance data.

[0235] The identification module 404 is used to identify karst landforms by using corrected hyperspectral reflectance data, topographic spatial reference data and ground measured benchmark data, and to obtain the karst landform identification results of the area to be detected.

[0236] Since the above is a system corresponding to a karst landform identification method based on UAV hyperspectral imaging, its implementation principle is the same as that of a karst landform identification method based on UAV hyperspectral imaging. For the sake of convenience and brevity, those skilled in the art can clearly understand that the specific working process of the system and modules described above can be referred to the corresponding process in the aforementioned method embodiments, and will not be repeated here.

[0237] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0238] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces, or indirect coupling or communication connection between apparatuses or units, and may be electrical, mechanical, or other forms.

[0239] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0240] Furthermore, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0241] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0242] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying karst landforms based on UAV hyperspectral imaging, characterized in that, include: Acquire the range data of the area to be explored and conduct survey and deployment. Based on the survey and deployment results, determine the take-off and landing points and flight route map of the UAV. Based on the take-off and landing point of the UAV, hyperspectral data and ground-measured reference data of the area to be detected are collected according to the flight route map, and terrain spatial reference data are acquired simultaneously. Based on the ground-measured reference data and the topographic spatial reference data, the hyperspectral data is subjected to radiometric geometric joint correction to obtain corrected hyperspectral reflectance data. Karst landform identification is performed using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground measured benchmark data to obtain the karst landform identification results for the area to be detected.

2. The karst landform identification method based on UAV hyperspectral imaging according to claim 1, characterized in that, The process of acquiring range data of the area to be detected and conducting surveying and deployment, and determining the UAV take-off and landing points and flight path map based on the surveying and deployment results, includes: Obtain the range data of the area to be detected; Based on the range data, obtain the edge survey information of the area to be detected; Based on the on-site edge survey information, the location layout information of ground measurement base stations and image control points is determined within the area to be detected; The take-off and landing points of the UAV are determined by using the on-site edge survey information and the point layout information as constraints. A flight path map is generated based on the range data and the on-site survey information of the edge.

3. The karst landform identification method based on UAV hyperspectral imaging according to claim 1 or 2, characterized in that, The step of performing radiometric geometric joint correction on the hyperspectral data based on the ground-measured reference data and the topographic spatial reference data to obtain corrected hyperspectral reflectance data includes: The hyperspectral data were radiometrically calibrated to obtain hyperspectral radiance data; Based on the ground-measured benchmark data, the hyperspectral radiance data is converted to spectral reflectance to obtain hyperspectral reflectance data. Using the ground-measured benchmark data and the topographic spatial reference data, the hyperspectral reflectance data is geometrically corrected to obtain corrected hyperspectral reflectance data.

4. The karst landform identification method based on UAV hyperspectral imaging according to claim 3, characterized in that, The step of radiometrically calibrating the hyperspectral data to obtain hyperspectral radiance data includes: Extract the hyperspectral raw data cube from the hyperspectral data; Extract the original DN value matrix corresponding to each spectral band from the hyperspectral raw data cube; Based on the original DN value matrix and the radiation intensity of the preset standard plate, calculate the radiation calibration coefficient corresponding to each spectral band; A radiometric calibration transformation function is constructed based on the radiometric calibration coefficients, and the original DN value matrices are substituted into the function to obtain the two-dimensional radiance matrix corresponding to each spectral band. The two-dimensional radiance matrices are reorganized according to the spectral band order to obtain the hyperspectral radiance data of the region to be detected.

5. The karst landform identification method based on UAV hyperspectral imaging according to claim 4, characterized in that, The ground-based measured reference data includes ground object spectral data. The process of converting the hyperspectral radiance data to spectral reflectance based on the ground-based measured reference data to obtain hyperspectral reflectance data includes: Extract measured data of standard targets from the spectral data of the ground features; Using the measured data of the standard target and the preset standard target reference reflectance, the illumination correction coefficient corresponding to each of the spectral bands is calculated. The target reflectivity conversion function is obtained by correcting the preset basic reflectivity conversion model using the illumination correction coefficient. Substitute each two-dimensional radiance matrix in the hyperspectral radiance data into the target reflectance conversion function to obtain the two-dimensional reflectance matrix corresponding to each spectral band. The two-dimensional reflectance matrices are reorganized according to the spectral band order to obtain the hyperspectral reflectance data of the region to be detected.

6. The karst landform identification method based on UAV hyperspectral imaging according to claim 5, characterized in that, The topographic spatial reference data includes DEM data and orthophoto data, and the ground-measured reference data includes image control point data. The step of using the ground-measured reference data and the topographic spatial reference data to perform geometric correction on the hyperspectral reflectance data to obtain corrected hyperspectral reflectance data includes: Using the orthophoto data as a spatial reference, the initial pixel coordinates corresponding to each two-dimensional reflectance matrix within the hyperspectral reflectance data are determined; Calculate the actual ground elevation corresponding to each pixel in the area to be detected based on the DEM data; Based on the reference plane elevation corresponding to the orthophoto data and the actual ground elevation, the positional offset of each pixel is calculated as an elevation correction term. Using the measured planar coordinates of the control point data as the truth benchmark, and combining the initial pixel coordinates with the elevation correction term, a geometric correction mathematical model is constructed. The model coefficients of the geometric correction mathematical model are solved using the least squares method to obtain the target geometric correction model; The hyperspectral reflectance data is substituted into the target geometric correction model for spatial transformation to obtain the corrected hyperspectral reflectance data.

7. The karst landform identification method based on UAV hyperspectral imaging according to claim 6, characterized in that, The process of identifying karst landforms using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground-measured benchmark data, to obtain karst landform identification results for the area to be detected, includes: The corrected hyperspectral reflectance data is preprocessed, including interference masking, dimensionality reduction, and noise suppression. Extract the measured limestone spectrum from the aforementioned ground feature spectral data; Based on the spectral angle mapping method, using the measured limestone spectrum as the end element and setting an angle threshold, the angle between the spectral vector of each pixel in the preprocessed corrected hyperspectral reflectance data and the limestone reference spectral vector is calculated. Pixels with included angles below the specified angle threshold are marked as carbonate rock target pixels to obtain carbonate rock extraction results. Linear spectral unmixing was performed on the carbonate rock extraction results to determine the carbonate rock region; Calculate the terrain factor corresponding to each pixel in the area to be detected based on the DEM data; Based on the aforementioned topographic factors, topographic features are screened to determine candidate regions for peak forests and depressions, respectively. The candidate areas of peak forests and depressions are spatially superimposed with the carbonate rock areas to obtain the karst landform identification results of the area to be detected.

8. The karst landform identification method based on UAV hyperspectral imaging according to claim 7, characterized in that, Also includes: Morphological screening was performed on the candidate region set of peak forests and depressions to obtain a second-screened candidate region set of peak forests and depressions. The candidate areas of peak forests and depressions selected in the secondary screening are spatially superimposed with the carbonate rock area to obtain new karst landform identification results for the area to be explored.

9. The karst landform identification method based on UAV hyperspectral imaging according to claim 8, characterized in that, Also includes: Based on the corrected hyperspectral reflectance data, the relative surface water content data of the area to be detected is retrieved. Extract the relative surface water content values ​​of each candidate region within the secondary screening set of peak forest and depression candidate regions from the relative surface water content data. Determine whether the relative surface water content values ​​corresponding to each candidate region meet the associated preset hydrological determination conditions; Candidate areas that meet the preset hydrological judgment conditions are retained to obtain a set of candidate areas for peak forests and depressions after three rounds of screening; The candidate areas of peak forests and depressions selected through the three screenings are spatially superimposed with the carbonate rock area to obtain new karst landform identification results for the area to be explored.

10. A karst landform identification system based on UAV hyperspectral imaging, characterized in that, include: The acquisition module is used to acquire range data of the area to be detected and to conduct survey and deployment. Based on the survey and deployment results, the take-off and landing points and flight route map of the UAV are determined. The data acquisition module is used to acquire hyperspectral data and ground-measured reference data of the area to be detected based on the take-off and landing point of the UAV and according to the flight route map, and simultaneously acquire terrain spatial reference data. The correction module is used to perform radiometric geometric joint correction on the hyperspectral data based on the ground measured reference data and the topographic spatial reference data to obtain corrected hyperspectral reflectance data. The identification module is used to identify karst landforms using the corrected hyperspectral reflectance data, the topographic spatial reference data, and the ground measured benchmark data, and to obtain the karst landform identification results of the area to be detected.

Images