Airborne multi-modal reliability fusion and graph reasoning for ground collapse monitoring

Abstract

The application discloses a ground collapse monitoring method based on air-based multi-modal reliability fusion and graph reasoning, and belongs to the technical field of urban geological safety monitoring. The method comprises the following steps: collecting visible light images, thermal infrared images and radar point clouds / elevation data of a city road area; performing quality assessment on standardized input data; constructing a feature coding branch, inputting the standardized input data into the corresponding feature coding branch, and obtaining three-modal features; performing weighted interaction and fusion on the three-modal features to obtain fusion features; constructing a geometric-thermal anomaly topology graph according to the fusion features, and performing message passing and risk propagation updating through a graph reasoning mechanism; performing risk decoding and output on the graph reasoning result to generate a pixel-level risk probability graph, a collapse hotspot mask and a target area sorting result. The application can improve the accuracy and robustness of city ground collapse hazard identification, and provides technical support for city road collapse risk perception.

Description

Technical Field

[0001] This invention relates to the field of urban geological safety monitoring technology, and in particular to a ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning. Background Technology

[0002] Urban ground subsidence is a common geological safety hazard during the construction of urban roads, underground pipe networks, underground space development, and rail transit. Ground subsidence is typically related to factors such as underground cavities, roadbed voids, pipeline leaks, loose soil, groundwater disturbance, and long-term traffic loads. It is characterized by its high degree of concealment, complex development process, and sudden and severe consequences. Once a ground subsidence occurs, it can easily cause serious consequences such as traffic disruption, personal injury, vehicle damage, and damage to underground pipelines. Therefore, rapid, accurate, and large-scale monitoring of subsidence hazards in urban road areas is of great significance.

[0003] Currently, urban ground subsidence hazard monitoring mainly relies on methods such as manual inspections, vehicle-mounted detection, ground radar detection, the deployment of settlement monitoring points, and single remote sensing image recognition. Manual inspection methods are greatly affected by subjective experience, resulting in low detection efficiency and difficulty in timely detection of early micro-cracks and localized settlement phenomena. While refined detection methods such as ground radar can obtain information on underground structures, they typically require the deployment of survey lines on the ground area by area, leading to high detection costs. Furthermore, their operational efficiency is greatly affected by road traffic conditions, making it difficult to meet the needs of rapid screening of large-scale urban roads.

[0004] In recent years, UAV remote sensing technology has been increasingly applied to urban road inspection and disaster hazard investigation due to its advantages such as high mobility, wide coverage, and high data acquisition efficiency. By equipping UAVs with visible light cameras, thermal infrared cameras, and radar sensors, information on road surface texture, temperature anomalies, and three-dimensional elevation structure can be acquired. Visible light images can reflect cracks, potholes, and road surface morphology; thermal infrared images can reflect temperature changes caused by water seepage, differences in thermal conductivity due to cavities, or material anomalies; and radar point cloud / elevation data can reflect road surface elevation depressions, roughness changes, and local geometric deformations. The strong complementarity among these multi-source data provides a new technological foundation for identifying urban ground subsidence hazards.

[0005] However, existing space-based multimodal monitoring methods still have the following problems: On the one hand, different modal data are easily affected by different interferences in complex urban road environments. For example, visible light images are easily affected by shadows, water stains, strong light, and occlusion; thermal infrared images are easily affected by heat reflection and changes in ambient temperature; and radar point clouds are easily affected by sparse point clouds, missing data, and noise. On the other hand, existing methods mostly adopt simple stitching, channel fusion, or post-decision fusion methods, lacking explicit evaluation of the local reliability of different modalities, which can easily lead to low-quality modalities misleading the identification results. In addition, urban ground subsidence risk usually has spatial continuity and regional propagation characteristics, and simply relying on pixel-level detection or local convolution features is difficult to fully express the spatial correlation between cracks, thermal anomalies, and elevation depressions.

[0006] Therefore, how to conduct reliability assessments of visible light, thermal infrared, and radar data based on UAV multimodal data, and how to adaptively fuse different modes according to their credibility in different regions, while combining geometric and thermal anomaly spatial topological relationships to achieve ground subsidence risk inference, are urgent technical problems that need to be solved in current urban ground subsidence monitoring. Summary of the Invention

[0007] Purpose of the invention: To address the above problems, the purpose of this invention is to provide a ground subsidence monitoring method based on spaceborne multimodal reliability fusion and graph reasoning. This method uses UAVs to simultaneously collect visible light images, thermal infrared images, and radar point cloud / elevation data, performs quality assessment on each modal data, and generates a modal reliability map. Furthermore, it utilizes a reliability-guided cross-modal attention fusion mechanism and a geometric-thermal anomaly topological graph reasoning mechanism to achieve rapid identification of urban ground subsidence hazards, risk propagation analysis, and target area ranking output.

[0008] Technical solution: The ground subsidence monitoring method based on spaceborne multimodal reliability fusion and graph reasoning of the present invention includes the following steps:

[0009] Step 1: Collect visible light images, thermal infrared images, and radar point cloud / elevation data of urban road areas, and preprocess the multimodal data to obtain standardized input data;

[0010] Step 2: Perform quality assessment on the standardized input data, calculate modal quality terms, registration error terms, and environmental interference terms, and generate reliability maps for each modality;

[0011] Step 3: Construct three feature coding branches with independent parameters, and input the standardized input data into the corresponding feature coding branches to obtain multi-scale visible light texture features, thermal anomaly features, and radar geometric features;

[0012] Step 4: Based on the cross-modal attention fusion mechanism guided by modal reliability, the three-modal features are weighted and fused to obtain the fused features;

[0013] Step 5: Construct a geometric-thermal anomaly topology map based on the fusion features, and update it through message passing and risk propagation via graph reasoning mechanism;

[0014] Step 6: Decode and output the graph reasoning results to generate a pixel-level risk probability map, a collapsed hotspot mask, and a target area ranking result.

[0015] Preferably, step 1 includes:

[0016] Using a drone platform equipped with a visible light camera, a thermal infrared camera, and a radar sensor, simultaneous aerial data collection was conducted on the urban road area to be monitored, obtaining the first... Multimodal data set corresponding to this flight ,in, Indicates the first Visible light images obtained during this flight. Indicates the first Thermal infrared images obtained during this flight. Indicates the first Radar point cloud / elevation data obtained during this flight , and These represent the visible light mode, the thermal infrared mode, and the lidar mode, respectively. This indicates the drone's pose, sensor intrinsic and extrinsic parameters, and georeferenced information. This represents the data acquisition timestamp, which is used to synchronize and match the three-modal data and associate data acquired from the same flight with the same flight sequence number. Indicates the flight sequence number;

[0017] For a point in three-dimensional space In modal The projection position in the image is represented as:

[0018] ;

[0019] in, Represents a set of modal identifiers. Indicates modal identifier, , Represents the homogeneous coordinate vector of a point in three-dimensional space. , and Indicates three-dimensional coordinate components, superscript Indicates transpose. This indicates that the point is in mode The two-dimensional projection position in , and Representing modes The corresponding intrinsic parameter matrix, rotation matrix, and translation vector of the sensor. Representing modes The projection function;

[0020] Using the visible light image coordinate system as the reference coordinate system, thermal infrared imagery and radar point cloud / elevation data are mapped to a unified road grid coordinate system, represented as follows:

[0021] ;

[0022] ;

[0023] in, This represents the spatial transformation function from the thermal infrared image coordinate system to the visible light reference coordinate system. This represents the spatial transformation function from the radar coordinate system to the visible light reference coordinate system. and These represent the mapped thermal infrared image and radar point cloud / elevation data, respectively.

[0024] The mapped multimodal data is then subjected to denoising, normalization, scale unification, and image segmentation to obtain a standardized input data set.

[0025] ;

[0026] in, Represents a standardized set of input data. , and These represent visible light images, thermal infrared images, and radar point cloud / elevation data after registration, denoising, normalization, and scale unification, respectively.

[0027] Preferably, step 2 includes:

[0028] For the visible light mode, its quality term is jointly determined by image sharpness, local texture information, and exposure anomalies, and is expressed as:

[0029] ;

[0030] in, Indicates the visible light mode at position The mass value at that location, Indicates pixel position or road grid position. This indicates that the input is normalized to... Functions on intervals Represents local variance. Represents the Laplace operator. Indicates the location of standardized visible light images Pixel value at that location, Indicates overexposure or underexposure interference;

