Karst area linear structure identification method and system based on multi-source remote sensing

Abstract

The invention provides a karst area linear structure identification method and system based on multi-source remote sensing, and the method comprises the steps: firstly obtaining an optical remote sensing image, radar remote sensing data and terrain elevation model data of a target area, constructing a target remote sensing data set, carrying out the ground surface coverage feature enhancement processing of the target remote sensing data set, and carrying out the recognition of a karst area linear structure. The method comprises the following steps: generating an enhanced multi-source feature set containing features such as a karst area earth surface coverage type, topographic relief and radar scattering characteristics, and then calling a pre-trained linear structure recognition model to recognize the enhanced set to obtain an initial result containing linear structure spatial distribution data and a structure type label; and finally, performing spatial topological relation verification processing on the initial result based on the geological background data of the target area, and generating an optimized linear structure identification result for outputting the fault and fracture spatial position and the structure confidence coefficient parameter of the karst area, thereby effectively identifying the linear structure of the karst area.

Description

Technical Field

[0001] The present invention relates to the field of artificial intelligence technology, and in particular, to a method and system for identifying linear structures in karst areas based on multi-source remote sensing. Background Technique

[0002] In the fields of geological exploration and geological disaster warning, the identification of linear structures in karst areas has always been a highly challenging task. The geological conditions in karst areas are complex, with diverse surface coverage types, large topographic fluctuations, and well-developed underground karst caves, fissures and other geological structures. These characteristics make it difficult for traditional linear structure identification methods to face many difficulties.

[0003] In the related technologies, the mutual correlation and complementarity between different features are ignored, and the multi-dimensional information of linear structures in karst areas cannot be fully exploited. Moreover, the related technologies mostly adopt methods based on rules or simple statistical models, and the above methods have poor adaptability to the complex and changeable geological conditions in karst areas and are difficult to process high-dimensional and non-linear feature data. Although deep learning models have achieved remarkable results in the field of image recognition in recent years, their application in the identification of linear structures in karst areas is relatively less, and there is a lack of a dedicated model for multi-source feature fusion and geological background constraint in karst areas. Summary of the Invention

[0004] In view of the problems mentioned above, in combination with the first aspect of the present invention, an embodiment of the present invention provides a method for identifying linear structures in karst areas based on multi-source remote sensing, and the method includes:

[0005] Obtain optical remote sensing images, radar remote sensing data and terrain elevation model data of the target area, and construct a target remote sensing data set;

[0006] Perform surface coverage feature enhancement processing on the target remote sensing data set to generate an enhanced multi-source feature set; the enhanced multi-source feature set includes surface coverage type features, terrain undulation features and radar scattering characteristic features in karst areas;

[0007] Call a pre-trained linear structure identification model to perform linear structure identification on the enhanced multi-source feature set to generate an initial linear structure identification result; the initial linear structure identification result includes spatial distribution data of linear structures and structure type labels;

[0008] Perform spatial topological relationship verification processing on the initial linear structure identification result based on the geological background data of the target area to generate an optimized linear structure identification result; the optimized linear structure identification result is used to output the spatial positions and structure confidence parameters of faults and fissures in karst areas.

[0009] In another aspect, an embodiment of the present invention further provides a linear structure recognition system for karst areas based on multi-source remote sensing, including a processor and a machine-readable storage medium. The machine-readable storage medium is connected to the processor. The machine-readable storage medium is used to store programs, instructions or codes, and the processor is used to execute the programs, instructions or codes in the machine-readable storage medium to implement the above method.

[0010] Based on the above aspects, in the embodiment of the present application, by acquiring optical remote sensing images, radar remote sensing data and terrain elevation model data of the target area and constructing a target remote sensing data set, and performing surface coverage feature enhancement processing on the target remote sensing data set to generate an enhanced multi-source feature set including karst area surface coverage type features, terrain undulation features and radar scattering characteristic features, not only the detail expression ability of surface information is improved, but also the distinguishability between different geological bodies is enhanced through multi-feature fusion. Further, a pre-trained linear structure recognition model is called to perform linear structure recognition on the enhanced multi-source feature set, and an initial linear structure recognition result including linear structure spatial distribution data and structure type labels is generated. By using the powerful feature learning and classification ability of the deep learning model, the automatic and intelligent recognition of linear structures is realized, greatly improving the recognition efficiency and accuracy. Particularly importantly, based on the geological background data of the target area, spatial topological relationship verification processing is performed on the initial linear structure recognition result to generate an optimized linear structure recognition result. By introducing geological prior knowledge to constrain and correct the recognition result, the phenomena of misrecognition and missed recognition are effectively eliminated, and the reliability and geological rationality of the recognition result are significantly improved. Thus, not only the accuracy and efficiency of linear structure recognition are improved, but also the geological interpretability and practicality of the recognition result are enhanced. The optimized linear structure recognition result can accurately output the spatial positions and structure confidence parameters of faults and fractures in the karst area. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 is a schematic execution flowchart of a method for identifying linear structures in karst areas based on multi-source remote sensing provided by an embodiment of the present invention.

[0012] Figure 2 is a schematic diagram of exemplary hardware and software components of a linear structure recognition system for karst areas based on multi-source remote sensing provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0013] The present invention will be specifically described below with reference to the accompanying drawings of the specification. Figure 1 is a schematic flowchart of a method for identifying linear structures in karst areas based on multi-source remote sensing provided by an embodiment of the present invention. The method for identifying linear structures in karst areas based on multi-source remote sensing will be introduced in detail below.

[0014] Step S110: Obtain the optical remote sensing image, radar remote sensing data, and terrain elevation model data of the target area, and construct a target remote sensing data set.

[0015] In this embodiment, a region with typical karst geomorphic features is used as an example for illustration. In this target area, an optical remote sensing device carried by a remote sensing satellite can be used to observe the area at a specific time point, so as to obtain the optical remote sensing image covering the area. For example, the resolution of the optical remote sensing image is set to 10 meters, and the image data is presented as a multi-band image matrix, covering multiple band information such as red, green, blue, and near-infrared. In the red band, the numerical range is between 0 and 255, and different numerical values represent the reflection degree of different objects on the surface of the target area to the red-band light.

[0016] At the same time, a radar remote sensing device can be used to scan the target area at a 45-degree incident angle using the C band (wavelength of about 5 cm) to obtain radar remote sensing data. The above radar remote sensing data is stored in the form of complex numbers and contains the backscattering information on the surface of the target area. For example, in the radar scattering intensity data of a certain area, the numerical values fluctuate between -20 dB and 0 dB, and different numerical values reflect the different scattering characteristics of the area.

[0017] In addition, the terrain elevation model data of the target area can be obtained through means such as aerial photogrammetry or ground survey. The terrain elevation model data is presented in the form of a regular grid, and each grid point represents the altitude at a specific location. For example, the altitude range of this area is between 100 meters and 500 meters, and different altitude values correspond to different numerical values in the grid data.

[0018] Then, the obtained optical remote sensing image, radar remote sensing data, and terrain elevation model data are integrated together to construct a target remote sensing data set.

[0019] Step S120: Perform surface cover feature enhancement processing on the target remote sensing data set to generate an enhanced multi-source feature set. The enhanced multi-source feature set includes karst area surface cover type features, terrain undulation features, and radar scattering characteristic features.

[0020] In this embodiment, for optical remote sensing images, radiation correction processing is first carried out. With the help of specific radiation correction algorithms, based on the sensor parameters of remote sensing satellites and atmospheric correction models, the optical remote sensing images are processed to obtain standardized optical images after radiation correction. For example, using the MODTRAN atmospheric radiative transfer model and combining the atmospheric parameters (such as atmospheric aerosol concentration, water vapor content, etc.) when the remote sensing image is acquired, radiation correction calculations are performed for each pixel point in the optical remote sensing image. After that, the characteristics of the surface cover types in the karst area are extracted from the standardized optical image. Taking the vegetation coverage index as an example, the NDVI (Normalized Difference Vegetation Index) algorithm is used and calculated through the values of the near-infrared band and the red band. Assuming the value of the near-infrared band is NIR and the value of the red band is R, then NDVI = (NIR - R) / (NIR + R). In this target area, the calculated NDVI values are between -0.2 and 0.8, and different values represent different vegetation coverage degrees.

[0021] For the rock exposure rate parameter, it can be determined by identifying and counting the rock areas in the image. For example, using the spectral characteristics of rocks in specific bands, the pixel points in the image are classified into rock and non-rock areas, and the proportion of the pixel points in the rock area to the total pixel points is calculated to obtain the rock exposure rate parameter. In this area, the rock exposure rate is between 10% and 30%. For the water body distribution density parameter, it is determined by identifying the water body areas in the image and calculating their area ratios. For example, using the low reflectivity characteristics of water bodies in the near-infrared band, the water body areas in the image are segmented, and the proportion of the water body area to the area of the entire target area is calculated. The water body distribution density in this area is between 5% and 15%.

[0022] Then, slope calculation and processing are performed on the terrain elevation model data. Specifically, a slope calculation algorithm based on the Digital Elevation Model (DEM) can be adopted. With each grid point as the center, the elevation change rate between it and the surrounding adjacent grid points is calculated to obtain the slope value of this point. For example, for a certain grid point, the elevations of its eight surrounding adjacent grid points are h1, h2, ..., h8, and the slope value of this point is calculated through a specific formula (such as the slope calculation formula for a 3×3 window). In the entire target area, the slope value ranges from 0 degrees to 60 degrees, and thus the terrain undulation characteristics are generated. Polarization decomposition processing is performed on the radar remote sensing data to extract the radar scattering characteristic features. For example, a specific polarization decomposition algorithm (such as the Cloude-Pottier decomposition algorithm) is used to decompose the polarization information in the radar remote sensing data to obtain the surface scattering intensity parameter and the secondary scattering phase difference parameter. In this target area, the surface scattering intensity parameter ranges from -15 dB to 5 dB, and the secondary scattering phase difference parameter ranges from 0 degrees to 180 degrees. Finally, the surface cover type characteristics, terrain undulation characteristics, and radar scattering characteristic features in the karst area are subjected to spatial registration and scale unification processing to generate an enhanced multi-source feature set. For example, different feature data are spatially registered through a geographic coordinate system to make them correspond one by one in terms of spatial position; methods such as resampling are used to unify the scales of feature data with different resolutions to ensure the consistency and compatibility of the data.