[0031] For the thermal infrared mode, its mass term is determined by the contrast of local thermal anomalies, and is expressed as:

[0032] ;

[0033] in, Indicates the thermal infrared mode at position The mass value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates position The neighborhood, Indicates the number of pixels in the neighborhood. Representing the neighborhood Pixel position index within;

[0034] For radar modes, the mass term is determined by both point cloud density and elevation stability, and is expressed as:

[0035] ;

[0036] in, Indicates the radar mode at position The mass value at that location, Indicates position The number of point clouds within a corresponding raster or voxel. This represents the maximum number of point clouds within a local statistical range. This represents the standard deviation of point cloud elevation within the location's neighborhood. This represents the positive scaling parameter that controls the rate of decay of elevation dispersion. Represents the natural exponential function;

[0037] modality In position The registration error term at the location is defined as:

[0038] ;

[0039] in, Representing modes In position Registration error at the location, Representing modes Center and reference position The corresponding position Indicates modality Spatial transformation function from the coordinate system to the visible light reference coordinate system. Represents the Euclidean norm;

[0040] modality The environmental disturbance term is represented as:

[0041] ;

[0042] in, Representing modes In position The amount of environmental disturbance at the location , , and These represent shadows, reflections, sparse point clouds, and noise interference, respectively. , , and These are the non-negative weight coefficients for the corresponding interference terms;

[0043] Calculate the modal quality term, registration error term, and environmental disturbance term. In position The reliability at the location is given by the formula:

[0044] ;

[0045] in, Representing modes In position Reliability before normalization at the location, , and These represent the modal quality term, registration error term, and environmental disturbance term, respectively. This represents the Sigmoid function. , and It is a non-negative scaling factor. These are bias parameters;

[0046] The reliability of each mode is normalized to obtain a normalized reliability map, which is represented as follows:

[0047] ;

[0048] in, Represents the normalized mode reliability, Indicates the modal summation index. This represents a very small positive number that prevents the denominator from being zero.

[0049] Preferably, step 3 includes:

[0050] Standardize visible light images The input is an RGB texture encoding branch; the RGB texture encoding branch includes a series of convolutional feature extraction units connected in sequence, with a downsampling unit set between adjacent convolutional feature extraction units. Each convolutional feature extraction unit outputs a visible light feature at one scale, represented as:

[0051] ;

[0052] in, Indicates the first Visible light characteristics at various scales This represents a multi-scale feature hierarchy index. This indicates that the RGB texture encoding branch is in the first... Feature extraction functions at various scales;

[0053] Each convolutional feature extraction unit includes a convolutional layer, a normalization layer, and a nonlinear activation layer. The downsampling unit employs stride convolution or pooling operations.

[0054] To enhance the information on road crack boundaries and surface morphology, an edge enhancement feature is introduced, represented as:

[0055] ;

[0056] in, Indicates the first Edge enhancement features at various scales Represents a non-linear activation function. This represents the convolution operation. This indicates the edge extraction operator. This indicates feature splicing based on channels;

[0057] Edge enhancement features are fused with visible light multi-scale features to obtain updated visible light texture features. , is represented as:

[0058] ;

[0059] Standardized thermal infrared images With local thermal anomaly map The thermal anomaly coding branch is input after concatenation by channel; the thermal anomaly coding branch includes a series of convolutional feature extraction units connected in sequence, with a downsampling unit set between adjacent convolutional feature extraction units, and each convolutional feature extraction unit outputs thermal anomaly features at one scale, represented as follows:

[0060] ;

[0061] in, Indicates the first Thermal anomaly characteristics at various scales Indicates the thermal anomaly coding branch at the 1st Feature extraction functions at various scales Indicates the first The local thermal anomaly map corresponding to the second flight, wherein the local thermal anomaly map is defined as:

[0062] ;

[0063] in, Indicates position The difference between the temperature response at a location and the average temperature response of its neighborhood is the local thermal anomaly map.

[0064] Standardized radar point cloud / elevation data Perform voxelization or road rasterization processing, and calculate local geometric features; for the first... Individual element units, with their center points represented as:

[0065] ;

[0066] in, Indicates the first A collection of point clouds within an individual element or road grid cell. Indicates voxel index, This represents the number of points in the set. Represents the set of elements. The three-dimensional coordinate vector of each radar point Represents a point index. , and These represent the horizontal and vertical coordinates of the point, respectively. Indicates the first The center point of the individual element;

[0067] No. The relative center point vector of an individual element is represented as:

[0068] ;

[0069] The standard deviation of point cloud is expressed as:

[0070] ;

[0071] The extreme value difference is expressed as:

[0072] ;

[0073] in, This represents the average offset vector of points within a voxel relative to the center point. This represents the standard deviation of point cloud coordinates. This represents the difference between the extreme values ​​of each coordinate component. and These represent taking the maximum and minimum values, respectively.

[0074] For the The local plane fitted to the neighborhood of a voxel is represented as:

[0075] ;

[0076] in, and Represents two horizontal coordinates within a local plane. Indicates the corresponding elevation. and This represents the plane slope coefficient obtained by fitting using the least squares method. Represents the plane intercept;

[0077] And calculate the local plane residual, using the formula:

[0078] ;

[0079] in, Indicates the first The average distance residual from each point within an individual element to the fitted plane. , , They represent the first Individual factors Inner Points The first horizontal coordinate component, the second horizontal coordinate component, and the elevation coordinate component in the local coordinate system; The index of a point within a voxel is used; let the eigenvalues ​​of the local covariance matrix of the point cloud be... Calculate curvature characteristics:

[0080] ;

[0081] in, , and These are the eigenvalues ​​of the local covariance matrix of the point cloud, arranged in descending order. Indicates the first Normalized curvature of individual elements;

[0082] Elevation depression features are calculated based on radar point cloud / elevation data, for location elevation value at the location The neighborhood average elevation is expressed as:

[0083] ;

[0084] Location The elevation depression characteristics at the location are represented as follows:

[0085] ;

[0086] in, Indicates neighborhood location The elevation value at that location, Indicates position The elevation value at that location, Indicates position The neighborhood average elevation, This represents the normalized local elevation depression obtained from the difference between the neighborhood average elevation and the current location elevation;

[0087] For the Individual elements or road grid cells, the elevation depression characteristics are represented as follows:

[0088] ;

[0089] in, Indicates the first The average elevation depression characteristics of individual elements or road grid cells. Indicates the first Position of each grid cell within an individual element Find the arithmetic mean;

[0090] By concatenating the geometric features, the radar voxel-enhanced geometric features are obtained, represented as follows:

[0091] ;

[0092] Enhance the geometric features of each voxel Input radar geometry coding branch: First, the enhanced geometric features of each voxel are mapped to a unified feature dimension through a multilayer perceptron mapping unit. Then, the features of spatially adjacent voxels are aggregated by a neighborhood feature aggregation unit, either by maximum value aggregation or average value aggregation. Finally, radar geometric features at different scales are formed step by step through a downsampling unit, as shown below:

[0093] ;

[0094] in, Indicates the radar geometry coding branch in the first... Feature extraction functions at various scales The radar geometric features at the corresponding scale are represented; the multilayer perceptron mapping unit includes a fully connected layer and a nonlinear activation layer, and the neighborhood feature aggregation unit adopts maximum value aggregation or average value aggregation.

[0095] Preferably, step 4 includes:

[0096] Map the three modalities of features to a unified feature space:

[0097] ;

[0098] ;

[0099] ;

[0100] in, Representing modes In the Input features at various scales; when The updated visible light texture features are used. ; , and Representing modes The query matrix, key matrix, and value matrix. , and For the corresponding learnable mapping matrix; when hour, Using thermal anomaly characteristics ;when hour, Using radar geometric features ;

[0101] Adjust the normalized reliability map to the corresponding scale, and align it with the first... The modal feature space has a consistent size across all scales, represented as:

[0102] ;

[0103] in, Representing modes In the Normalized reliability plots at various scales This indicates an interpolation scaling operation. This indicates the spatial dimensions of the corresponding scale feature map;

[0104] The normalized reliability map is incorporated as an attention bias into the cross-modal attention calculation:

[0105] ;

[0106] in, Indicates the first Modal at each scale For modes Attention contribution weight, This is the reliability adjustment factor. To query the feature dimensions of the key matrix and the key matrix, Represents the normalized exponential function, Represents the natural logarithm function;

[0107] Based on cross-modal attention weights, information from different modalities is fused to obtain the modality. Enhanced features , is represented as:

[0108] ;