[0023] Step S130, call the pre-trained linear structure recognition model to perform linear structure recognition on the enhanced multi-source feature set, and generate an initial linear structure recognition result. The initial linear structure recognition result includes the spatial distribution data of the linear structure and the structure type label.

[0024] In this embodiment, the enhanced multi-source feature set can be input into the feature encoding layer of the pre-trained linear structure recognition model to respectively generate an optical feature encoding vector, a terrain feature encoding vector, and a radar feature encoding vector. For example, the optical feature encoding vector is a vector with a dimension of 512, and its value is obtained by the feature encoding layer performing a non-linear transformation on the surface cover type characteristics in the karst area; the terrain feature encoding vector has a dimension of 256 and is generated by encoding the terrain undulation characteristics; the radar feature encoding vector has a dimension of 384 and is the result of encoding the radar scattering characteristic features.

[0025] Next, a first attention weight matrix is constructed based on the spatial correlation between the optical feature encoding vector and the terrain feature encoding vector. For example, by calculating the cosine similarity between the two vectors in each dimension, a 512×256 similarity matrix is obtained, and then through a series of non-linear transformations (such as the Softmax function), the first attention weight matrix is generated. A second attention weight matrix is constructed based on the scattering characteristic correlation between the radar feature encoding vector and the terrain feature encoding vector. For example, by analyzing the elements related to the scattering characteristics in the two vectors and using a specific algorithm (such as an algorithm based on distance metric) to calculate the correlation, a 384×256 matrix is obtained, and after processing, the second attention weight matrix is generated.

[0026] Then, the optical feature encoding vector and the terrain feature encoding vector are weighted and fused according to the first attention weight matrix to generate a first fused feature. For example, each element of the optical feature encoding vector is multiplied by the corresponding row element in the first attention weight matrix, and then summed by column to obtain a first fused feature vector with a dimension of 256. The radar feature encoding vector and the terrain feature encoding vector are weighted and fused according to the second attention weight matrix to generate a second fused feature. For example, in the same way, the radar feature encoding vector is operated with the second attention weight matrix to obtain a second fused feature vector with a dimension of 256.

[0027] Finally, the first fused feature and the second fused feature are subjected to cross-validation processing to generate an initial linear structure recognition result containing the spatial distribution data of the linear structure and the structure type label. For example, through a specific cross-validation algorithm (such as concatenating the two fused feature vectors and inputting them into a fully connected layer for classification and localization prediction), the spatial distribution data of the linear structure is obtained, which is represented in the form of the coordinate points of a series of line segments; at the same time, the structure type labels are obtained, such as different type labels like "fault" and "fracture".

[0028] Step S140, perform a spatial topological relationship verification process on the initial linear structure recognition result based on the geological background data of the target area to generate an optimized linear structure recognition result. The optimized linear structure recognition result is used to output the spatial positions and structure confidence parameters of the faults and fractures in the karst area.

[0029] In this embodiment, first, the known fault zone distribution data and the rock formation attitude parameters in the geological background data are extracted to construct a geological prior knowledge base. For example, the known fault zone distribution data is stored in the form of vector graphics, including the position coordinate information of the fault zone; the rock formation attitude parameters include information such as the dip angle and the dip, the dip angle range is between 0 degrees and 360 degrees, and the dip range is between 0 degrees and 90 degrees.

[0030] Next, perform a spatial overlay analysis on the initial linear structure recognition results and the geological prior knowledge base to determine the spatial consistency parameter and the formation attitude matching degree parameter between the spatial distribution data of the linear structure and the known fault zones. Conduct buffer analysis on the known fault zone distribution data in the geological prior knowledge base to generate buffer range data for each known fault zone. For example, with each fault zone as the center, set a buffer with a radius of 50 meters to obtain the buffer range data. Perform spatial overlay processing on the spatial distribution data of the linear structures in the initial linear structure recognition results and the buffer range data, extract the linear structure line segments that intersect with the buffer range data, calculate the ratio of the length of the intersecting line segments to the total length of the corresponding linear structure line segments, and generate the spatial consistency parameter. For example, if the total length of a linear structure line segment is 200 meters and the length of the line segment intersecting with the buffer range is 100 meters, then the spatial consistency parameter of this linear structure is 0.5.

[0031] Extract the strike angle data of each linear structure line segment from the initial linear structure recognition results to form a set of linear structure strike angles. For example, by calculating the coordinate changes of the endpoints of the linear structure line segment, obtain its strike angle, which ranges from 0 degrees to 360 degrees. Calculate the angular difference between the set of linear structure strike angles and the dip angle data of the formation attitude parameters in the geological prior knowledge base to generate the corresponding strike-dip angle difference for each linear structure line segment. For example, if the strike angle of a linear structure is 45 degrees and the dip angle of a formation is 60 degrees, then the strike-dip angle difference is 15 degrees. Perform normalization processing based on the absolute value of the strike-dip angle difference to generate the formation attitude matching degree parameter. For example, divide the absolute value of the strike-dip angle difference by 180 degrees (the maximum possible difference) to obtain the formation attitude matching degree parameter, and the smaller the absolute value, the higher the formation attitude matching degree parameter.

[0032] Then, adjust the structural confidence parameter in the initial linear structure recognition results according to the spatial consistency parameter and the formation attitude matching degree parameter to generate an optimized linear structure recognition result. Normalize the spatial consistency parameter and the formation attitude matching degree parameter to obtain the normalized spatial consistency parameter and the normalized formation attitude matching degree parameter. For example, normalize the spatial consistency parameter by dividing it by 1 (the maximum value), and perform corresponding normalization processing on the formation attitude matching degree parameter. Perform weighted summation on the normalized spatial consistency parameter and the normalized formation attitude matching degree parameter to generate a comprehensive adjustment coefficient. For example, set the weight of the spatial consistency parameter to 0.6 and the weight of the formation attitude matching degree parameter to 0.4, and perform weighted summation to obtain the comprehensive adjustment coefficient.

[0033] The comprehensive adjustment coefficient is multiplied by the structural confidence parameter in the initial linear structure recognition result to generate the adjusted structural confidence parameter. For example, if the initial structural confidence parameter is 0.8 and the comprehensive adjustment coefficient is 0.9, the adjusted structural confidence parameter is 0.72. The adjusted structural confidence parameter is truncated, and the parameter value exceeding the preset confidence threshold (such as 0.9) is set as the preset confidence threshold to generate the optimized structural confidence parameter. For example, if the adjusted structural confidence parameter is 0.72 and does not exceed the threshold, the optimized structural confidence parameter is still 0.72. Finally, the optimized structural confidence parameter is merged with the linear structure spatial distribution data and structural type label in the initial linear structure recognition result to generate the optimized linear structure recognition result, which is used to output the spatial position and structural confidence parameters of faults and fissures in the karst area.

[0034] Based on the above steps, the embodiment of the present application obtains the optical remote sensing image, radar remote sensing data and terrain elevation model data of the target area and constructs the target remote sensing data set, performs surface coverage feature enhancement processing on the target remote sensing data set, generates an enhanced multi-source feature set including the surface coverage type feature of the karst area, the terrain relief feature and the radar scattering characteristic feature, not only improves the detail expression ability of the surface information, but also enhances the distinction between different geological bodies through multi-feature fusion, further, calls the pre-trained linear structure recognition model to recognize the enhanced multi-source feature set for linear structure, generates an initial linear structure recognition result including linear structure spatial distribution data and structure type label, and utilizes the powerful feature learning and classification capabilities of the deep learning model to realize the automatic and intelligent recognition of linear structure, greatly improving the recognition efficiency and accuracy. More importantly, the spatial topological relationship verification processing is performed on the initial linear structure recognition result based on the geological background data of the target area to generate an optimized linear structure recognition result, and the recognition result is constrained and corrected by introducing geological prior knowledge, effectively eliminating the phenomenon of misidentification and missed recognition, and significantly improving the reliability and geological rationality of the recognition result. Therefore, not only the accuracy and efficiency of linear structure identification are improved, but also the geological interpretation and practicality of the identification results are enhanced. The optimized linear structure identification results can accurately output the spatial location and structural confidence parameters of faults and cracks in the karst area.

[0035] In a possible implementation, step S120 includes:

[0036] Step S121, performing radiation correction processing on the optical remote sensing image to obtain a standardized optical image after radiation correction, and extracting the surface cover type characteristics of the karst area of the standardized optical image. The surface cover type characteristics of the karst area include vegetation coverage index, rock layer exposure rate parameter and water body distribution density parameter.

[0037] In this embodiment, when acquiring optical remote sensing images, due to the influence of factors such as atmospheric scattering and absorption, there will be radiation errors in the images. In order to obtain a standardized optical image after radiometric correction, a correction method based on the radiative transfer model can be used. Taking the MODTRAN radiative transfer model as an example, it takes into account various factors affecting radiative transfer, such as atmospheric components (such as oxygen, carbon dioxide, water vapor, etc.), aerosol characteristics (such as concentration, particle size distribution), solar altitude angle, observation angle, etc. When processing the optical remote sensing image of the target area, according to the satellite observation parameters recorded when acquiring the image, including satellite altitude, sensor angle, etc., and the meteorological data of this area (such as the atmospheric aerosol concentration is 0.2mg / m 3 , and the water vapor content is 5g / m 3 , etc.), they are input into the MODTRAN model. Then, each pixel in the image is calculated. By simulating the radiative transfer process in the atmosphere, the interference of the atmosphere on radiation is removed, and thus a standardized optical image after radiometric correction is obtained.

[0038] For the vegetation coverage index, the NDVI (Normalized Difference Vegetation Index) algorithm is used for calculation. In the standardized optical image of the target area, the near-infrared band data presents a series of values. Suppose the near-infrared band value of a certain pixel is NIR = 0.6, and the red band value is R = 0.2. According to the NDVI algorithm, first calculate the numerator part, that is, NIR - R = 0.6 - 0.2 = 0.4; then calculate the denominator part, NIR + R = 0.6 + 0.2 = 0.8; finally, the NDVI value of this pixel is 0.4 ÷ 0.8 = 0.5. The same calculation is performed for all pixels in the entire image to obtain the vegetation coverage index distribution of the entire target area.

[0039] For the rock exposure rate parameter, the spectral characteristics of rocks in specific bands are used to identify rock areas. Through a large number of field investigations and spectral measurements in the target area, it is determined that the reflectance characteristics of rocks in the visible light band are significantly different from those of other ground objects. For example, in the blue band, the reflectance of rocks is relatively low, generally between 0.1 and 0.3. By setting an appropriate threshold, threshold segmentation is performed on the blue band data of the standardized optical image. Suppose the threshold is set to 0.3, and the pixel points in the blue band with values less than 0.3 are marked as pixel points in the rock area. The number of pixel points in the rock area is counted as 10,000, and the total number of pixel points in the entire image is 100,000, then the rock exposure rate parameter is 10,000 ÷ 100,000 × 100% = 10%.

[0040] For the water body distribution density parameter, utilize the low reflectivity characteristic of the water body in the near-infrared band. In the near-infrared band data of the standardized optical image, set a relatively low threshold, such as 0.1. Identify the pixel points with near-infrared band values less than 0.1 as the pixel points in the water body area. The number of pixel points in the water body area is counted as 5000. The area of the entire target area image is calculated through the number of pixel points and the actual area represented by each pixel point (assuming that each pixel point represents 100 square meters and the total number of pixel points in the image is 100000, then the total area is 100×100000 = 10000000 square meters). The area of the water body area is 5000×100 = 500000 square meters, so the water body distribution density parameter is 500000÷10000000×100% = 5%.

[0041] Step S122, perform slope calculation processing on the terrain elevation model data to generate terrain undulation characteristics, and perform polarization decomposition processing on the radar remote sensing data to extract radar scattering characteristic features. The radar scattering characteristic features include surface scattering intensity parameters and secondary scattering phase difference parameters.

[0042] Taking a certain grid point as an example, the elevation heights of the eight adjacent grid points around this grid point are h1 = 150 meters, h2 = 148 meters, h3 = 152 meters, h4 = 145 meters, h5 = 155 meters, h6 = 147 meters, h7 = 153 meters, and h8 = 149 meters respectively. Using the slope calculation formula based on a 3×3 window, first calculate the slope component in the east-west direction, that is, (h3 + 2×h5 + h7 - h1 - 2×h4 - h6)÷(8×dx). Assuming dx is the grid spacing, here dx = 10 meters, then the slope component in the east-west direction is (152 + 2×155 + 153 - 150 - 2×145 - 147)÷(8×10) = (152 + 310 + 153 - 150 - 290 - 147)÷80 = 28÷80 = 0.35. Then calculate the slope component in the north-south direction, that is, (h7 + 2×h8 + h6 - h1 + 2×h2 + h4)÷(8×dy). Assuming dy = 10 meters, then the slope component in the north-south direction is (153 + 2×149 + 147 - 150 + 2×148 + 145)÷(8×10) = (153 + 298 + 147 - 150 + 296 + 145)÷80 = 799÷80 = 9.9875. Finally, calculate the slope value of this grid point according to the Pythagorean theorem, that is, the slope value = √(0.35² + 9.9875²)≈9.993. Perform such calculations on all grid points of the entire terrain elevation model data to generate terrain undulation characteristics.

[0043] On this basis, the Cloude-Pottier decomposition algorithm can be adopted. The Cloude-Pottier decomposition algorithm is based on the polarization characteristics of radar echo signals. For the radar remote sensing data of the target area, the signal intensities and phase information of different polarization channels (such as HH, HV, VH, VV) are first analyzed. Suppose in a certain area, the signal intensity of the HH channel is \(I_{HH} = 10 dB\), the signal intensity of the HV channel is \(I_{HV} = 5 dB\), the signal intensity of the VH channel is \(I_{VH} = 3 dB\), and the signal intensity of the VV channel is \(I_{VV} = 8 dB\). Through a series of calculation steps in the Cloude-Pottier decomposition algorithm, first calculate the covariance matrix, and then perform operations such as eigenvalue decomposition. For example, calculate the element \(C_{11}=I_{HH}\), \(C_{12}=\sqrt{I_{HH}\times I_{HV}}\times exp(j\times\varphi_{HV})\) (where \(\varphi_{HV}\) is the phase difference between the HV channel signal and the HH channel signal), \(C_{22}=I_{VV}\), etc. After a series of matrix operations and eigenvalue decompositions, the surface scattering intensity parameter and the secondary scattering phase difference parameter are obtained. Suppose the surface scattering intensity parameter of this area is finally calculated to be -10 dB, and the secondary scattering phase difference parameter is 90 degrees.

[0044] Step S123: Perform spatial registration and scale unification processing on the karst area surface coverage type feature, the terrain undulation degree feature, and the radar scattering characteristic feature to generate the enhanced multi-source feature set.

[0045] In this embodiment, the geographic information system (GIS) technology can be used. According to the geographic coordinate information corresponding to the karst area surface coverage type feature, the terrain undulation degree feature, and the radar scattering characteristic feature, they are accurately aligned in space to ensure that the karst area surface coverage type feature, the terrain undulation degree feature, and the radar scattering characteristic feature accurately correspond geographically. For scale unification, since the resolutions of different feature data may be different, a resampling method is adopted. For example, the resolution of the karst area surface coverage type feature data is 10 meters, the resolution of the terrain undulation degree feature data is 20 meters, and the resolution of the radar scattering characteristic feature data is 15 meters. Resample the terrain undulation degree feature data and the radar scattering characteristic feature data to the same 10-meter resolution as the karst area surface coverage type feature data. Process the data through resampling algorithms such as bilinear interpolation, and finally generate the enhanced multi-source feature set.

[0046] In a possible implementation manner, step S130 includes:

[0047] Step S131: Input the enhanced multi-source feature set into the feature encoding layer of the pre-trained linear structure recognition model to generate an optical feature encoding vector, a terrain feature encoding vector, and a radar feature encoding vector respectively.

[0048] In this embodiment, the enhanced multi-source feature set includes the karst area surface cover type feature, the terrain undulation degree feature, and the radar scattering characteristic feature. Taking the karst area surface cover type feature as an example, the vegetation coverage index is calculated by the NDVI algorithm and shows different numerical distributions in the target area. For example, in a certain area, after calculating the NDVI values of multiple pixels, the average vegetation coverage index of this area is 0.6; the rock exposure rate parameter is obtained by identifying and counting the rock areas in the standardized optical image, and the rock exposure rate in this area is calculated to be 15%; the water body distribution density parameter is obtained by identifying the water body areas in the image and calculating the area ratio, and the water body distribution density in this area is 8%. The terrain undulation degree feature is generated by calculating the slope of the terrain elevation model data. For example, the slope value of a certain area is calculated to be between 10 degrees and 30 degrees. The radar scattering characteristic feature includes the surface scattering intensity parameter and the secondary scattering phase difference parameter. In this area, the surface scattering intensity parameter is calculated to be -12 dB, and the secondary scattering phase difference parameter is 100 degrees.

[0049] After inputting the above features into the feature encoding layer of the pre-trained linear structure recognition model, an optical feature encoding vector, a terrain feature encoding vector, and a radar feature encoding vector are generated respectively. The optical feature encoding vector is obtained by performing a non-linear transformation on the karst area surface cover type feature. Assuming the dimension of this vector is 512, the value of each of its elements is obtained by performing matrix multiplication on the specific weight matrix in the feature encoding layer and the input karst area surface cover type feature, and then processed by an activation function (such as the ReLU function). For example, for the vegetation coverage index 0.6 in the karst area surface cover type feature, after weighted summation of the corresponding elements in the weight matrix, and then processed by the ReLU function (if the weighted sum is negative, output 0, if it is positive, output the positive number), the value of an element at a certain position in the optical feature encoding vector is obtained. The dimension of the terrain feature encoding vector is 256, which is also generated by encoding the terrain undulation degree feature. Through similar matrix operations and activation function processing, information such as the slope value in the terrain undulation degree feature is converted into vector element values. The dimension of the radar feature encoding vector is 384, which is used to encode the radar scattering characteristic feature. For example, information such as the surface scattering intensity parameter -12 dB and the secondary scattering phase difference parameter 100 degrees is converted into vector elements through specific operations.

[0050] Step S132: Construct a first attention weight matrix based on the spatial correlation between the optical feature encoding vector and the terrain feature encoding vector, and construct a second attention weight matrix based on the scattering characteristic correlation between the radar feature encoding vector and the terrain feature encoding vector.

[0051] In this embodiment, the correlation between two vectors can be measured by calculating the cosine similarity in each dimension. For example, for the first-dimensional element value a1 of the optical feature encoding vector and the first-dimensional element value b1 of the terrain feature encoding vector, the calculation process of the cosine similarity is as follows: First, calculate the dot product of vector a (optical feature encoding vector) and vector b (terrain feature encoding vector), that is, a1*b1 + a2*b2 +... + a512*b256 (where a2, a3, etc. are the other-dimensional element values of the optical feature encoding vector, and b2, b3, etc. are the other-dimensional element values of the terrain feature encoding vector), and the dot product result is S1. Then calculate the norm of vector a, that is, √(a1² + a2² +... + a512²), and the norm is |a|; calculate the norm of vector b, that is, √(b1² + b2² +... + b256²), and the norm is |b|. Finally, the cosine similarity is S1 / (|a|*|b|). Such calculations are performed for all dimensions to obtain a 512×256 similarity matrix. After being processed by the Softmax function, the role of the Softmax function is to convert each element in the similarity matrix into a probability value, so that the sum of the elements in each row is 1. For example, for a row of elements [x1, x2, x3] in the similarity matrix, after being processed by the Softmax function, the first element becomes exp(x1) / (exp(x1) + exp(x2) + exp(x3)), the second element becomes exp(x2) / (exp(x1) + exp(x2) + exp(x3)), and the third element becomes exp(x3) / (exp(x1) + exp(x2) + exp(x3)), and finally the first attention weight matrix is generated.

[0052] Furthermore, the correlation can be calculated by analyzing the elements related to the scattering characteristics in two vectors and using an algorithm based on distance measurement. For example, for the element c1 related to the surface scattering intensity in the radar feature encoding vector and the element d1 related to the influence of the terrain on scattering in the terrain feature encoding vector, calculate the Euclidean distance between them, that is, √((c1 - d1)² + (c2 - d2)² +... + (c384 - d256)²) (where c2, c3, etc. are the other related element values of the radar feature encoding vector, and d2, d3, etc. are the other related element values of the terrain feature encoding vector), and the distance value D1 is obtained. Such calculations are performed for all related elements to obtain a 384×256 distance matrix. Then, through a specific conversion function (such as converting the distance value into a similarity value, for example, similarity = 1 / (1 + D1)), the distance matrix is converted into a similarity matrix. After being processed by a function similar to the Softmax function, the second attention weight matrix is generated.