[0109] Modality-level fusion weights are generated based on the global context and modality reliability, and are expressed as follows:

[0110] ;

[0111] in, Representing modes In the Fusion weights at various scales This indicates normalization in the modal dimension. This represents a multilayer perceptron. Indicates global average pooling. Indicates spatial average. This represents the reliability weighting adjustment coefficient;

[0112] The three modality enhancement features are weighted and summed to obtain the fused feature, which is represented as:

[0113] ;

[0114] By aggregating the fused features from multiple scales, a fused feature at a unified scale is obtained, represented as:

[0115] ;

[0116] in, This represents a multi-scale feature convergence operation. Indicates the total number of feature scales. Indicates the first Fusion features at various scales This represents the final fusion feature under a unified scale.

[0117] Preferably, step 5 includes:

[0118] Final fusion features The topology graph node set is obtained by dividing the area into several nodes according to the road grid or candidate risk area. ,in, Indicates the first 1 node Indicates the node index. Indicates the total number of nodes;

[0119] No. The node characteristics of each node are represented as follows:

[0120] ;

[0121] in, Indicates the first The feature vector of each node This indicates a region pooling operation. Indicates the first Each road grid area or candidate risk area , and These represent the normalized reliability of the three modes, respectively.

[0122] Calculate the first The formulas for the spatial center location, average elevation depression characteristics, and average thermal anomaly characteristics of the region corresponding to each node are as follows:

[0123] ;

[0124] ;

[0125] ;

[0126] in, Indicates the first Each node corresponds to the spatial center of the region. This indicates the number of locations within the corresponding area. This indicates the degree of average elevation depression in the corresponding area. Indicates the first The average thermal anomaly intensity of the corresponding region during this flight. Indicates Calculate the arithmetic mean for the statistical domain;

[0127] Any two nodes and The edge weights between components are determined by spatial distance, semantic similarity, thermal anomaly consistency, elevation depression continuity, and crack connectivity, and are expressed as follows:

[0128] ;

[0129] ;

[0130] ;

[0131] ;

[0132] in, Indicates spatial proximity. Indicates the similarity of node features. Indicates the continuity of elevation depressions. Indicates consistency of thermal anomalies; The second norm of a vector, also known as the Euclidean norm. Represents the square of the Euclidean norm; Represents cosine similarity. , and A positive scaling parameter for measuring correspondence; and They represent the first The node and the first One node; , , and They represent the first The spatial center location of each node, node characteristics, average elevation depression degree, and the first Average thermal anomaly intensity during the second flight;

[0133] The crack skeleton map is generated from the edge enhancement features and is represented as follows:

[0134] ;

[0135] in, This represents a crack skeleton diagram. This indicates a threshold segmentation operation. This represents the skeletonization process; the generated crack skeleton map is used to eliminate crack width differences and preserve centerline connectivity, providing a topological basis for subsequent calculations of the shortest crack path and connectivity between nodes.

[0136] Node calculation based on crack skeleton diagram With nodes The connectivity of the cracks between them is represented as:

[0137] ;

[0138] in, Represents a node With nodes The connectivity of the cracks between them This indicates that the two are in the crack skeleton diagram. The shortest path distance on the [path]. This represents the scale parameter indicating the connectivity of the crack;

[0139] The topological adjacency edge weights are generated as follows:

[0140] ;

[0141] in, Represents a node To the node The topological adjacency edge weights after Softmax normalization based on the neighbor dimension. , , , and The non-negative weight parameters are for the five types of relationships. Represents a node Graph neighbor set, Index the neighboring nodes;

[0142] Introduce self-connects and normalize the adjacency matrix:

[0143] ;

[0144] in, This represents the normalized adjacency matrix obtained after adding self-joins. Indicates by The topological adjacency matrix that has been constructed but not yet normalized for symmetry. Representation matrix The degree matrix, express An identity matrix of order 1;

[0145] Node features are updated using a graph neural network, with the following formula:

[0146] ;

[0147] in, and They represent the first The layer's input and output node feature matrices, and the initial node feature matrix. Based on the characteristics of each node composition, Indicates the first The learnable parameter matrix of the layer, Indicates the layer index of a graph neural network;

[0148] Mapping the node features obtained after graph inference back to the road grid or pixel space yields graph-enhanced risk features, represented as:

[0149] ;

[0150] in, This represents the total number of layers in a graph neural network. This represents the inverse pooling operation that backfills node features into the corresponding road grid or pixel region. Indicates position The graph at that location enhances the risk characteristics.

[0151] Preferably, step 6 includes:

[0152] A pixel-level risk probability map is generated based on graph-enhanced risk features, represented as follows:

[0153] ;

[0154] in, Indicates position The probability of collapse at the location Indicates the risk decoding weight. Indicates risk decoding bias;

[0155] Generate a collapse hotspot mask based on the risk probability map: when hour, , indicating position It belongs to the collapse hotspot area; when hour, , indicating position It is not located in a collapse hotspot area; among them, Indicates position Binary hotspot markers at the location, Indicates the risk threshold;

[0156] Connectivity analysis was performed on the collapse hotspot areas to obtain the target region set. ,in, Represents the target region set. Indicates the first A connected target region, Indicates the target region index. Indicates the total number of target areas;

[0157] And calculate the first Risk score for each target area:

[0158] ;

[0159] in, Indicates the first Risk score for each target area Indicates the target area Internal locations The arithmetic mean of the risk probabilities. Indicate target area The area;

[0160] Based on the risk score All target areas are sorted to obtain target area sorting results, and the risk priority of different suspected collapse areas is characterized by the target area sorting results.

[0161] Beneficial effects: Compared with the prior art, the significant advantages of this invention are:

[0162] 1. This invention calculates the mass terms, registration error terms, and environmental interference terms of the visible light, thermal infrared, and radar / elevation modes respectively to form a position-related modal reliability map, enabling the system to explicitly identify low-confidence areas such as shadows, water stains, thermal reflections, imaging noise, and sparse point clouds, thereby reducing the interference of abnormal observations on monitoring results.

[0163] 2. In this invention, the normalized reliability map is added as an attention bias to the cross-modal attention calculation, and the three-modal features are dynamically weighted by combining modal-level fusion weights. Compared with fixed weights or direct feature splicing, it can adaptively enhance high-confidence modalities and suppress low-confidence modalities at different road locations, thereby improving the stability of fusion results in complex environments.

[0164] 3. In this invention, RGB texture coding branches and Sobel edge enhancement are used to extract crack and surface morphology information, local thermal contrast is used to extract thermal anomaly information, and voxel centers, local planar residuals, covariance curvature and elevation depression are used to extract radar geometric information. This can form complementary evidence from three aspects: texture, temperature and three-dimensional deformation, and enhance the ability to identify early collapse signs.

[0165] 4. This invention constructs topological edges based on node spatial proximity, feature similarity, elevation continuity, thermal anomaly consistency, and crack skeleton connectivity, and propagates risk information through a graph neural network with residual connections, enabling discrete local anomalies to form continuous risk areas according to road spatial relationships, improving the integrity of hazard boundaries and reducing isolated misjudgments.

[0166] 5. This invention, through risk decoding, threshold segmentation, connected component analysis, and target area scoring, simultaneously outputs pixel-level risk probability maps, collapse hotspot masks, and target area ranking results. It can transform large-scale UAV inspection results into a priority review list, thereby improving the efficiency of urban road hazard screening and on-site review. Attached Figure Description

[0167] Figure 1 This is a flowchart of the method of the present invention;

[0168] Figure 2 This is a schematic diagram of the method of the present invention;

[0169] Figure 3 This is a comparison chart of PothRGBD elevation / geometric proxy experimental indicators;

[0170] Figure 4 This is a comparison chart of experimental indicators for IRT-Crack RGB-thermal infrared crack recognition;

[0171] Figure 5 This is a comparison chart of the key performance improvements in the two sets of experiments;

[0172] Figure 6 This is a qualitative visualization result of PothRGBD and IRT-Crack. Detailed Implementation

[0173] The embodiments of the present invention will be further described in detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the present invention and not intended to limit the scope of the invention. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the embodiments of the present invention, and not all structures.

[0174] In the following description, specific details such as target system architecture and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods are omitted so as not to obscure the description of this application with unnecessary detail.

[0175] It should be understood that, when used in this application specification and the appended claims, the term "comprising" indicates the presence of the described features, integrals, steps, operations, elements and / or components, but does not exclude the presence or addition of one or more other features, integrals, steps, operations, elements, components and / or a collection thereof.