[0053] Step S133, perform weighted fusion on the optical feature encoding vector and the terrain feature encoding vector according to the first attention weight matrix to generate a first fusion feature, and perform weighted fusion on the radar feature encoding vector and the terrain feature encoding vector according to the second attention weight matrix to generate a second fusion feature.

[0054] For example, for the first element a1 of the optical feature encoding vector, the first element b1 of the terrain feature encoding vector, and the weight w11 at the corresponding position in the first attention weight matrix, the calculation process of weighted fusion is as follows: First, calculate the weighted value of the element of the optical feature encoding vector, that is, a1 * w11, and the weighted value of the element of the terrain feature encoding vector, that is, b1 * w11. Then add these two weighted values to obtain the value of the first element after fusion. Such calculations are performed for all dimensions, and finally a first fusion feature vector with a dimension of 256 is generated. Perform weighted fusion on the radar feature encoding vector and the terrain feature encoding vector according to the second attention weight matrix to generate a second fusion feature. In the same way, for the first element c1 of the radar feature encoding vector, the first element b1 of the terrain feature encoding vector, and the weight w21 at the corresponding position in the second attention weight matrix, calculate the weighted value of the element of the radar feature encoding vector c1 * w21, and the weighted value of the element of the terrain feature encoding vector b1 * w21. Add the two weighted values to obtain the value of the first element after fusion, and so on, to generate a second fusion feature vector with a dimension of 256.

[0055] Step S134, perform cross - validation processing on the first fusion feature and the second fusion feature to generate the initial linear structure recognition result including linear structure spatial distribution data and structure type labels.

[0056] For example, the two fusion feature vectors can be concatenated and then input into a fully - connected layer for classification and localization prediction. For example, the first fusion feature vector and the second fusion feature vector are concatenated in sequence into a vector with a dimension of 512 and input into the fully - connected layer. The weight matrix in the fully - connected layer performs matrix multiplication with the concatenated vector. Assume the weight matrix is W with a dimension of 512×1024, and the concatenated vector is V with a dimension of 512×1, then the operation result is a vector with a dimension of 1024×1. Then it is processed by an activation function (such as the Sigmoid function). The Sigmoid function converts each element in the vector into a probability value between 0 and 1. Classification judgment is made according to the above probability values. For example, if the probability value at a certain position is greater than 0.5, it is judged as a certain structure type (such as "fault"), and if it is less than 0.5, it is judged as other types (such as "fracture"). At the same time, through the analysis of some specific elements in the vector, the spatial distribution data of the linear structure is obtained, which is represented in the form of coordinate points of a series of line segments, so as to generate the initial linear structure recognition result.

[0057] In a possible implementation, step S140 includes:

[0058] Step S141, extract the known fault zone distribution data and rock formation attitude parameters in the geological background data, and construct a geological prior knowledge base.

[0059] For example, the known fault zone distribution data is stored in the form of vector graphics, including the position coordinate information of the fault zone. For example, in this target area, there is a known fault zone, whose starting point coordinates are (100, 200) and the ending point coordinates are (300, 400). The rock formation attitude parameters include information such as dip angle and inclination angle. Suppose the dip angle of a certain rock formation is 120 degrees and the inclination angle is 30 degrees. Integrate the above data together to construct a geological prior knowledge base.

[0060] Step S142, perform a spatial overlay analysis on the initial linear structure recognition result and the geological prior knowledge base to determine the spatial consistency parameter between the linear structure spatial distribution data and the known fault zone and the rock formation attitude matching degree parameter.

[0061] Step S143, adjust the structure confidence parameter in the initial linear structure recognition result according to the spatial consistency parameter and the rock formation attitude matching degree parameter to generate the optimized linear structure recognition result. Among them, the structure confidence parameter is used to characterize the matching degree between the linear structure recognition result and the geological background data.

[0062] For example, in a possible implementation, step S142 includes:

[0063] Step S1421, perform buffer analysis processing on the known fault zone distribution data in the geological prior knowledge base to generate buffer range data for each known fault zone.

[0064] For example, in the geological prior knowledge base of the target area, there is a known fault zone, whose spatial position is represented by a series of coordinate points. The starting point coordinates are (1000, 2000), and the ending point coordinates are (3000, 4000). The above coordinates are determined based on a specific geographical coordinate system of this area. To generate buffer zone range data, the buffer zone radius is set to 100 meters. Taking the starting point (1000, 2000) as the center and 100 meters as the radius to determine the buffer zone range. For the points on the buffer zone boundary, the distance from them to the center of the circle needs to satisfy the distance formula. For example, for the point (x, y), its distance to the center of the circle (1000, 2000) is the square root of [(x - 1000)^2 + (y - 2000)^2]. When this distance is equal to 100, this point is on the inner circular boundary of the buffer zone with the starting point as the center. By calculating a series of points that meet this condition, a set of discrete points of the circular boundary of the buffer zone with the starting point as the center is obtained. Using the same method, taking the ending point (3000, 4000) as the center and 100 meters as the radius to determine another set of discrete points of the circular boundary of the buffer zone. At the same time, considering the buffer zone range on both sides of the fault zone line segment, it is determined by calculating two parallel lines that are parallel to the fault zone line segment and at a distance of 100 meters. Suppose the slope of the fault zone line segment is k = (4000 - 2000) ÷ (3000 - 1000) = 1, then the equations of the two parallel lines that are parallel to it and at a distance of 100 meters can be determined by the point-to-line distance formula. For the line Ax + By + C = 0 (here A = -1, B = 1, C = 2000 - 1000), the point-to-line distance formula for the point (x0, y0) is |Ax0 + By0 + C| ÷ the square root of (A^2 + B^2). By setting the distance to 100, the points on the two parallel lines are obtained. Finally, all the above points are integrated to generate the buffer zone range data of the known fault zone.

[0065] Step S1422: Perform spatial overlay processing on the linear structure spatial distribution data in the initial linear structure recognition result and the buffer zone range data, extract the linear structure line segments that intersect with the buffer zone range data, calculate the ratio of the length of the intersecting line segments to the total length of the corresponding linear structure line segments, and generate the spatial consistency parameter.

[0066] For example, in the initial linear structure recognition result, there is a linear structure segment with a starting coordinate of (1500, 2500) and an ending coordinate of (2500, 3500). To determine whether it intersects with the buffer range data and the length of the intersecting segment, it is necessary to judge one by one whether the points on this linear structure segment are within the buffer range. For any point (x, y) on the linear structure segment, its coordinates can be represented by the parametric equations x = 1500 + t×(2500 - 1500), y = 2500 + t×(3500 - 2500) (0 ≤ t ≤ 1). Substitute the above points into the judgment conditions of the buffer range, that is, judge the distance relationship between the above points and each boundary of the known fault zone buffer. After calculation, it is found that there is an intersecting part between this linear structure segment and the buffer range. Suppose the starting coordinate of the intersecting part is (1800, 2800) and the ending coordinate is (2200, 3200). Calculate the length of the intersecting segment. According to the distance formula between two points, the length of the intersecting segment is the square root of [(2200 - 1800) squared + (3200 - 2800) squared], that is, the square root of [400 squared + 400 squared] = the square root of (160000 + 160000) = the square root of 320000 ≈ 565.69 meters. The total length of the original linear structure segment is the square root of [(2500 - 1500) squared + (3500 - 2500) squared] = the square root of [1000 squared + 1000 squared] = the square root of (1000000 + 1000000) = the square root of 2000000 ≈ 1414.21 meters. Then the spatial consistency parameter is 565.69 ÷ 1414.21 ≈ 0.4.

[0067] Step S1423, extract the trend angle data of each linear structure segment from the initial linear structure recognition result to form a linear structure trend angle set.

[0068] For example, for the above-mentioned linear structure segment, the calculation method of its trend angle is as follows: First, calculate the change amounts in the x-axis and y-axis directions of this segment. The change amount in the x-axis direction is 2500 - 1500 = 1000, and the change amount in the y-axis direction is 3500 - 2500 = 1000. The trend angle is calculated through the arctangent function, that is, the trend angle = arctangent (y-axis change amount ÷ x-axis change amount), which is arctangent (1000 ÷ 1000) = arctangent (1) = 45 degrees. Perform such calculations on all linear structure segments in the initial linear structure recognition result, collect the obtained trend angle data, and form a linear structure trend angle set. For example, there is another linear structure segment. After calculation, its trend angle is 60 degrees. Put these two angle values together, and the linear structure trend angle set contains these two values of 45 degrees and 60 degrees.

[0069] Step S1424: Calculate the angular difference between the set of linear structure strike angles and the dip angle data of the rock formation attitude parameters in the geological prior knowledge base to generate the strike-dip angle differences corresponding to each linear structure segment.

[0070] For example, in the geological prior knowledge base, the dip angle of a certain rock formation attitude parameter is 135 degrees. For the 45-degree strike angle in the set of linear structure strike angles, its strike-dip angle difference is 135 - 45 = 90 degrees; for the 60-degree strike angle, its strike-dip angle difference is 135 - 60 = 75 degrees. In this way, the strike-dip angle differences corresponding to each linear structure segment are obtained.

[0071] Step S1425: Perform normalization processing according to the absolute value of the strike-dip angle difference to generate the rock formation attitude matching degree parameter. Among them, the smaller the absolute value, the higher the rock formation attitude matching degree parameter.

[0072] In this embodiment, the purpose of normalization is to convert the above difference into the range of 0 to 1 to better measure the matching degree. The maximum possible angular difference is 180 degrees. For the 90-degree strike-dip angle difference, its normalized rock formation attitude matching degree parameter is 90÷180 = 0.5; for the 75-degree strike-dip angle difference, the normalized rock formation attitude matching degree parameter is 75÷180≈0.42. It can be seen here that the smaller the absolute value, the higher the rock formation attitude matching degree parameter, indicating that the matching degree between the linear structure and the rock formation attitude is better.

[0073] For example, in a possible implementation manner, step S143 includes:

[0074] Step S1431: Perform normalization processing on the spatial consistency parameter and the rock formation attitude matching degree parameter to obtain the normalized spatial consistency parameter and the normalized rock formation attitude matching degree parameter.

[0075] In this embodiment, first perform normalization processing on the spatial consistency parameter and the rock formation attitude matching degree parameter. Since the spatial consistency parameter 0.4 and the rock formation attitude matching degree parameter 0.42 are already within the range of 0 to 1, the normalized spatial consistency parameter is still 0.4, and the normalized rock formation attitude matching degree parameter is still 0.42.

[0076] Step S1432: Perform weighted summation on the normalized spatial consistency parameter and the normalized rock formation attitude matching degree parameter to generate a comprehensive adjustment coefficient.

[0077] For example, assume that the weight of the spatial consistency parameter is set to 0.6, and the weight of the formation attitude matching degree parameter is 0.4. Then the comprehensive adjustment coefficient is 0.4×0.6 + 0.42×0.4. First, calculate the multiplication parts: 0.4×0.6 = 0.24, 0.42×0.4 = 0.168. Then add the two parts together: 0.24 + 0.168 = 0.408. The obtained comprehensive adjustment coefficient is 0.408.

[0078] Step S1433: Multiply the comprehensive adjustment coefficient by the structural confidence parameter in the initial linear structure recognition result to generate an adjusted structural confidence parameter.

[0079] For example, assume that the structural confidence parameter in the initial linear structure recognition result is 0.8. Then the adjusted structural confidence parameter is 0.408×0.8. The calculation process is 0.408×0.8 = (408×8)÷1000 = 3264÷1000 = 0.3264.

[0080] Step S1434: Perform a truncation process on the adjusted structural confidence parameter, set the parameter value exceeding the preset confidence threshold to the preset confidence threshold, and generate an optimized structural confidence parameter.

[0081] For example, assume that the preset confidence threshold is 0.7. Since 0.3264 is less than 0.7, the optimized structural confidence parameter remains 0.3264.

[0082] Step S1435: Merge the optimized structural confidence parameter with the linear structure spatial distribution data and the structural type label in the initial linear structure recognition result to generate the optimized linear structure recognition result.

[0083] For example, the linear structure spatial distribution data in the initial linear structure recognition result is represented in the form of coordinate points as the coordinates of the two linear structure segments mentioned above, and the structural type labels are "fault" and "fracture". Integrate the optimized structural confidence parameter 0.3264 with the above linear structure spatial distribution data and the structural type label to form the final optimized linear structure recognition result. This result can more accurately reflect the spatial positions and structural confidence parameters of faults and fractures in the karst area, providing more reliable data support for subsequent geological analysis and research.

[0084] In a possible implementation manner, the training method of the pre-trained linear structure recognition model includes:

[0085] Step S210: Obtain a multi-source remote sensing training data set of the sample karst area and the corresponding linear structure annotation data set. The multi-source remote sensing training data set includes sample optical remote sensing images, sample radar remote sensing data, and sample terrain elevation model data.

[0086] In this embodiment, within the sample karst area, a multi-source remote sensing training data set can be obtained through remote sensing equipment. Among them, the sample optical remote sensing images are taken by high-resolution optical satellites, covering the entire sample area, with a resolution set at 5 meters. The data is presented as an image matrix containing multiple bands, such as including four bands: blue, green, red, and near-infrared. In the blue band, different ground objects reflect different intensities of blue light, with the numerical range between 0 - 255. For example, the pixel values of the blue band in a certain area mostly concentrate between 50 - 100, representing the reflection characteristics of specific ground objects in this area to blue light. The sample radar remote sensing data is obtained by scanning the sample area using a radar device with a specific frequency band (such as the X band, with a wavelength of about 3 cm). The above data is stored in complex form, recording the backscattering information of the target area surface. For example, the radar scattering intensity in a certain area fluctuates between -15 dB and -5 dB. The sample terrain elevation model data is obtained through aerial photogrammetry combined with ground surveys, presented in the form of a regular grid. Each grid point represents the altitude of a specific location, and the altitude range of this sample area is between 80 meters and 450 meters.

[0087] The corresponding linear structure annotation data set can be obtained by professional geologists through on-site investigation and analysis to mark the linear structures in the sample area. For example, the above linear structures include faults, fractures, etc. The annotation data details the spatial position coordinates and structure types of each linear structure. For example, the starting coordinates of a certain fault are (500, 800), the ending coordinates are (1200, 1500), and the structure type is annotated as "fault".

[0088] Step S220: Perform surface cover feature enhancement processing on the multi-source remote sensing training data set to generate a training multi-source feature set. The training multi-source feature set includes karst area surface cover type features, terrain undulation features, and radar scattering characteristic features.

[0089] In a possible implementation manner, step S220 includes:

[0090] Step S221: Perform radiometric correction processing on the sample optical remote sensing images to obtain standardized sample optical images, and extract the karst area surface cover type features of the standardized sample optical images.

[0091] In this embodiment, a processing flow of setting a specific radiation correction algorithm can be adopted, and the radiation error of the image can be corrected by combining the parameters of the remote sensing device and the atmospheric condition data of the sample area. For example, according to the meteorological report of the sample area, the atmospheric aerosol concentration at that time is obtained as 0.15 mg / m 3 , the water vapor content is 4 g / m 3 and other information, and the above parameters and the sensor parameters of the image itself (such as gain, offset, etc.) are input into the radiation correction algorithm. The algorithm calculates each pixel point in the image, and through a series of mathematical operations (such as the model established according to the atmospheric radiation transfer principle), removes the interference of the atmosphere on the radiation, and obtains a standardized sample optical image.

[0092] Then, the characteristics of the surface cover type in the karst area are extracted from the standardized sample optical image. Taking the vegetation coverage index as an example, the NDVI (Normalized Difference Vegetation Index) algorithm is used for calculation. For a certain pixel point in the image, the near-infrared band value is obtained as 0.7, and the red band value is 0.3. First, calculate the numerator part, that is, the near-infrared band value minus the red band value: 0.7 - 0.3 = 0.4; then calculate the denominator part, the near-infrared band value plus the red band value: 0.7 + 0.3 = 1. Then the NDVI value of this pixel point is 0.4 ÷ 1 = 0.4. Repeat this calculation process for all pixel points in the entire image to obtain the distribution of the vegetation coverage index in the entire sample area. For the rock exposure rate parameter, it is determined by identifying and counting the rock areas in the image. Using the spectral characteristics of rocks in a specific band (such as the red band), a suitable threshold (such as 0.4) is set. Pixel points with values greater than the threshold in the red band are initially judged as non-rock areas, and pixel points less than the threshold are marked as rock area pixel points. The number of rock area pixel points is counted as 8000, and the total number of pixel points in the entire image is 80000, then the rock exposure rate parameter is 8000 ÷ 80000 × 100% = 10%. For the water body distribution density parameter, using the low reflectivity characteristic of the water body in the near-infrared band, a lower threshold (such as 0.2) is set. Pixel points with near-infrared band values less than the threshold are identified as water body area pixel points. The number of water body area pixel points is counted as 6000, and the area of the entire sample area image is calculated by the number of pixel points and the actual area represented by each pixel point (assuming that each pixel point represents 25 square meters and the total number of pixel points in the image is 80000, then the total area is 25 × 80000 = 2000000 square meters). The water body area is 6000 × 25 = 150000 square meters, then the water body distribution density parameter is 150000 ÷ 2000000 × 100% = 7.5%.

[0093] Step S222, perform slope calculation processing on the sample terrain elevation model data to generate the terrain undulation characteristics.

[0094] Taking a certain grid point as an example, the elevation heights of the eight adjacent grid points around this grid point are h1 = 120 meters, h2 = 118 meters, h3 = 125 meters, h4 = 115 meters, h5 = 130 meters, h6 = 117 meters, h7 = 128 meters, and h8 = 119 meters respectively. Using the slope calculation method based on a 3×3 window, first calculate the slope component in the east-west direction. The calculation of the slope component in the east-west direction is carried out by calculating the difference in elevation heights of adjacent grid points in the east-west direction. That is, (h3 + 2×h5 + h7 - h1 - 2×h4 - h6)÷(8×dx), where dx is the grid spacing, and assume dx = 5 meters. First calculate the value inside the parentheses: h3 + 2×h5 + h7 - h1 - 2×h4 - h6 = 125 + 2×130 + 128 - 120 - 2×115 - 117 = 125 + 260 + 128 - 120 - 230 - 117 = 48. Then the slope component in the east-west direction is 48÷(8×5) = 48÷40 = 1.2. Then calculate the slope component in the north-south direction, that is, (h7 + 2×h8 + h6 - h1 + 2×h2 + h4)÷(8×dy), and assume dy = 5 meters. First calculate the value inside the parentheses: h7 + 2×h8 + h6 - h1 + 2×h2 + h4 = 128 + 2×119 + 117 - 120 + 2×118 + 115 = 128 + 238 + 117 - 120 + 236 + 115 = 714. Then the slope component in the north-south direction is 714÷(8×5) = 714÷40 = 17.85. Finally, calculate the slope value of this grid point according to the Pythagorean theorem, that is, slope value = √(1.2² + 17.85²). First calculate the squared values: 1.2² = 1.44, 17.85² = 318.6225. Add the two: 1.44 + 318.6225 = 318.6225 + 1.44 = 319.0625. Then take the square root to get the slope value approximately 17.86. Such calculations are performed on all grid points of the entire sample topographic elevation model data to generate the terrain undulation characteristics.

[0095] Step S223, perform polarization decomposition processing on the sample radar remote sensing data, and extract the radar scattering characteristic features.