[0176] It should also be understood that the term “and / or” as used in this application specification and the appended claims means any combination of one or more of the associated listed items and all possible combinations, and includes such combinations.

[0177] Furthermore, in the description of this application and the appended claims, the terms "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0178] References to "one embodiment" or "some embodiments" in this specification mean that one or more embodiments of this application include the target features, structures, or characteristics described in connection with that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this specification do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized.

[0179] Combination Figure 1 and Figure 2 As shown in this embodiment, the ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning includes the following steps:

[0180] Step 1: Collect visible light images, thermal infrared images, and radar point cloud / elevation data of urban road areas, and preprocess the multimodal data to obtain standardized input data.

[0181] Further, step 1 includes:

[0182] Using a drone platform equipped with a visible light camera, a thermal infrared camera, and a radar sensor, simultaneous aerial data collection was conducted on the urban road area to be monitored, obtaining the first... Multimodal data set corresponding to this flight ,in, Indicates the first Visible light images obtained during this flight. Indicates the first Thermal infrared images obtained during this flight. Indicates the first Radar point cloud / elevation data obtained during this flight , and These represent the visible light mode, the thermal infrared mode, and the lidar mode, respectively. This indicates the drone's pose, sensor intrinsic and extrinsic parameters, and georeferenced information. This represents the data acquisition timestamp. The timestamp is used to synchronize and match the three-modal data and associate data acquired from the same flight with the same flight sequence number. Indicates the flight sequence number;

[0183] For a point in three-dimensional space In modal The projection position in the image is represented as:

[0184] ;

[0185] in, Represents a set of modal identifiers. Indicates modal identifier, , Represents the homogeneous coordinate vector of a point in three-dimensional space. , and Indicates three-dimensional coordinate components, superscript Indicates transpose. This indicates that the point is in mode The two-dimensional projection position in , and Representing modes The corresponding intrinsic parameter matrix, rotation matrix, and translation vector of the sensor. Representing modes The projection function;

[0186] Using the visible light image coordinate system as the reference coordinate system, thermal infrared imagery and radar point cloud / elevation data are mapped to a unified road grid coordinate system, represented as follows:

[0187] ;

[0188] ;

[0189] in, This represents the spatial transformation function from the thermal infrared image coordinate system to the visible light reference coordinate system. This represents the spatial transformation function from the radar coordinate system to the visible light reference coordinate system, both of which are derived from the UAV's pose and sensor intrinsic and extrinsic parameters. Sure; and These represent the mapped thermal infrared image and radar point cloud / elevation data, respectively.

[0190] The mapped multimodal data is then subjected to denoising, normalization, scale unification, and image segmentation to obtain a standardized input data set.

[0191] ;

[0192] in, Represents a standardized set of input data. , and These represent visible light images, thermal infrared images, and radar point cloud / elevation data after registration, denoising, normalization, and scale unification, respectively.

[0193] Step 2: Perform quality assessment on the standardized input data, calculate modal quality terms, registration error terms, and environmental interference terms, and generate reliability maps for each modality.

[0194] Further, step 2 includes:

[0195] For the visible light mode, its quality term is jointly determined by image sharpness, local texture information, and exposure anomalies, and is expressed as:

[0196] ;

[0197] in, Indicates the visible light mode at position The mass value at that location, Indicates pixel position or road grid position. This indicates that the input is normalized to... Functions on intervals Represents local variance. Represents the Laplace operator. Indicates the location of standardized visible light images Pixel value at that location, Indicates overexposure or underexposure interference;

[0198] For the thermal infrared mode, its mass term is determined by the contrast of local thermal anomalies, and is expressed as:

[0199] ;

[0200] in, Indicates the thermal infrared mode at position The mass value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates position The neighborhood, Indicates the number of pixels in the neighborhood. Representing the neighborhood Pixel position index within;

[0201] For radar modes, the mass term is determined by both point cloud density and elevation stability, and is expressed as:

[0202] ;

[0203] in, Indicates the radar mode at position The mass value at that location, Indicates position The number of points within a corresponding raster or voxel, the number of points, and the standard deviation of elevation are both derived from standardized radar point cloud / elevation data. Calculations show that This represents the maximum number of point clouds within a local statistical range. This represents the standard deviation of point cloud elevation within the location's neighborhood. This represents the positive scaling parameter that controls the rate of decay of elevation dispersion. Represents the natural exponential function;

[0204] modality In position The registration error term at the location is defined as:

[0205] ;

[0206] in, Representing modes In position Registration error at the location, Representing modes Center and reference position The corresponding position Indicates modality Spatial transformation function from the coordinate system to the visible light reference coordinate system. Represents the Euclidean norm;

[0207] modality The environmental disturbance term is represented as:

[0208] ;

[0209] in, Representing modes In position The amount of environmental disturbance at the location , , and These represent shadows, reflections, sparse point clouds, and noise interference, respectively. , , and For the mode, the non-negative weight coefficients corresponding to the interference term are... For interference items that are not applicable, the weight coefficient is set to 0.

[0210] Calculate the modal quality term, registration error term, and environmental disturbance term. In position The reliability at the location is given by the formula:

[0211] ;

[0212] in, Representing modes In position Reliability before normalization at the location, , and These represent the modal quality term, registration error term, and environmental disturbance term, respectively. This represents the Sigmoid function. , and It is a non-negative scaling factor. These are bias parameters;

[0213] The reliability of each mode is normalized to obtain a normalized reliability map, which is represented as follows:

[0214] ;

[0215] in, Represents the normalized mode reliability, Indicates the modal summation index. This represents a very small positive number that prevents the denominator from being zero;

[0216] Step 3: Construct three feature coding branches with independent parameters, and input the standardized input data into the corresponding feature coding branches to obtain multi-scale visible light texture features, thermal anomaly features, and radar geometric features.

[0217] The three feature coding branches, which do not share parameters, include the RGB texture coding branch, the thermal anomaly coding branch, and the radar geometry coding branch.

[0218] Furthermore, step 3 includes:

[0219] Standardize visible light images The input is the RGB texture encoding branch; the RGB texture encoding branch consists of multiple convolutional feature extraction units connected in sequence, with downsampling units set between adjacent convolutional feature extraction units. Each convolutional feature extraction unit outputs a visible light feature at one scale, represented as:

[0220] ;

[0221] in, Indicates the first Visible light characteristics at various scales This represents a multi-scale feature hierarchy index. This indicates that the RGB texture encoding branch is in the first... Feature extraction functions at various scales;

[0222] Each convolutional feature extraction unit includes a convolutional layer, a normalization layer, and a non-linear activation layer. The downsampling unit uses stride convolution or pooling operations.

[0223] To enhance the information on road crack boundaries and surface morphology, an edge enhancement feature is introduced, represented as:

[0224] ;

[0225] in, Indicates the first Edge enhancement features at various scales Represents a non-linear activation function. This represents the convolution operation. This indicates the edge extraction operator. This indicates feature splicing based on channels;

[0226] Edge enhancement features are fused with visible light multi-scale features to obtain updated visible light texture features. , is represented as:

[0227] ;

[0228] Standardized thermal infrared images With local thermal anomaly map The thermal anomaly coding branch is input after concatenation by channel. The thermal anomaly coding branch consists of multiple sequentially connected convolutional feature extraction units. Downsampling units are set between adjacent convolutional feature extraction units. Each convolutional feature extraction unit outputs thermal anomaly features at one scale, represented as follows:

[0229] ;

[0230] in, Indicates the first Thermal anomaly characteristics at various scales Indicates the thermal anomaly coding branch at the 1st Feature extraction functions at various scales; Indicates the first The local thermal anomaly map corresponding to the second flight is defined as follows:

[0231] ;

[0232] in, Indicates position The difference between the temperature response at a location and the average temperature response of its neighborhood is the local thermal anomaly map.

[0233] Standardized radar point cloud / elevation data Perform voxelization or road rasterization processing, and calculate local geometric features; for the first... Individual element units, with their center points represented as:

[0234] ;

[0235] in, Indicates the first A collection of point clouds within an individual element or road grid cell. Indicates voxel index, This represents the number of points in the set. Represents the set of elements. The three-dimensional coordinate vector of each radar point Represents a point index. , and These represent the horizontal and vertical coordinates of the point, respectively. Indicates the first The center point of the individual element;

[0236] No. The relative center point vector of an individual element is represented as:

[0237] ;

[0238] The standard deviation of point cloud is expressed as:

[0239] ;