[0096] In this embodiment, a specific polarization decomposition algorithm (such as the Freeman-Durden decomposition algorithm) can be adopted, and this algorithm is based on the polarization characteristics of radar echo signals. For the sample radar remote sensing data, analyze the signal intensity and phase information of different polarization channels (such as HH, HV, VH, VV). Assume that in a certain area, the signal intensity of the HH channel is \(I_{HH} = 12 dB\), the signal intensity of the HV channel is \(I_{HV} = 6 dB\), the signal intensity of the VH channel is \(I_{VH} = 4 dB\), and the signal intensity of the VV channel is \(I_{VV} = 10 dB\). First, according to the steps of the Freeman-Durden decomposition algorithm, calculate the power contributions of different scattering mechanisms. For example, the surface scattering power \(P_s = I_{HH}+I_{VV}-2\times\sqrt{I_{HH}\times I_{VV}}\times\cos(\varphi_{HH}-\varphi_{VV})\), where \(\varphi_{HH}\) and \(\varphi_{VV}\) are the phases of the signals in the HH channel and the VV channel respectively. Assume that \(\varphi_{HH}-\varphi_{VV} = 30\) degrees, and \(\cos(30\) degrees\()\approx0.866\). First calculate the square root part: \(\sqrt{I_{HH}\times I_{VV}}=\sqrt{12\times10}\approx10.95\). Then calculate \(2\times\sqrt{I_{HH}\times I_{VV}}\times\cos(\varphi_{HH}-\varphi_{VV}) = 2\times10.95\times0.866\approx19.04\). Then the surface scattering power \(P_s = 12 + 10 - 19.04 = 2.96 dB\). Calculate the double-bounce scattering power and the volume scattering power in the same way. Through the above calculations, radar scattering characteristic features such as the surface scattering intensity parameter and the double-bounce scattering phase difference parameter are obtained. For example, finally, the surface scattering intensity parameter of this area is calculated to be \(-10 dB\), and the double-bounce scattering phase difference parameter is \(110\) degrees.

[0097] Step S224: Perform spatial registration and scale unification processing on the karst area surface coverage type feature, the terrain undulation degree feature, and the radar scattering characteristic feature to generate the training multi-source feature set.

[0098] For example, the Geographic Information System (GIS) technology can be utilized to precisely align them spatially based on the geographic coordinate information corresponding to the above-mentioned feature data, ensuring that each feature data accurately corresponds to its geographical location. For scale unification, since the resolutions of different feature data may vary, resampling methods are adopted. For example, the resolution of the surface cover type feature data in the karst area is 5 meters, the resolution of the terrain undulation feature data is 10 meters, and the resolution of the radar scattering characteristic feature data is 7.5 meters. Resample the terrain undulation feature data and the radar scattering characteristic feature data to the same 5-meter resolution as the surface cover type feature data in the karst area. Taking the terrain undulation feature data as an example, the bilinear interpolation algorithm is used for resampling. For the new grid points at the target resolution, the new altitude value is calculated by weighted averaging the altitude values of the four original grid points around it. Assume that the new grid point is located between the original grid points A (altitude 120 meters), B (altitude 125 meters), C (altitude 130 meters), and D (altitude 128 meters), and the altitude value of the new grid point is calculated through the bilinear interpolation formula. After such processing, all feature data are unified in space and scale to generate a multi-source feature set for training.

[0099] Step S230, construct an initial recognition model, where the initial recognition model includes a feature encoding layer, an attention fusion layer, and a prediction layer.

[0100] Step S240, iteratively train the initial recognition model based on the training multi-source feature set and the linear structure annotation data set until the error between the prediction result output by the prediction layer and the linear structure annotation data set reaches the convergence threshold, and generate the pre-trained linear structure recognition model.

[0101] In a possible implementation manner, step S240 includes:

[0102] Step S241, input the training multi-source feature set into the feature encoding layer to respectively generate a training optical feature encoding vector, a training terrain feature encoding vector, and a training radar feature encoding vector.

[0103] In this embodiment, the feature encoding layer may adopt a convolutional neural network (CNN) structure, which contains multiple convolutional kernels inside and is used to extract and encode the input multi-source feature set for training. For example, for the karst area surface coverage type feature in the multi-source feature set for training, the size of the convolutional kernel in the feature encoding layer is 3×3, and the feature data is scanned through convolutional operations to extract feature information at different scales. Suppose the first convolutional kernel performs a convolutional operation with a 3×3 area in the karst area surface coverage type feature data, and the information such as the vegetation coverage index, the rock exposure rate parameter, and the water body distribution density parameter in this area is weighted and summed, and then processed through an activation function (such as the ReLU function) to obtain a new feature value. Such convolutional operations are performed on the entire karst area surface coverage type feature data to generate a training optical feature encoding vector. Using the same method, the terrain undulation degree feature and the radar scattering characteristic feature are processed to generate a training terrain feature encoding vector and a training radar feature encoding vector respectively.

[0104] Step S242: Invoke the attention fusion layer to perform dynamic weight allocation processing on the training optical feature encoding vector, the training terrain feature encoding vector, and the training radar feature encoding vector to generate a training fusion feature vector.

[0105] In a possible implementation manner, step S242 includes:

[0106] Step S2421: Construct a training first attention weight matrix based on the spatial correlation between the training optical feature encoding vector and the training terrain feature encoding vector.

[0107] In this embodiment, the training optical feature encoding vector is obtained by processing the optical remote sensing image of the sample karst area through the feature encoding layer, and the dimension is set to 512. Each element therein carries the encoding information of the karst area surface coverage type feature in different feature dimensions. For example, the element value of the 10th dimension is 0.6, which is obtained through a series of operations of the feature encoding layer on various information such as the vegetation coverage index, the rock exposure rate parameter, and the water body distribution density parameter, and represents the comprehensive feature representation of the above surface coverage type feature in this feature dimension. The dimension of the training terrain feature encoding vector is 256, which is the result of processing the sample terrain elevation model data through the feature encoding layer. Taking the 15th dimension as an example, its element value is 0.8, which reflects the feature information of the terrain undulation degree of a certain area in this feature dimension and is encoded through processing such as slope calculation and terrain change analysis of the terrain elevation data.

[0108] When calculating the correlation between these two vectors in each dimension, for the $i$-th element of the optical feature dimension and the $j$-th element of the terrain feature dimension. First, find the product of the corresponding dimension elements of vector $\mathbf{a}$ (the training optical feature encoding vector) and vector $\mathbf{b}$ (the training terrain feature encoding vector), that is, perform the operation of multiplying the $i$-th element of $\mathbf{a}$ by the $j$-th element of $\mathbf{b}$. For example, when $i = 1$ and $j = 1$, the value of the first element $a_1$ of $\mathbf{a}$ is $0.5$, and the value of the first element $b_1$ of $\mathbf{b}$ is $0.3$. Then $a_1$ is multiplied by $b_1$, that is, $0.5$ multiplied by $0.3$, resulting in $0.15$. Then, perform such multiplication operations for all dimension combinations, and finally form a $512\times256$ product matrix.

[0109] After that, calculate the norm of vector $\mathbf{a}$. First, square each element of vector $\mathbf{a}$. For example, for the first three elements $a_1 = 0.5$, $a_2 = 0.6$, and $a_3 = 0.7$ of vector $\mathbf{a}$, then the square of $a_1$ is $0.5$ multiplied by $0.5$ equal to $0.25$, the square of $a_2$ is $0.6$ multiplied by $0.6$ equal to $0.36$, and the square of $a_3$ is $0.7$ multiplied by $0.7$ equal to $0.49$. Add up the squared values of all elements, and then take the square root of this sum to obtain the norm of vector $\mathbf{a}$. Calculate the norm of vector $\mathbf{b}$ in the same way.

[0110] Next, calculate the cosine similarity. Divide each element in the previously obtained product matrix by the product of the norm of vector $\mathbf{a}$ and the norm of vector $\mathbf{b}$. For example, for the product $0.15$ of $a_1$ and $b_1$ calculated above, assuming the norm of vector $\mathbf{a}$ is calculated as $1.5$ and the norm of vector $\mathbf{b}$ is calculated as $1.2$, then $0.15$ is divided by the product of $1.5$ and $1.2$, that is, $0.15$ divided by ($1.5$ multiplied by $1.2$). $1.5$ multiplied by $1.2$ equals $1.8$, and $0.15$ divided by $1.8$ is approximately equal to $0.083$. In this way, a $512\times256$ cosine similarity matrix is obtained.

[0111] Finally, after processing by the Softmax function, the Softmax function will transform each element in the cosine similarity matrix so that the sum of the elements in each row is $1$. For example, the cosine similarity values of a certain row are $[0.083, 0.12, 0.09]$. First, take the exponential of each element. The exponential value of $0.083$, $e^{0.083}$, is approximately $1.087$, the exponential value of $0.12$, $e^{0.12}$, is approximately $1.127$, and the exponential value of $0.09$, $e^{0.09}$, is approximately $1.094$. Add up these three exponential values, $1.087 + 1.127 + 1.094 = 3.308$. Then, the first element becomes $1.087$ divided by $3.308$, approximately equal to $0.329$; the second element becomes $1.127$ divided by $3.308$, approximately equal to $0.341$; the third element becomes $1.094$ divided by $3.308$, approximately equal to $0.330$. After such processing, the training first attention weight matrix is generated.

[0112] Step S2422: Construct a training second attention weight matrix based on the scattering characteristic correlation between the training radar feature encoding vector and the training terrain feature encoding vector.

[0113] For example, the dimension of the training radar feature encoding vector is 384, which contains the scattering characteristic feature encoding information of the sample radar remote sensing data processed by the feature encoding layer. Taking the 20th dimension as an example, its element value is 0.7, which is the feature representation in this feature dimension after processing information such as surface scattering intensity and secondary scattering phase difference in radar scattering characteristics.

[0114] When calculating the correlation between these two vectors for elements related to scattering characteristics, an algorithm based on distance metric is adopted. Taking the m-th radar feature dimension element and the n-th terrain feature dimension element related to scattering characteristics as an example. Calculating their Euclidean distance is to calculate the square root of the sum of the squares of the differences between these two elements. For example, when m = 1 and n = 1, the m-th radar feature dimension element c1 is 0.6, and the n-th terrain feature dimension element d1 is 0.5. Then their difference is 0.6 minus 0.5 equals 0.1, the square of 0.1 is 0.01, and the square root of 0.01 is 0.1, that is, the Euclidean distance is 0.1. Such calculations are performed for all relevant dimension combinations to obtain a 384×256 Euclidean distance matrix.

[0115] Then, the distance value is converted into a similarity value through a specific conversion function. Here, the similarity is equal to 1 divided by (1 plus the Euclidean distance). For the previously calculated Euclidean distance of 0.1, 1 plus 0.1 equals 1.1, and 1 divided by 1.1 is approximately 0.909, thus obtaining a 384×256 similarity matrix.

[0116] After further processing similar to the Softmax function, the elements in the similarity matrix are converted into probability values so that the sum of the elements in each row is 1, thereby generating a training second attention weight matrix.