[0240] The extreme value difference is expressed as:

[0241] ;

[0242] in, This represents the average offset vector of points within a voxel relative to the center point. This represents the standard deviation of point cloud coordinates. This represents the difference between the extreme values ​​of each coordinate component. and These represent taking the maximum and minimum values, respectively.

[0243] For the The local plane fitted to the neighborhood of a voxel is represented as:

[0244] ;

[0245] in, and Represents two horizontal coordinates within a local plane. Indicates the corresponding elevation. and This represents the plane slope coefficient obtained by fitting using the least squares method. Represents the plane intercept;

[0246] And calculate the local plane residual, using the formula:

[0247] ;

[0248] in, Indicates the first The average distance residual from each point within an individual element to the fitted plane; , , They represent the first Individual factors Inner Points The first horizontal coordinate component, the second horizontal coordinate component, and the elevation coordinate component in the local coordinate system; Indices representing points within a voxel;

[0249] Let the eigenvalues ​​of the local covariance matrix of the point cloud be... Calculate curvature characteristics:

[0250] ;

[0251] in, , and These are the eigenvalues ​​of the local covariance matrix of the point cloud, arranged in descending order. Indicates the first Normalized curvature of individual elements;

[0252] Elevation depression features are calculated based on radar point cloud / elevation data, for location elevation value at the location The neighborhood average elevation is expressed as:

[0253] ;

[0254] Location The elevation depression characteristics at the location are represented as follows:

[0255] ;

[0256] in, Indicates neighborhood location The elevation value at that location, Indicates position The elevation value at that location, Indicates position The neighborhood average elevation, This represents the normalized local elevation depression obtained from the difference between the neighborhood average elevation and the current location elevation;

[0257] For the Individual elements or road grid cells, the elevation depression characteristics are represented as follows:

[0258] ;

[0259] in, Indicates the first The average elevation depression characteristics of individual elements or road grid cells. Indicates the first Position of each grid cell within an individual element Find the arithmetic mean;

[0260] By concatenating the geometric features, the radar voxel-enhanced geometric features are obtained, represented as follows:

[0261] ;

[0262] Enhance the geometric features of each voxel Input radar geometry coding branch: First, the enhanced geometric features of each voxel are mapped to a unified feature dimension through a multilayer perceptron mapping unit. Then, the features of spatially adjacent voxels are aggregated by a neighborhood feature aggregation unit, either by maximum value aggregation or average value aggregation. Finally, radar geometric features at different scales are formed step by step through a downsampling unit, as shown below:

[0263] ;

[0264] in, Indicates the radar geometry coding branch in the first... Feature extraction functions at various scales The radar geometric features at the corresponding scale are represented; the multilayer perceptron mapping unit includes a fully connected layer and a nonlinear activation layer, and the neighborhood feature aggregation unit adopts maximum value aggregation or average value aggregation.

[0265] Step 4: Based on the cross-modal attention fusion mechanism guided by modal reliability, the three modal features are weighted and fused to obtain the fused features.

[0266] Furthermore, step 4 includes:

[0267] Map the three modalities of features to a unified feature space:

[0268] ;

[0269] ;

[0270] ;

[0271] in, Representing modes In the Input features at various scales; when The updated visible light texture features are used. ; , and Representing modes The query matrix, key matrix, and value matrix. , and For the corresponding learnable mapping matrix; when hour, Using thermal anomaly characteristics ;when hour, Using radar geometric features ;

[0272] Adjust the normalized reliability map to the corresponding scale, and align it with the first... The modal feature space has a consistent size across all scales, represented as:

[0273] ;

[0274] in, Representing modes In the Normalized reliability plots at various scales This indicates an interpolation scaling operation. This indicates the spatial dimensions of the corresponding scale feature map;

[0275] The normalized reliability map is incorporated as an attention bias into the cross-modal attention calculation:

[0276] ;

[0277] in, Indicates the first Modal at each scale For modes Attention contribution weight, This is the reliability adjustment factor. To query the feature dimensions of the key matrix and the key matrix, Represents the normalized exponential function, Represents the natural logarithm function;

[0278] Based on cross-modal attention weights, information from different modalities is fused to obtain the modality. Enhanced features , is represented as:

[0279] ;

[0280] Modality-level fusion weights are generated based on the global context and modality reliability, and are expressed as follows:

[0281] ;

[0282] in, Representing modes In the Fusion weights at various scales This indicates normalization in the modal dimension. This represents a multilayer perceptron. Indicates global average pooling. Indicates spatial average. This represents the reliability weighting adjustment coefficient;

[0283] The three modality enhancement features are weighted and summed to obtain the fused feature, which is represented as:

[0284] ;

[0285] By aggregating the fused features from multiple scales, a fused feature at a unified scale is obtained, represented as:

[0286] ;

[0287] in, This represents the final fusion features under a unified scale. This represents a multi-scale feature convergence operation. Indicates the total number of feature scales. Indicates the first Fusion features at various scales.

[0288] Step 5: Construct a geometry-thermal anomaly topology graph based on the fusion features, and update it through message passing and risk propagation via graph reasoning mechanism.

[0289] Furthermore, step 5 includes:

[0290] Final fusion features The topology graph node set is obtained by dividing the area into several nodes according to the road grid or candidate risk area. ,in, Indicates the first 1 node Indicates the node index. Indicates the total number of nodes;

[0291] No. The node characteristics of each node are represented as follows:

[0292] ;

[0293] in, Indicates the first The feature vector of each node This indicates a region pooling operation. Indicates the first Each road grid area or candidate risk area , and These represent the normalized reliability of the three modes, respectively.

[0294] Calculate the first The formulas for the spatial center location, average elevation depression characteristics, and average thermal anomaly characteristics of the region corresponding to each node are as follows:

[0295] ;

[0296] ;

[0297] ;

[0298] in, Indicates the first Each node corresponds to the spatial center of the region. This indicates the number of locations within the corresponding area. This indicates the degree of average elevation depression in the corresponding area. Indicates the first The average thermal anomaly intensity of the corresponding region during this flight. Indicates Calculate the arithmetic mean for the statistical domain;

[0299] Any two nodes and The edge weights between components are determined by spatial distance, semantic similarity, thermal anomaly consistency, elevation depression continuity, and crack connectivity, and are expressed as follows:

[0300] ;

[0301] ;

[0302] ;

[0303] ;

[0304] in, Indicates spatial proximity. Indicates the similarity of node features. Indicates the continuity of elevation depressions. Indicates consistency of thermal anomalies; The second norm of a vector, also known as the Euclidean norm. Represents the square of the Euclidean norm; Represents cosine similarity. , and A positive scaling parameter for measuring correspondence; and They represent the first The node and the first One node; , , and They represent the first The spatial center location of each node, node characteristics, average elevation depression degree, and the first Average thermal anomaly intensity during the second flight;

[0305] The crack skeleton map is generated from the edge enhancement features and is represented as follows:

[0306] ;

[0307] in, This represents a crack skeleton diagram. This indicates the skeletonization operation. This indicates a threshold segmentation operation; a crack skeleton map is generated to eliminate crack width differences and preserve centerline connectivity, providing a topological basis for subsequent calculations of the shortest crack path and connectivity between nodes;

[0308] Node calculation based on crack skeleton diagram With nodes The connectivity of the cracks between them is represented as:

[0309] ;

[0310] in, Represents a node With nodes The connectivity of the cracks between them This indicates that the two are in the crack skeleton diagram. The shortest path distance on the [path]. This represents the scale parameter indicating the connectivity of the crack;

[0311] Generate topological adjacency edge weights to control the message passing strength between adjacent nodes in the graph neural network, enabling spatial, feature, elevation, thermal anomaly, and crack connectivity evidence to participate in risk propagation, represented as:

[0312] ;

[0313] in, Represents a node To the node The topological adjacency edge weights after Softmax normalization based on the neighbor dimension. , , , and The non-negative weight parameters are for the five types of relationships. Represents a node Graph neighbor set, Index the neighboring nodes;

[0314] Introduce self-connects and normalize the adjacency matrix:

[0315] ;

[0316] in, This represents the normalized adjacency matrix obtained after adding self-joins. Indicates by The topological adjacency matrix that has been constructed but not yet normalized for symmetry. Representation matrix The degree matrix, express An identity matrix of order 1;

[0317] Node features are updated using a graph neural network, with the following formula:

[0318] ;

[0319] in, and They represent the first The layer's input and output node feature matrices, and the initial node feature matrix. Based on the characteristics of each node composition, Indicates the first The learnable parameter matrix of the layer, This represents the layer index of a graph neural network. (The last part of the equation is missing.) Convolution calculation with graph This represents the feature matrix of the same input node. It is element-wise added to the graph convolution result to form a residual connection, used to preserve the original node features and mitigate feature degradation caused by increasing the number of network layers. In this embodiment, the first... Layer graph convolution results and They have the same feature dimension;

[0320] Mapping the node features obtained after graph inference back to the road grid or pixel space yields graph-enhanced risk features, represented as:

[0321] ,

[0322] in, This represents the total number of layers in a graph neural network. This represents the inverse pooling operation that backfills node features into the corresponding road grid or pixel region. Indicates position The graph at that location enhances the risk characteristics.

[0323] Step 6: Decode and output the graph reasoning results to generate a pixel-level risk probability map, a collapsed hotspot mask, and a target area ranking result.

[0324] Furthermore, step 6 includes:

[0325] A pixel-level risk probability map is generated based on graph-enhanced risk features, represented as follows:

[0326] ;

[0327] in, Indicates position The probability of collapse at the location Indicates the risk decoding weight. Indicates risk decoding bias;

[0328] Generate a collapse hotspot mask based on the risk probability map; when hour, , indicating position It belongs to the collapse hotspot area; when hour, , indicating position It is not located in a collapse hotspot area; among them, Indicates position Binary hotspot markers at the location, Indicates the risk threshold;

[0329] Connectivity analysis was performed on the collapse hotspot areas to obtain the target region set. ,in, Represents the target region set. Indicates the first A connected target region, Indicates the target region index. Indicates the total number of target areas;

[0330] And calculate the first Risk score for each target area:

[0331] ,

[0332] in, Indicates the first Risk score for each target area Indicates the target area Internal locations The arithmetic mean of the risk probabilities. Indicate target area The area;

[0333] Based on risk score All target areas are sorted to obtain target area sorting results, and the risk priority of different suspected collapse areas is characterized by the target area sorting results.

[0334] To verify the supporting role of elevation / geometric anomaly modeling, RGB-thermal infrared fusion, and topological reasoning in road hazard identification in this invention, this embodiment uses two sets of publicly available datasets for surrogate experiments. Since currently available datasets rarely contain RGB, thermal infrared, radar / elevation data and pixel-level collapse annotations for urban road collapse scenarios, the verification process is split into elevation / geometric surrogate experiments and RGB-thermal infrared crack identification experiments.

[0335] The first set of experiments used the PothRGBD dataset (Yurdakul M, Tasdemir Ş. An enhanced yolov8 model for real-time and accurate pothole detection and measurement[J]. arXiv preprint arXiv:2505.04207, 2025.). The original publicly available size of this dataset is 1000 pairs of RGB-depth paired images. In this example, after data integrity checks, 998 pairs of paired samples with valid depth maps and pixel-level annotations were selected, with 698, 149, and 151 pairs for training, validation, and test sets, respectively. In the experiments, the depth map was used as a proxy input for the radar / elevation grid, and local elevation depression, positive and negative residuals, curvature response, and validity reliability channels were constructed to verify the effectiveness of elevation geometric anomaly modeling and topological risk propagation in this invention. The experimental results are shown in Table 1 and... Figure 3 As shown.

[0336] The comparison methods in Table 1 are set as follows: Geometry heuristic is a non-learned geometric baseline constructed by this example, which generates segmentation results based on local elevation depressions, planar residuals, and curvature responses after thresholding; Geometry heuristic + topology adds neighborhood topological smoothing to the output of this geometric baseline; Depth U-Net is a single-channel depth segmentation baseline constructed by this example based on the U-Net encoder-decoder structure (Ronneberger O, Fischer P, Brox T. U-net: Convolutional networks for biomedical imagesegmentation[C] / / International Conference on Medical image computing and computer-assisted intervention. Cham: Springer international publishing,2015: 234-241.); Reliability geometry U-Net adds elevation depression, planar residual, curvature, and validity reliability channels to the above U-Net input; Reliability geometry U-Net + The topology (Reliability Geometry U-Net + Topology) further incorporates the topology risk propagation of this invention into the output of the Reliability Geometry U-Net, providing the corresponding implementation of the elevation / geometry branch of this invention on the proxy dataset. Apart from the deep U-Net infrastructure, the aforementioned geometric input combinations, reliability channels, and topology processing are all control or ablation configurations set according to the technical solution of this invention.

[0337] Tables 1 and 2 use the same pixel-level evaluation metrics: Precision represents the proportion of pixels predicted as potential hazards to actual hazard pixels; Recall represents the proportion of pixels correctly detected among actual hazard pixels; Dice / F1 represents the harmonic mean of Precision and Recall, reflecting the degree of overlap between the predicted and actual regions; IoU represents the ratio of the intersection area to the union area of ​​the predicted and actual regions; Accuracy represents the proportion of correctly classified pixels; AUC represents the area under the receiver operating characteristic curve, used to measure the overall discriminative ability at different decision thresholds. Higher values ​​for these metrics indicate better recognition performance.

[0338] Table 1. Results of PothRGBD Elevation / Geometric Proxy Experiment

[0339]

[0340] Table 1 and Figure 3 The results show that the geometric heuristic method achieves a Dice / F1 ratio of 0.3986 and an IoU of 0.2489. After adding topological smoothing, the Dice / F1 ratio improves to 0.4021, and the IoU improves to 0.2516. Using a deep U-Net, the Dice / F1 ratio improves to 0.7101, and the IoU improves to 0.5505. Further introducing reliability geometric features results in a Dice / F1 ratio of 0.7532 and an IoU of 0.6041. Combining topological inference further improves the Dice / F1 ratio to 0.7546, the IoU to 0.6059, and the AUC to 0.9277. These results demonstrate that elevation depressions, geometric reliability, and topological continuity can effectively enhance the identification ability of local road depression areas.

[0341] The second set of experiments used the IRT-Crack Segmentation dataset (Liu F, Liu J, Wang L. Asphalt pavement crack detection based on convolutional neural network and infrared thermography[J]. IEEE Transactions on Intelligent Transportation Systems, 2022, 23(11): 22145-22155.), which contained 448 sets of visible light images, thermal infrared images, pre-fused images, and pixel-level crack annotation samples. The training set, validation set, and test set were 358, 44, and 46 sets, respectively. This experiment was used to verify the supporting role of the RGB texture coding branch, thermal infrared anomaly branch, and RGB-T fusion output in the present invention for the identification of road cracks and potential collapse precursors. The experimental results are shown in Table 2 and... Figure 4 As shown.

[0342] Table 2 uses the same U-Net backbone network (Yurdakul M, Tasdemir Ş. An enhanced yolov8 model for real-time and accurate pothole detection and measurement[J].arXiv preprint arXiv:2505.04207, 2025.) and evaluation metrics as Table 1, with only the input modalities and fusion strategies differing: RGB U-Net (Visible Light U-Net) in Table 2 uses only visible light images as input; Infrared U-Net (Thermal Infrared U-Net) uses only thermal infrared images as input; Pre-fused Image U-Net uses pre-fused images provided by the input dataset; RGB-T Early Fusion U-Net (RGB-Thermal Infrared Early Fusion U-Net) concatenates visible light and thermal infrared images by channel before inputting them into the network; Reliability-guided RGB-T U-Net (Reliability-guided RGB-Thermal Infrared U-Net) adds local reliability weights during the dual-modal fusion process; the last row in the table shows the corresponding implementation of the RGB-Thermal Infrared branch of this invention on the proxy dataset. The first four methods are all input modalities or fusion configurations set on the same U-Net backbone in this experiment, and do not correspond to the same-named algorithms in external papers; the latter two methods are ablation methods and complete corresponding implementations set according to the technical solution of this invention.

[0343] Table 2. Experimental Results of IRT-Crack RGB-Thermal Infrared Crack Identification

[0344]

[0345] Table 2 and Figure 4 The results show that RGB U-Net has a Dice / F1 ratio of 0.8120, a recall of 0.8114, and an AUC of 0.9912. RGB-T Early Fusion U-Net's recall is improved to 0.8268, and its AUC to 0.9923, indicating that thermal infrared modal and visible light texture information have complementary effects. The Reliability-guided RGB-T U-Net + topology achieves a recall of 0.8295, an improvement of 0.0181 compared to RGB U-Net, demonstrating that reliability guidance and topology propagation help reduce missed detections in road inspection tasks that emphasize the recall of suspected hazards.