[0117] Step S2423: Perform weighted fusion on the training optical feature encoding vector and the training terrain feature encoding vector according to the training first attention weight matrix to generate a training first fusion feature.

[0118] For example, the first element a1 of the trained optical feature encoding vector is 0.5, the first element b1 of the trained terrain feature encoding vector is 0.3, and the weight w11 at the corresponding position in the trained first attention weight matrix is 0.4. When performing weighted fusion, first calculate the weighted value of the elements of the trained optical feature encoding vector, that is, a1 multiplied by w11, 0.5 multiplied by 0.4 equals 0.2; then calculate the weighted value of the elements of the trained terrain feature encoding vector, that is, b1 multiplied by w11, 0.3 multiplied by 0.4 equals 0.12. Then add these two weighted values, 0.2 plus 0.12 equals 0.32, which is the value of the first element after fusion. Calculate all dimensions in this way, and finally generate a trained first fusion feature vector with a dimension of 256.

[0119] Step S2424: Perform weighted fusion on the trained radar feature encoding vector and the trained terrain feature encoding vector according to the trained second attention weight matrix to generate a trained second fusion feature.

[0120] For example, the first element c1 of the trained radar feature encoding vector is 0.6, the first element b1 of the trained terrain feature encoding vector is 0.3, and the weight w11' at the corresponding position in the trained second attention weight matrix is 0.5. First calculate the weighted value of the elements of the trained radar feature encoding vector, that is, c1 multiplied by w11', 0.6 multiplied by 0.5 equals 0.3; then calculate the weighted value of the elements of the trained terrain feature encoding vector, that is, b1 multiplied by w11', 0.3 multiplied by 0.5 equals 0.15. Add the two weighted values, 0.3 plus 0.15 equals 0.45, which is the value of the first element after fusion. Perform the same calculation for all dimensions to generate a trained second fusion feature vector with a dimension of 256.

[0121] Step S2425: Perform cross-validation processing on the trained first fusion feature and the trained second fusion feature to generate the trained fusion feature vector.

[0122] For example, concatenate the two fusion feature vectors in sequence into a vector with a dimension of 512. Then, through some verification algorithms, such as performing multiple classification predictions on the concatenated vector, comparing the prediction results with the known correct classifications, and analyzing the accuracy of the predictions, etc., process the concatenated vector, and finally generate a trained fusion feature vector.

[0123] Step S243: Input the trained fusion feature vector into the prediction layer to generate predicted linear structure spatial distribution data and predicted structure type labels.

[0124] In this embodiment, the prediction layer can adopt a fully connected neural network structure, and the internal weight matrix performs matrix multiplication operations with the training fusion feature vector. For example, the dimension of the weight matrix is 256×512 (assuming the dimension of the training fusion feature vector is 256), multiplying the training fusion feature vector by the weight matrix to obtain an output vector with a dimension of 512. Then, through the processing of an activation function (such as the Sigmoid function), each element in the output vector is converted into a probability value between 0 and 1. Then, classification judgment is made according to the above probability values. For example, if the probability value at a certain position is greater than 0.5, it is judged as a certain structural type (such as "fault"), and if it is less than 0.5, it is judged as other types (such as "fracture"). At the same time, by analyzing certain specific elements in the output vector, the predicted linear structure spatial distribution data is obtained, which is represented in the form of coordinate points of a series of line segments.

[0125] Step S244, construct a joint loss function according to the spatial position error between the predicted linear structure spatial distribution data and the linear structure annotation data set and the classification error between the predicted structure type label and the annotation type, and update the parameters of the initial recognition model based on the gradient descent algorithm until the joint loss function reaches the convergence threshold.

[0126] In a possible implementation manner, step S244 includes:

[0127] Step S2441, calculate the average Euclidean distance between the endpoints of each linear structure line segment in the predicted linear structure spatial distribution data and the corresponding annotation line segment endpoints in the linear structure annotation data set, and generate a spatial position error term.

[0128] Suppose there is a line segment in the predicted linear structure spatial distribution data, with a starting point coordinate of (700, 1000) and an ending point coordinate of (1400, 1700); the starting point coordinate of the corresponding line segment in the linear structure annotation data set is (720, 1020), and the ending point coordinate is (1380, 1680).

[0129] First, calculate the Euclidean distance of the starting point. According to the calculation method of the Euclidean distance, it is to calculate the square root of the sum of the square of the difference in abscissa and the square of the difference in ordinate. The difference in abscissa is 720 minus 700 equals 20, and the square of 20 is 400; the difference in ordinate is 1020 minus 1000 equals 20, and the square of 20 is 400. Adding the two together, 400 plus 400 equals 800, and the square root of 800 is approximately 28.28, which is the Euclidean distance of the starting point.

[0130] Then calculate the Euclidean distance of the end point. The difference in the abscissa is 1380 minus 1400 which equals -20, and the square of -20 is 400. The difference in the ordinate is 1680 minus 1700 which equals -20, and the square of -20 is 400. Adding the two together, 400 plus 400 equals 800, and the square root of 800 is approximately 28.28, which is the Euclidean distance of the end point.

[0131] Add the Euclidean distances of the starting point and the end point. 28.28 plus 28.28 equals 56.56, and then divide by 2. 56.56 divided by 2 equals 28.28. This result is the average Euclidean distance, which is also the spatial position error term.

[0132] Step S2442: Calculate the cross-entropy loss between the predicted structure type label and the annotation type to generate a classification error term.

[0133] Suppose the probability that the predicted structure type label is judged as "fault" is 0.7, while the actual probability that the annotation type is "fault" is 1. According to the calculation method of cross-entropy loss, it is the sum of the negative of the actual probability multiplied by the natural logarithm of the predicted probability. For the "fault" category, the actual probability is 1 and the predicted probability is 0.7. First, calculate the natural logarithm of 0.7, which is approximately -0.357. 1 multiplied by -0.357 equals -0.357, and then take the negative, -(-0.357) equals 0.357, which is the classification error term.

[0134] Step S2443: Weightedly sum the spatial position error term and the classification error term according to a preset ratio to generate the joint loss function.

[0135] Suppose the weight of the preset spatial position error term is 0.6 and the weight of the classification error term is 0.4. First, calculate the spatial position error term multiplied by its weight. 28.28 multiplied by 0.6, 28.28 multiplied by 0.6 equals 16.968. Then calculate the classification error term multiplied by its weight. 0.357 multiplied by 0.4, 0.357 multiplied by 0.4 equals 0.1428. Finally, add the two together. 16.968 plus 0.1428 equals 17.1108, which is the value of the joint loss function.

[0136] Step S2444: During each iteration training process, calculate the gradients of the joint loss function with respect to the parameters in the feature encoding layer, the attention fusion layer, and the prediction layer.

[0137] Taking the parameter p of a convolutional kernel in the feature encoding layer as an example, through the backpropagation algorithm, starting from the joint loss function, the error is gradually propagated to the feature encoding layer. Suppose the partial derivative of the joint loss function L with respect to a certain feature value z output by the feature encoding layer is -0.1. This means that the rate of change of the joint loss function L with respect to the feature value z is -0.1, that is, for every one-unit increase in the feature value z, the joint loss function L will decrease by 0.1. And the partial derivative of the feature value z with respect to the parameter p is 0.2, that is, for every one-unit increase in the parameter p, the feature value z will increase by 0.2.

[0138] According to the chain rule, the gradient of the joint loss function with respect to the parameter p is the partial derivative of the joint loss function with respect to the feature value z multiplied by the partial derivative of the feature value z with respect to the parameter p. That is, -0.1 multiplied by 0.2, -0.1 multiplied by 0.2 equals -0.02, thus obtaining the gradient of the joint loss function with respect to the parameter p. Such gradient calculations are performed for all parameters in the feature encoding layer, the attention fusion layer, and the prediction layer.

[0139] Step S2445, adjust the numerical values of each parameter according to the gradient direction, and update the learning rate parameter using a learning rate decay strategy.

[0140] For example, assume the current learning rate is 0.005. For the parameter p with the calculated gradient of -0.02 above, the original parameter value is 0.3. The method to adjust the parameter value is to subtract the product of the learning rate and the gradient from the original parameter value, that is, 0.3 minus (0.005 multiplied by -0.02). 0.005 multiplied by -0.02 equals -0.0001, 0.3 minus -0.0001 equals 0.3 plus 0.0001 equals 0.3001, thus obtaining the adjusted parameter value.

[0141] At the same time, adopt a learning rate decay strategy. For example, every 80 times of iterative training, multiply the learning rate by a decay factor of 0.9. Assume that 80 times of iterative training have been performed currently, then the new learning rate is 0.005 multiplied by 0.9. 0.005 multiplied by 0.9 equals 0.0045, thus obtaining the new learning rate.

[0142] Step S2446, when the decrease amplitude of the joint loss function in consecutive set numbers of iterative training is less than a preset threshold, determine that the convergence threshold is reached and terminate the training, and generate the pre-trained linear structure recognition model.

[0143] For example, when the decrease in the combined loss function during consecutive setting iterations (e.g., 60 times) is less than a preset threshold (e.g., 0.002), it is determined that the convergence threshold has been reached and the training is terminated, generating a pre-trained linear structure recognition model. For example, in a certain round of iterative training, the value of the combined loss function decreases from 17.1108 to 17.109, and the decrease amplitude is 17.1108 minus 17.109 which equals 0.0018; in the next round of iterative training, the value of the combined loss function decreases from 17.109 to 17.1075, and the decrease amplitude is 17.109 minus 17.1075 which equals 0.0015. When the decrease amplitude is less than 0.002 in 60 consecutive such iterative trainings, it is considered that the combined loss function has reached the convergence threshold. At this time, after multiple rounds of parameter adjustment and optimization, the initial recognition model can already process the input multi-source feature set for training well, generating accurate predicted linear structure spatial distribution data and predicted structure type labels, thus completing the training process of the pre-trained linear structure recognition model. This pre-trained model can be used for subsequent linear structure recognition tasks in the target area.

[0144] Figure 2 FIG. shows a schematic diagram of exemplary hardware and software components of a karst area linear structure recognition system 100 based on multi-source remote sensing that can implement the ideas of the present application provided by some embodiments of the present application. For example, the processor 120 can be used on the karst area linear structure recognition system 100 based on multi-source remote sensing and is used to execute the functions in the present application.

[0145] The karst area linear structure recognition system 100 based on multi-source remote sensing can be a general-purpose server or a special-purpose server, both of which can be used to implement the karst area linear structure recognition method of the present application. Although only one server is shown in the present application, for convenience, the functions described in the present application can be implemented in a distributed manner on multiple similar platforms to balance the processing load.

[0146] For example, the karst area linear structure recognition system 100 based on multi-source remote sensing can include a network port 110 connected to the network, one or more processors 120 for executing program instructions, a communication bus 130, and different forms of storage media 140, such as disks, ROM, or RAM, or any combination thereof. Exemplarily, the karst area linear structure recognition system 100 based on multi-source remote sensing can also include program instructions stored in ROM, RAM, or other types of non-transitory storage media, or any combination thereof. The method of the present application can be implemented according to the above program instructions. The karst area linear structure recognition system 100 based on multi-source remote sensing also includes an input / output (I / O) interface 150 between the computer and other input / output devices.

[0147] For ease of description, only one processor is described in the multi-source remote sensing-based linear structure recognition system 100 for karst areas. However, it should be noted that the multi-source remote sensing-based linear structure recognition system 100 in this application may also include multiple processors. Therefore, the steps performed by one processor described in this application can also be jointly performed or separately performed by multiple processors. For example, if the processor of the multi-source remote sensing-based linear structure recognition system 100 for karst areas performs steps A and B, it should be understood that steps A and B can also be jointly performed by two different processors or separately performed in one processor. For example, the first processor performs step A, the second processor performs step B, or the first processor and the second processor jointly perform steps A and B.

[0148] In addition, an embodiment of the present invention further provides a readable storage medium, in which computer-executable instructions are preset. When the processor executes the computer-executable instructions, the above multi-source remote sensing-based linear structure recognition method for karst areas is implemented.

[0149] It should be noted that, in order to simplify the description of the present invention disclosure and thus help the understanding of one or more embodiments of the present invention, in the previous description of the embodiments of the present invention, sometimes multiple features are merged into one embodiment, drawing or description thereof.

Claims

1. A method for identifying linear structures in karst areas based on multi-source remote sensing, characterized in that, The method includes: Obtaining optical remote sensing images, radar remote sensing data, and terrain elevation model data of the target area, and constructing a target remote sensing data set; Performing surface coverage feature enhancement processing on the target remote sensing data set to generate an enhanced multi-source feature set; the enhanced multi-source feature set includes karst area surface coverage type features, terrain undulation features, and radar scattering characteristic features; Invoking a pre-trained linear structure recognition model to perform linear structure recognition on the enhanced multi-source feature set to generate an initial linear structure recognition result; the initial linear structure recognition result includes spatial distribution data of linear structures and structure type labels; Performing spatial topological relationship verification processing on the initial linear structure recognition result based on the geological background data of the target area to generate an optimized linear structure recognition result; the optimized linear structure recognition result is used to output the spatial positions and structure confidence parameters of faults and fractures in the karst area.

2. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 1, characterized in that, The performing surface coverage feature enhancement processing on the target remote sensing data set to generate an enhanced multi-source feature set includes: Performing radiometric correction processing on the optical remote sensing image to obtain a radiometrically corrected standardized optical image, and extracting karst area surface coverage type features of the standardized optical image; the karst area surface coverage type features include vegetation coverage index, rock exposure rate parameter, and water body distribution density parameter; Performing slope calculation processing on the terrain elevation model data to generate terrain undulation features, and performing polarization decomposition processing on the radar remote sensing data to extract radar scattering characteristic features; the radar scattering characteristic features include surface scattering intensity parameter and secondary scattering phase difference parameter; Performing spatial registration and scale unification processing on the karst area surface coverage type features, the terrain undulation features, and the radar scattering characteristic features to generate the enhanced multi-source feature set.

3. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 1, wherein The invoking a pre-trained linear structure recognition model to perform linear structure recognition on the enhanced multi-source feature set to generate an initial linear structure recognition result includes: Inputting the enhanced multi-source feature set into the feature encoding layer of the pre-trained linear structure recognition model to respectively generate an optical feature encoding vector, a terrain feature encoding vector, and a radar feature encoding vector; Constructing a first attention weight matrix based on the spatial correlation between the optical feature encoding vector and the terrain feature encoding vector, and constructing a second attention weight matrix based on the scattering characteristic correlation between the radar feature encoding vector and the terrain feature encoding vector; Performing weighted fusion on the optical feature encoding vector and the terrain feature encoding vector according to the first attention weight matrix to generate a first fusion feature, and performing weighted fusion on the radar feature encoding vector and the terrain feature encoding vector according to the second attention weight matrix to generate a second fusion feature; Performing cross-validation processing on the first fusion feature and the second fusion feature to generate the initial linear structure recognition result including spatial distribution data of linear structures and structure type labels.

4. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 3, wherein Performing spatial topological relationship verification processing on the initial linear structure recognition result based on the geological background data of the target area to generate an optimized linear structure recognition result, including: Extracting the known fault zone distribution data and rock stratum attitude parameters in the geological background data to construct a geological prior knowledge base; Performing spatial overlay analysis on the initial linear structure recognition result and the geological prior knowledge base to determine the spatial consistency parameter between the linear structure spatial distribution data and the known fault zone and the rock stratum attitude matching degree parameter; Adjusting the structural confidence parameter in the initial linear structure recognition result according to the spatial consistency parameter and the rock stratum attitude matching degree parameter to generate the optimized linear structure recognition result; wherein, the structural confidence parameter is used to characterize the matching degree between the linear structure recognition result and the geological background data.

5. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 1, characterized in that, The training method of the pre-trained linear structure recognition model includes: Obtaining a multi-source remote sensing training data set of a sample karst area and a corresponding linear structure annotation data set; the multi-source remote sensing training data set includes sample optical remote sensing images, sample radar remote sensing data, and sample terrain elevation model data; Performing surface coverage feature enhancement processing on the multi-source remote sensing training data set to generate a training multi-source feature set; the training multi-source feature set includes karst area surface coverage type features, terrain undulation features, and radar scattering characteristic features; Constructing an initial recognition model, the initial recognition model includes a feature encoding layer, an attention fusion layer, and a prediction layer; Performing iterative training on the initial recognition model based on the training multi-source feature set and the linear structure annotation data set until the error between the prediction result output by the prediction layer and the linear structure annotation data set reaches a convergence threshold to generate the pre-trained linear structure recognition model.

6. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 5, characterized in that, The performing surface coverage feature enhancement processing on the multi-source remote sensing training data set to generate a training multi-source feature set includes: Performing radiometric correction processing on the sample optical remote sensing image to obtain a standardized sample optical image, and extracting the karst area surface coverage type features of the standardized sample optical image; Performing slope calculation processing on the sample terrain elevation model data to generate the terrain undulation features; Performing polarization decomposition processing on the sample radar remote sensing data to extract the radar scattering characteristic features; Performing spatial registration and scale unification processing on the karst area surface coverage type features, the terrain undulation features, and the radar scattering characteristic features to generate the training multi-source feature set.

7. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 5, wherein The performing iterative training on the initial recognition model based on the training multi-source feature set and the linear structure annotation data set includes: Inputting the training multi-source feature set into the feature encoding layer to respectively generate a training optical feature encoding vector, a training terrain feature encoding vector, and a training radar feature encoding vector; Invoking the attention fusion layer to perform dynamic weight allocation processing on the training optical feature encoding vector, the training terrain feature encoding vector, and the training radar feature encoding vector to generate a training fusion feature vector; Input the training fusion feature vector into the prediction layer to generate predicted linear structure spatial distribution data and predicted structure type labels; Construct a joint loss function based on the spatial position error between the predicted linear structure spatial distribution data and the linear structure annotation data set and the classification error between the predicted structure type label and the annotation type, and update the parameters of the initial recognition model based on the gradient descent algorithm until the joint loss function reaches the convergence threshold.

8. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 7, characterized in that, The step of calling the attention fusion layer to perform dynamic weight allocation processing on the training optical feature encoding vector, the training terrain feature encoding vector, and the training radar feature encoding vector to generate a training fusion feature vector includes: Construct a training first attention weight matrix based on the spatial correlation between the training optical feature encoding vector and the training terrain feature encoding vector; Construct a training second attention weight matrix based on the scattering characteristic correlation between the training radar feature encoding vector and the training terrain feature encoding vector; Perform weighted fusion on the training optical feature encoding vector and the training terrain feature encoding vector according to the training first attention weight matrix to generate a training first fusion feature; Perform weighted fusion on the training radar feature encoding vector and the training terrain feature encoding vector according to the training second attention weight matrix to generate a training second fusion feature; Perform cross-validation processing on the training first fusion feature and the training second fusion feature to generate the training fusion feature vector.

9. The method for identifying linear structures in karst areas based on multi-source remote sensing according to claim 7, characterized in that The step of constructing a joint loss function based on the spatial position error between the predicted linear structure spatial distribution data and the linear structure annotation data set and the classification error between the predicted structure type label and the annotation type includes: Calculate the average Euclidean distance between the endpoints of each linear structure segment in the predicted linear structure spatial distribution data and the endpoints of the corresponding annotation segments in the linear structure annotation data set to generate a spatial position error term; Calculate the cross-entropy loss between the predicted structure type label and the annotation type to generate a classification error term; Perform weighted summation of the spatial position error term and the classification error term according to a preset ratio to generate the joint loss function; The step of updating the parameters of the initial recognition model based on the gradient descent algorithm until the joint loss function reaches the convergence threshold includes: In each iteration training process, calculate the gradients of the joint loss function with respect to the parameters in the feature encoding layer, the attention fusion layer, and the prediction layer; Adjust the numerical values of the parameters according to the gradient direction and update the learning rate parameter using a learning rate decay strategy; When the decrease amplitude of the joint loss function in continuously set number of iteration trainings is less than the preset threshold, determine that the convergence threshold is reached and terminate the training to generate the pre-trained linear structure recognition model.

10. A linear structure identification system for karst areas based on multi-source remote sensing, characterized in that, The linear structure recognition system for karst areas based on multi-source remote sensing includes a processor and a memory. The memory is connected to the processor. The memory is used to store programs, instructions or codes, and the processor is used to execute the programs, instructions or codes in the memory to implement the linear structure recognition method for karst areas based on multi-source remote sensing according to any one of claims 1-9 above.

Images