[0346] Figure 5Key improvements corresponding to the core technology of this invention are summarized as follows: In the PothRGBD experiment, the reliability geometric U-Net improves the IoU by 0.0536 compared to the deep U-Net, and further improves it to 0.6059 after adding topology inference; In the IRT-Crack experiment, the RGB-T fusion model improves both Recall and AUC compared to the RGB single-modal model, and the topology inference model achieves the highest Recall. Figure 6 The qualitative results further show that the elevation / depth risk heat map can cover pits or local depressions, and the RGB-T model can complement the visible crack boundary and the thermal infrared response.

[0347] The combined results of the two sets of experiments demonstrate that this invention does not rely solely on single image appearance features, but rather employs a joint assessment based on elevation geometric anomalies, thermal infrared anomalies, visible light textures, and topological spatial relationships. The reliability-guided mechanism enhances the robustness of multimodal fusion, while the topological inference mechanism strengthens the continuity and recall capability of risk areas. Therefore, the experimental results provide evidence of the effectiveness of the method described in this invention in monitoring urban ground subsidence hazards.

Claims

1. A ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning, characterized in that, Includes the following steps: Step 1: Collect visible light images, thermal infrared images, and radar point cloud / elevation data of urban road areas, and preprocess the multimodal data to obtain standardized input data; Step 2: Perform quality assessment on the standardized input data, calculate modal quality terms, registration error terms, and environmental interference terms, and generate reliability maps for each modality; Step 3: Construct three feature coding branches with independent parameters, and input the standardized input data into the corresponding feature coding branches to obtain multi-scale visible light texture features, thermal anomaly features, and radar geometric features; Step 4: Based on the cross-modal attention fusion mechanism guided by modal reliability, the three-modal features are weighted and fused to obtain the fused features; Step 5: Construct a geometric-thermal anomaly topology map based on the fusion features, and update it through message passing and risk propagation via graph reasoning mechanism; Step 6: Decode and output the graph reasoning results to generate a pixel-level risk probability map, a collapsed hotspot mask, and a target area ranking result.

2. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 1, characterized in that, Step 1 includes: Using a drone platform equipped with a visible light camera, a thermal infrared camera, and a radar sensor, simultaneous aerial data collection was conducted on the urban road area to be monitored, obtaining the first... Multimodal data set corresponding to this flight ,in, Indicates the first Visible light images obtained during this flight. Indicates the first Thermal infrared images obtained during this flight. Indicates the first Radar point cloud / elevation data obtained during this flight , and These represent the visible light mode, the thermal infrared mode, and the lidar mode, respectively. This indicates the drone's pose, sensor intrinsic and extrinsic parameters, and georeferenced information. This represents the data acquisition timestamp, which is used to synchronize and match the three-modal data and associate data acquired from the same flight with the same flight sequence number. Indicates the flight sequence number; For a point in three-dimensional space In modal The projection position in the image is represented as: ； in, Represents a set of modal identifiers. Indicates modal identifier, , Represents the homogeneous coordinate vector of a point in three-dimensional space. , and Indicates three-dimensional coordinate components, superscript Indicates transpose. This indicates that the point is in mode The two-dimensional projection position in , and Representing modes The corresponding intrinsic parameter matrix, rotation matrix, and translation vector of the sensor. Representing modes The projection function; Using the visible light image coordinate system as the reference coordinate system, thermal infrared imagery and radar point cloud / elevation data are mapped to a unified road grid coordinate system, as follows: ； ； in, This represents the spatial transformation function from the thermal infrared image coordinate system to the visible light reference coordinate system. This represents the spatial transformation function from the radar coordinate system to the visible light reference coordinate system. and These represent the mapped thermal infrared image and radar point cloud / elevation data, respectively. The mapped multimodal data is then subjected to denoising, normalization, scale unification, and image segmentation to obtain a standardized input data set. ;in, Represents a standardized set of input data. , and These represent visible light images, thermal infrared images, and radar point cloud / elevation data after registration, denoising, normalization, and scale unification, respectively.

3. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 2, characterized in that, Step 2 includes: For the visible light mode, its quality term is jointly determined by image sharpness, local texture information, and exposure anomalies, and is expressed as: ； in, Indicates the visible light mode at position The mass value at that location, Indicates pixel position or road grid position. This indicates that the input is normalized to... Functions on intervals Represents local variance. Represents the Laplace operator. Indicates the location of standardized visible light images Pixel value at that location, Indicates overexposure or underexposure interference; For the thermal infrared mode, its mass term is determined by the contrast of local thermal anomalies, and is expressed as: ； in, Indicates the thermal infrared mode at position The mass value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates the location of standardized thermal infrared images Temperature response value at that location, Indicates position The neighborhood, Indicates the number of pixels in the neighborhood. Representing the neighborhood Pixel position index within; For radar modes, the mass term is determined by both point cloud density and elevation stability, and is expressed as: ； in, Indicates the radar mode at position The mass value at that location, Indicates position The number of point clouds within a corresponding raster or voxel. This represents the maximum number of point clouds within a local statistical range. This represents the standard deviation of point cloud elevation within the location's neighborhood. This represents the positive scaling parameter that controls the rate of decay of elevation dispersion. Represents the natural exponential function; modality In position The registration error term at the location is defined as: ； in, Representing modes In position Registration error at the location, Representing modes Center and reference position The corresponding position Indicates modality Spatial transformation function from the coordinate system to the visible light reference coordinate system. Represents the Euclidean norm; modality The environmental disturbance term is represented as: ； in, Representing modes In position The amount of environmental disturbance at the location , , and These represent shadows, reflections, sparse point clouds, and noise interference, respectively. , , and These are the non-negative weight coefficients for the corresponding interference terms; Calculate the modal quality term, registration error term, and environmental disturbance term. In position The reliability at the location is given by the formula: ； in, Representing modes In position Reliability before normalization at the location, , and These represent the modal quality term, registration error term, and environmental disturbance term, respectively. This represents the Sigmoid function. , and It is a non-negative scaling factor. These are bias parameters; The reliability of each mode is normalized to obtain a normalized reliability map, which is represented as follows: ； in, Represents the normalized mode reliability, Indicates the modal summation index. This represents a very small positive number that prevents the denominator from being zero.

4. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 3, characterized in that, Step 3 includes: Standardize visible light images The input is an RGB texture encoding branch; the RGB texture encoding branch includes multiple levels of convolutional feature extraction units connected in sequence, with downsampling units set between adjacent levels of convolutional feature extraction units, and each level of convolutional feature extraction unit outputs a visible light feature at one scale, represented as: ； in, Indicates the first Visible light characteristics at various scales This represents a multi-scale feature hierarchy index. This indicates that the RGB texture encoding branch is in the first... Feature extraction functions at various scales; Each convolutional feature extraction unit includes a convolutional layer, a normalization layer, and a nonlinear activation layer. The downsampling unit employs stride convolution or pooling operations. To enhance the information on road crack boundaries and surface morphology, an edge enhancement feature is introduced, represented as: ； in, Indicates the first Edge enhancement features at various scales Represents a non-linear activation function. This represents the convolution operation. This indicates the edge extraction operator. This indicates feature splicing based on channels; Edge enhancement features are fused with visible light multi-scale features to obtain updated visible light texture features. , represented as: ； Standardized thermal infrared images With local thermal anomaly map The thermal anomaly coding branch is input after concatenation by channel; the thermal anomaly coding branch includes a series of convolutional feature extraction units connected in sequence, with a downsampling unit set between adjacent convolutional feature extraction units, and each convolutional feature extraction unit outputs thermal anomaly features at one scale, represented as follows: ； in, Indicates the first Thermal anomaly characteristics at various scales Indicates the thermal anomaly coding branch at the 1st Feature extraction functions at various scales Indicates the first The local thermal anomaly map corresponding to the second flight, wherein the local thermal anomaly map is defined as: ； in, Indicates position The difference between the temperature response at a location and the average temperature response of its neighborhood is the local thermal anomaly map. Standardized radar point cloud / elevation data Perform voxelization or road rasterization processing, and calculate local geometric features; for the first... Individual element units, with their center points represented as: ； in, Indicates the first The center point of the individual element, Indicates the first A collection of point clouds within an individual element or road grid cell. Indicates voxel index, This represents the number of points in the set. Represents the set of elements. The three-dimensional coordinate vector of each radar point Represents a point index. , and These represent the horizontal and vertical coordinates of the point, respectively. No. The relative center point vector of an individual element is represented as: ； The standard deviation of point cloud is expressed as: ； The extreme value difference is expressed as: ； in, This represents the average offset vector of points within a voxel relative to the center point. This represents the standard deviation of point cloud coordinates. This represents the difference between the extreme values ​​of each coordinate component. and These represent taking the maximum and minimum values, respectively. For the The local plane fitted to the neighborhood of a voxel is represented as: ； in, and Represents two horizontal coordinates within a local plane. Indicates the corresponding elevation. and This represents the plane slope coefficient obtained by fitting using the least squares method. Represents the plane intercept; And calculate the local plane residual, using the formula: ； in, Indicates the first The average distance residual from each point within an individual element to the fitted plane. , , They represent the first Individual factors Inner Points The first horizontal coordinate component, the second horizontal coordinate component, and the elevation coordinate component in the local coordinate system; Let the eigenvalues ​​of the local covariance matrix of the point cloud be... Calculate curvature characteristics: ； in, , and These are the eigenvalues ​​of the local covariance matrix of the point cloud, arranged in descending order. Indicates the first Normalized curvature of individual elements; Elevation depression features are calculated based on radar point cloud / elevation data, for location Elevation value at the location The neighborhood average elevation is expressed as: ； Location The elevation depression characteristics at the location are represented as follows: ； in, Indicates neighborhood location The elevation value at that location, Indicates position The elevation value at that location, Indicates position The neighborhood average elevation, This represents the normalized local elevation depression obtained from the difference between the neighborhood average elevation and the current location elevation; For the Individual elements or road grid cells, the elevation depression characteristics are represented as follows: ； in, Indicates the first The average elevation depression characteristics of individual elements or road grid cells. Indicates the first Position of each grid cell within an individual element Find the arithmetic mean; By concatenating the geometric features, the radar voxel-enhanced geometric features are obtained, represented as follows: ； Enhance the geometric features of each voxel Input radar geometry coding branch: First, the enhanced geometric features of each voxel are mapped to a unified feature dimension through a multilayer perceptron mapping unit. Then, the features of spatially adjacent voxels are aggregated by a neighborhood feature aggregation unit, either by maximum value aggregation or average value aggregation. Finally, radar geometric features at different scales are formed step by step through a downsampling unit, as shown below: ； in, Indicates the radar geometry coding branch in the first... Feature extraction functions at various scales The radar geometric features are represented at the corresponding scale; the multilayer perceptron mapping unit includes a fully connected layer and a nonlinear activation layer, and the neighborhood feature aggregation unit adopts maximum value aggregation or average value aggregation.

5. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 4, characterized in that, Step 4 includes: Map the three modalities of features to a unified feature space: ； ； ； in, Representing modes In the Input features at various scales; when The updated visible light texture features are used. ; , and Representing modes The query matrix, key matrix, and value matrix. , and For the corresponding learnable mapping matrix; when hour, Using thermal anomaly characteristics ;when hour, Using radar geometric features ; Adjust the normalized reliability map to the corresponding scale, and align it with the first... The modal feature space has a consistent size across all scales, represented as: ； in, Representing modes In the Normalized reliability plots at various scales This indicates an interpolation scaling operation. This indicates the spatial dimensions of the corresponding scale feature map; The normalized reliability map is incorporated as an attention bias into the cross-modal attention calculation: ； in, Indicates the first Modal at each scale For modes Attention contribution weight, This is the reliability adjustment factor. To query the feature dimensions of the key matrix and the key matrix, Represents the normalized exponential function, Represents the natural logarithm function; Based on cross-modal attention weights, information from different modalities is fused to obtain the modality. Enhanced features , represented as: ； Modality-level fusion weights are generated based on the global context and modality reliability, and are expressed as follows: ； in, Representing modes In the Fusion weights at various scales This indicates normalization in the modal dimension. This represents a multilayer perceptron. Indicates global average pooling. Indicates spatial average. This represents the reliability weighting adjustment coefficient; The three modality enhancement features are weighted and summed to obtain the fused feature, which is represented as: ； By aggregating the fused features from multiple scales, a fused feature at a unified scale is obtained, represented as: ； in, This indicates a multi-scale feature convergence operation. Indicates the total number of feature scales. Indicates the first Fusion features at various scales This represents the final fusion feature under a unified scale.

6. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 5, characterized in that, Step 5 includes: Final fusion features The topology graph node set is obtained by dividing the area into several nodes according to the road grid or candidate risk area. ,in, Indicates the first 1 node Indicates the node index. Indicates the total number of nodes; No. The node characteristics of each node are represented as follows: ； in, Indicates the first The feature vector of each node This indicates a region pooling operation. Indicates the first Each road grid area or candidate risk area , and These represent the normalized reliability of the three modes, respectively. Calculate the first The formulas for the spatial center location, average elevation depression characteristics, and average thermal anomaly characteristics of the region corresponding to each node are as follows: ； ； ； in, Indicates the first Each node corresponds to the spatial center of the region. This indicates the number of locations within the corresponding area. This indicates the degree of average elevation depression in the corresponding area. Indicates the first The average thermal anomaly intensity of the corresponding region during this flight. Indicated by Calculate the arithmetic mean for the statistical domain; Any two nodes and The edge weights between components are determined by spatial distance, semantic similarity, thermal anomaly consistency, elevation depression continuity, and crack connectivity, and are expressed as follows: ； ； ； ； in, Indicates spatial proximity. Indicates the similarity of node features. Indicates the continuity of elevation depressions. Indicates consistency of thermal anomalies; The second norm of a vector, also known as the Euclidean norm. Represents the square of the Euclidean norm; Represents cosine similarity. , and A positive scaling parameter for measuring correspondence; and They represent the first The node and the first One node; , , and They represent the first The spatial center location of each node, node characteristics, average elevation depression degree, and the first Average thermal anomaly intensity during the second flight; The crack skeleton map is generated from the edge enhancement features and is represented as follows: ； in, This represents a crack skeleton diagram. This indicates a threshold segmentation operation. This indicates a skeletonization operation; Node calculation based on crack skeleton diagram With nodes The connectivity of the cracks between them is represented as: ； in, Represents a node With nodes The connectivity of the cracks between them This indicates that the two are in the crack skeleton diagram. The shortest path distance on the road, This represents the scale parameter indicating the connectivity of the crack; The topological adjacency edge weights are generated as follows: ； in, Represents a node To the node The topological adjacency edge weights after Softmax normalization based on the neighbor dimension. , , , and The non-negative weight parameters are for the five types of relationships. Represents a node Graph neighbor set, Index the neighboring nodes; Introduce self-connects and normalize the adjacency matrix: ； in, This represents the normalized adjacency matrix obtained after adding self-joins. Indicates by The topological adjacency matrix that has been constructed but not yet normalized for symmetry. Representation matrix The degree matrix, express An identity matrix of order 1; Node features are updated using a graph neural network, with the following update formula: ； in, and They represent the first The layer's input and output node feature matrices, and the initial node feature matrix. Based on the characteristics of each node composition, Indicates the first The learnable parameter matrix of the layer, Indicates the layer index of a graph neural network; Mapping the node features obtained after graph inference back to the road grid or pixel space yields graph-enhanced risk features, represented as: ； in, Indicates position The graph at that location enhances risk characteristics. This represents the total number of layers in a graph neural network. This represents the inverse pooling operation that backfills node features into the corresponding road grid or pixel region.

7. The ground subsidence monitoring method based on space-based multimodal reliability fusion and graph reasoning according to claim 6, characterized in that, Step 6 includes: A pixel-level risk probability map is generated based on graph-enhanced risk features, represented as follows: ； in, Indicates position The probability of collapse at the location Indicates the risk decoding weight. Indicates risk decoding bias; Generate a collapse hotspot mask based on the risk probability map: when hour, , indicating position It belongs to the collapse hotspot area; when hour, , indicating position It is not located in a collapse hotspot area; among them, Indicates position Binary hotspot markers at the location, Indicates the risk threshold; Connectivity analysis was performed on the collapse hotspot areas to obtain the target region set. ,in, Represents the target region set. Indicates the first A connected target region, Indicates the target region index. Indicates the total number of target areas; And calculate the first Risk score for each target area: ； in, Indicates the first Risk score for each target area Indicates the target area Internal locations The arithmetic mean of the risk probabilities. Indicate target area The area; Based on the risk score All target areas are sorted to obtain target area sorting results, and the risk priority of different suspected collapse areas is characterized by the target area sorting results.

Images