Unmanned aerial vehicle-based insar ground deformation information detection method and system

Abstract

The present application relates to the technical field of topographic survey, in particular to an InSAR ground deformation information detection method and system based on a UAV, comprising the following steps: S1: collecting point cloud data of a mining area by a laser radar system carried by a UAV, constructing a continuous terrain surface, and calculating a terrain gradient vector based on the surface to extract a main terrain gradient direction angle; S2: constructing a physical constraint field quantifying the terrain mutation risk and the interaction effect of mining and terrain based on the terrain gradient vector, the rock mass integrity coefficient of the mining area, and the mining direction angle. The present application generates a physical field quantifying the fault risk and type by fusing the UAV laser radar terrain data, the rock mass integrity coefficient, and the artificial mining direction; geological information is changed from post-explanation to prior constraint, providing accurate physical basis for subsequent phase processing and preventing misjudgment of the deformation direction, such as misjudging subsidence as uplift.

Description

Technical Field

[0001] This invention relates to the field of topographic surveying technology, and more specifically to an InSAR surface deformation information detection method and system based on unmanned aerial vehicles (UAVs). Background Technology

[0002] Surface deformation refers to the slow or sudden displacement, subsidence, or uplift of the Earth's surface caused by various natural factors (such as earthquakes, volcanic activity, glacial melting, and changes in groundwater) and human activities (such as mineral resource extraction, excessive urban groundwater extraction, large-scale engineering construction, and oil and gas extraction). Surface deformation not only directly threatens people's lives and property but also poses a serious challenge to the stability of major infrastructure (such as dams, bridges, railways, and high-rise buildings) and regional sustainable development. Therefore, high-precision, large-scale, and timely monitoring of surface deformation is a key technological requirement in fields such as disaster early warning, resource management, and engineering safety assessment.

[0003] As a significant breakthrough in the field of remote sensing, UAV-borne InSAR surface deformation monitoring technology achieves high-precision dynamic monitoring of deformation phenomena such as surface subsidence and landslides in mining areas by integrating a synthetic aperture radar interferometry system into a UAV platform. This technology analyzes millimeter-level surface displacement by receiving phase differences in radar echo signals, with its core objective being to provide real-time early warning support for mine safety production and geological disaster prevention. The existing standardized process includes steps such as periodic aerial surveying, interferogram generation, phase unwrapping, and deformation calculation, and has been successfully applied in the field of mine stability assessment.

[0004] However, in mining scenarios, when blasting in open-pit mining or the advancement of underground working faces induces step-like fault collapses, the surface undergoes discontinuous, abrupt deformation within a very short time. This abrupt change causes the displacement difference between adjacent observation points to exceed the radar wavelength's carrying capacity, triggering periodic ambiguity during phase unwrapping. The phase continuity assumption upon which existing algorithms rely completely fails in this scenario, causing the system to misjudge large-scale subsidence as reverse uplift and failing to accurately locate fault abrupt boundary. Consequently, misjudgments of deformation direction and magnitude can mask the true collapse risk, leading to the failure of the disaster early warning system. Furthermore, fault boundary identification errors prevent subsequent support engineering from accurately locating dangerous areas. These two deficiencies directly cause delays and inaccuracies in disaster prevention and control measures, significantly increasing the probability of mining accidents. This technological bottleneck has become a key obstacle restricting safe and efficient mining operations, necessitating the development of dedicated monitoring methods for discontinuous deformation scenarios.

[0005] Therefore, a method and system for detecting InSAR surface deformation information based on unmanned aerial vehicles (UAVs) are proposed to address the aforementioned problems. Summary of the Invention

[0006] Technical problems to be solved:

[0007] To address the aforementioned shortcomings of existing technologies, this invention provides an InSAR surface deformation information detection method and system based on unmanned aerial vehicles (UAVs). This effectively solves the problems of phase unwrapping failure and deformation misjudgment caused by severe deformation in mines, which are common in traditional InSAR technology. It achieves the goal of high-precision monitoring and accurate early warning of discontinuous deformations such as faults and collapses.

[0008] Technical solution:

[0009] To achieve the above objectives, the present invention provides the following technical solution:

[0010] This invention provides a method for detecting InSAR surface deformation information based on unmanned aerial vehicles (UAVs), comprising the following steps:

[0011] S1: Collect point cloud data of the mining area using the lidar system carried by the drone, construct a continuous terrain surface, calculate the terrain gradient vector based on the surface, and extract the main gradient direction angle of the terrain.

[0012] S2: Based on the topographic gradient vector, the rock mass integrity coefficient of the mining area, and the mining direction angle, a physical constraint field is constructed to quantify the risk of terrain abrupt change and the interaction effect between mining and terrain;

[0013] S3: Based on the initial interferometric phase diagram and the physical constraint field, calculate the compensation value and compensate the initial interferometric phase diagram to obtain the reconstructed continuous phase field;

[0014] S4: Based on the reconstructed continuous phase field and the physical constraint field, the phase is unwrapped using a direction-adaptive unwrapping algorithm to obtain a continuous phase field with a unified global reference.

[0015] S5: Transform the global benchmark unified continuous phase field into a three-dimensional spatiotemporal evolution deformation field, and perform accuracy verification and disaster risk assessment.

[0016] Furthermore, step S1 also includes:

[0017] S101: Obtain the original 3D point cloud dataset of the mining area;

[0018] S102: Preprocess the original 3D point cloud dataset to obtain a clean surface point cloud dataset;

[0019] S103: Construct a continuous terrain surface based on the aforementioned pure surface point cluster. ;

[0020] S104: Based on the continuous terrain surface, calculate the terrain gradient vector at each location, i.e.:

[0021]

[0022] In the formula, S is the topographic gradient vector; S is the elevation field, i.e., the continuous topographic surface. ; Let be the partial derivative of the elevation field in the x-direction; Let be the partial derivative of the elevation field in the y-direction;

[0023] S105: Based on the terrain gradient vector, calculate the principal gradient direction angle of the terrain, i.e.:

[0024]

[0025] In the formula, The direction angle of the principal gradient of the terrain; It is the arctangent function.

[0026] Further, step S2 includes:

[0027] S201: Obtain the distribution map of rock mass integrity coefficient and dynamic matrix of mining direction in the mining area, and perform data fusion;

[0028] S202: Based on the fused data, construct a deformation sensitivity function, namely:

[0029]

[0030] In the formula, Deformation sensitivity; The normalized rock series number;

[0031] S203: Based on the dynamic matrix of the mining direction and the principal gradient direction angle of the terrain, calculate the directional synergy factor, i.e.:

[0032]

[0033] In the formula, For directional coordination factors; This indicates the direction of mining operations;

[0034] S204: Based on the deformation sensitivity and the orientation coordination factor, calculate the physical constraint field to indicate the fault type and activity level, i.e.:

[0035]

[0036] In the formula, It is a physical constraint field.

[0037] Further, step S3 includes:

[0038] S301: Generate initial interferometric phase map ;

[0039] S302: Calculate the phase gradient magnitude of the initial interferometric phase pattern, i.e.:

[0040]

[0041] In the formula, The magnitude of the phase gradient;

[0042] S303: Identify the phase transition region based on the phase gradient magnitude, that is:

[0043]

[0044] In the formula, This is a distribution diagram of the phase jump domain;

[0045] S304: Based on mining plan data, calculate the estimated surface displacement, i.e.:

[0046]

[0047] In the formula, Estimated value of vertical displacement of the Earth's surface; The radar incident angle; This refers to the satellite's heading angle; This represents a horizontal displacement in the east-west direction. This represents a north-south horizontal displacement. For vertical pre-displacement field;

[0048] S305: Convert the estimated surface displacement into the theoretical number of compensation cycles, i.e.:

[0049]

[0050] In the formula, This represents the theoretical number of compensation cycles. It is a rounding function; The radar wavelength;

[0051] S306: Based on the physical constraint field, the theoretical compensation cycle number is constrained to obtain the final compensation value, i.e.:

[0052]

[0053] In the formula, k is the compensation value;

[0054] S307: The initial interferometric phase map is compensated using the final compensation value to obtain the reconstructed continuous phase field, i.e.:

[0055]

[0056] In the formula, This is the reconstructed continuous phase field.

[0057] Furthermore, in step S304, the estimated surface displacement is calculated as follows:

[0058] S3041: Calculate the radius of influence of mining based on the depth of the mining area and the fracture angle of the rock strata, i.e.:

[0059]

[0060] In the formula, The radius of influence of mining activities; The depth of the mining area; Rock strata fracture angle.

[0061] S3042: Calculate the vertical pre-displacement field based on the mining boundary and mining thickness of the mining area to predict the vertical subsidence of the ground surface, i.e.:

[0062]

[0063] In the formula, The vertical pre-displacement field is represented by q, where q is the subsidence coefficient. For the thickness of the ore; It is the error function; This refers to the corner point of the mining boundary; The radius of influence of the mining operation is the x-axis radius component; The radius of influence of the sampling is the radius of influence component on the y-axis.

[0064] S3043: Calculate the east-west and north-south horizontal displacements based on the vertical pre-displacement field and the radius of influence of mining, i.e.:

[0065]

[0066] In the formula, b is the horizontal movement coefficient based on experimental calibration.

[0067] Further, step S4 includes:

[0068] S401: Calculate the phase derivative field of the reconstructed continuous phase field, i.e.:

[0069]

[0070] In the formula, The phase gradient is in the x-axis direction; The phase gradient is in the y-axis direction; For phase gradient field; Let be the reconstructed continuous phase field value at position (x, y);

[0071] S402: Based on the phase derivative field and the terrain principal gradient direction angle, calculate the path safety weight coefficient to determine the optimal integration path and merge them into the optimal path set, i.e.:

[0072]

[0073] In the formula, For terrain-adaptive path safety weighting coefficients;

[0074] S403: Based on the optimal path set, calculate the unbiased phase difference between adjacent pixels, i.e.:

[0075]

[0076] In the formula, The unambiguous phase difference between the i-th and j-th points in the optimal path set; The i-th pixel in the optimal path set; The j-th pixel in the optimal path set;

[0077] S404: Perform phase integration along the optimal integration path to obtain the preliminary unwrapped phase field, i.e.:

[0078]

[0079] In the formula, The absolute phase value after untangling at the target point; The reference point phase value for the reconstructed continuous phase field; For directional accumulation operators along the optimal path set;

[0080] S405: Based on the aforementioned physical constraint field and phase jump boundary, construct the fault barrier function, namely:

[0081]

[0082] In the formula, This represents the intensity field of geological fault activity. The Euclidean distance from the current point (x, y) to the fault boundary; For the boundary of the jump domain; For physical constraint fields;

[0083] S406: Based on the fault barrier function, the entire domain is divided into multiple continuous sub-regions for independent unwrapping to obtain the multi-sub-region unwrapped phase field;

[0084] S407: Perform global continuity correction on the unwrapped phase field of the multi-sub-region to obtain a global reference-unified continuous phase field, i.e.:

[0085]

[0086] In the formula, A continuous phase field with a unified global reference. Unwrap the phase field for the multi-sub-region; To determine the optimal cycle offset.

[0087] Furthermore, the boundary of the transition domain is calculated as follows:

[0088] S4051: The curl field for calculating the phase gradient based on the phase derivative field, i.e.:

[0089]

[0090] In the formula, The curl field represents the phase gradient; The starting gradient is in the x-axis direction; The gradient is the upper right gradient along the y-axis. For gradient compensation in the x-axis direction; Gradient compensation for the starting point in the y-axis direction;

[0091] S4052: Set the boundary of the jump domain based on the curl field of the phase gradient, that is:

[0092]

[0093] In the formula, This is the boundary of the jump domain.

[0094] Further, step S5 includes:

[0095] S501: The continuous phase field of the global reference is transformed into a radar line-of-sight displacement field, that is:

[0096]

[0097] In the formula, This is the displacement vector along the radar line of sight; The gradient operator for a continuous phase field with a unified global reference. The radar wavelength;

[0098] S502: Decompose the radar line-of-sight displacement field into three-dimensional displacement components in the vertical, east-west, and north-south directions to construct a three-dimensional displacement field, namely:

[0099]

[0100] In the formula, This represents the vertical displacement component; This represents the east-west displacement component; This represents the displacement component in the north-south direction; It is the inverse of the satellite's geometric projection matrix; It is a three-dimensional displacement field;

[0101] S503: Construct a time-series three-dimensional displacement field set based on multi-period monitoring data;

[0102] S504: Calculate the displacement rate field and cumulative displacement field based on the aforementioned three-dimensional displacement field set, i.e.:

[0103]

[0104] In the formula, For displacement velocity field; and Represents the three-dimensional displacement field set in adjacent time phases and Displacement observations;

[0105]

[0106] In the formula, To accumulate the displacement field, the displacement from the initial time was quantized. The cumulative displacement up to the current time t;

[0107] S505: The cumulative displacement field is corrected using data from ground monitoring stations to obtain the final displacement field, i.e.:

[0108]

[0109] In the formula, The magnitude of the residual vector; This is the measured displacement field; This is the corrected final displacement field;

[0110] S506: Based on the final displacement field and the physical constraint field, calculate the fault differential motion and activation risk index, i.e.:

[0111]

[0112] In the formula, This represents the vertical displacement on both sides of the fault; This is a sign function used to determine normal / reverse faults;

[0113]

[0114] In the formula, This is the fault activation risk index.

[0115] The UAV-based InSAR surface deformation information detection system includes:

[0116] The terrain gradient extraction module is used to collect and process point cloud data of the mining area through the lidar system carried by the UAV, and to calculate the terrain gradient vector and the terrain principal gradient direction angle.

[0117] The geomechanical constraint field construction module is used to construct physical constraint fields based on topographic data, rock mass data, and mining data.

[0118] The phase jump compensation module is used to perform phase compensation based on the initial interference phase and the physical constraint field to obtain a reconstructed continuous phase field.

[0119] An orientation-adaptive phase unwrapping module is used to perform orientation-adaptive unwrapping on the reconstructed continuous phase field to obtain a globally unified continuous phase field.

[0120] The three-dimensional displacement field generation and verification module is used to convert the phase field into a three-dimensional spatiotemporal deformation field and perform accuracy verification and disaster risk assessment.

[0121] Beneficial effects

[0122] The technical solution provided by this invention has the following advantages compared with the prior art:

[0123] This invention generates a physical field that quantifies fault risk and type by integrating UAV lidar topographic data, rock mass integrity coefficient, and manual mining direction; it transforms geological information from post-hoc interpretation to pre-hoc constraints, providing a precise physical basis for subsequent phase processing and preventing misjudgment of deformation direction, such as misjudging subsidence as uplift.

[0124] Furthermore, by using the mining plan to estimate the displacement and guided by the physical constraint field, intelligent compensation is performed on the phase jump in the fracture region to restore the true phase continuity; the periodic ambiguity problem caused by drastic deformation is eliminated directly from the mathematical root, laying a correct mathematical foundation for subsequent untangling.

[0125] In particular, the unwrapping path no longer integrates blindly, but instead prioritizes the direction along the topographic contour lines (perpendicular to the principal gradient direction) to avoid high-risk fault areas; this greatly reduces the risk of error propagation caused by the unwrapping path crossing faults and ensures the reliability of results in continuous regions.

[0126] In addition, by actively identifying and respecting the fault boundary and using it as a barrier to separate the untangled sub-regions, and then performing global calibration, the physical discontinuity of deformation on both sides of the fault plane is maintained (without forced smoothing), and a unified and high-precision deformation result is finally obtained across the entire domain. Attached Figure Description

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

[0128] Figure 1 This is a schematic diagram of the detection method in an embodiment of the present invention. Detailed Implementation

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

[0130] The present invention will be further described below with reference to embodiments.

[0131] Example 1:

[0132] See appendix Figure 1 This case proposes an InSAR surface deformation information detection method based on unmanned aerial vehicles (UAVs), which includes the following steps:

[0133] S1: Raw 3D point cloud data of the mining area is collected using a lidar system mounted on a drone. The discrete points in the raw 3D point cloud data are transformed into a continuous terrain surface. Based on this surface, a terrain gradient vector representing the surface slope and steepness at each location is calculated. Finally, the principal gradient direction angle is extracted as a key guidance parameter. Specifically:

[0134] S101: Using a drone equipped with a lidar system to perform 3D scanning, obtain a dataset of spatial coordinates of discrete points on the surface of the mining area, i.e.:

[0135]

[0136] In the formula, P is the point cloud dataset, which represents the original data carrier of the three-dimensional spatial information of the mining area surface and is the basic object for all subsequent processing. Let be the spatial coordinates of the i-th point cloud, and be a spatial vector. The coordinates are represented by i; i is the index number of the point cloud; N is the total number of point clouds, reflecting the data density and coverage.

[0137] S102: Preprocess the point cloud dataset to generate a clean surface point cloud dataset. Preprocessing includes noise filtering to remove non-surface interference points such as vegetation and buildings; data normalization to unify the coordinate system and elevation benchmark; and density optimization to interpolate and complete sparse areas.

[0138] S103: Processing pure surface point aggregates based on the Delaunay triangulation algorithm to construct continuous terrain surfaces. That is, the elevation function defined in the planar domain.

[0139] S104: Calculate the topographic gradient vector at each spatial location based on the continuous topographic surface and extract the gradient magnitude to quantify the rate of change of landform and slope aspect characteristics, i.e.:

[0140]

[0141] In the formula, is the topographic gradient vector, representing the direction and intensity of the maximum slope of the ground surface at location (x, y); S is the elevation field, a continuous topographic surface generated from UAV lidar data, reflecting the elevation of the ground surface at any point; Let be the partial derivative of the elevation field in the x-direction; Let be the partial derivative of the elevation field in the y-direction.

[0142] S105: The maximum topographic change direction angle used to identify the dominant direction of rock strata fractures based on topographic gradient vector calculation, i.e.:

[0143]

[0144] In the formula, The principal gradient direction angle indicates the direction of the maximum slope of the surface at location (x, y), revealing the direction of the potential sliding tendency of the rock mass; The arctangent function is used to convert gradient ratios into angle values.

[0145] S2: Based on the topographic gradient vector, the rock mass integrity coefficient of the mining area, and the mining direction angle, the risk of abrupt topographic changes is quantified. Then, the interaction effect between artificial mining and natural topography is analyzed through a directional coupling function to obtain a physical constraint field with multi-dimensional control. Specifically:

[0146] S201: Obtaining the distribution map of rock mass integrity coefficient in the mining area based on geological exploration data. The data was then aligned and normalized to obtain the normalized rock series data. ; Obtain the dynamic matrix of mining direction in the mining area based on engineering plan data The rock mass index and mining direction dynamic matrix are unified to the same coordinate system and grid resolution, and a standardized fused data cube is output.

[0147] S202: Construct a response function for rock mass strength attenuation based on fused data cubes, and output the deformation sensitivity to determine the stability of the geological system, i.e.:

[0148]

[0149] In the formula, For deformation sensitivity, the closer the output value is to 1, the higher the risk of deformation; the closer it is to 0, the more stable it is. This is the normalized rock series number.

[0150] S203: Based on the dynamic matrix of the mining direction and the orientation angle of the topographic principal gradient, calculate the directional synergistic factor of the amplification effect of human mining on natural faults, i.e.:

[0151]

[0152] In the formula, For directional coordination factors; The direction of mining advance reflects the mechanical guidance imposed by human activities; the interpolation of the mining advance direction and the direction angle of the main topographic gradient reveals the spatial conflict between natural faults and artificial mining.

[0153] S204: Calculate the physical constraint field mapping potential fault types based on deformation sensitivity and orientation synergy factor, i.e.:

[0154]

[0155] In the formula, The physical constraint field has a value range of [0, 1], integrating the three-dimensional mechanical response of terrain, lithology, and mining; when and When the direction of mining is the same as the main gradient direction, manual mining will exacerbate gravitational instability, causing the rock mass to slide along the slope, resulting in surface subsidence and the formation of a sinkhole; when and When the direction of mining is nearly orthogonal to the direction of the main gradient, the lateral pressure of mining causes the rock strata to be squeezed and displaced, resulting in surface uplift and triggering rock bursts.

[0156] S3: Based on the initial interferometric phase map and physical constraint field, calculate the compensated phase field to eliminate mathematical jumps, thereby compressing the abrupt phase difference in the fault region to... Within this range, it lays the mathematical foundation for continuous differentiability for subsequent untangling. Specifically:

[0157] S301: Accurately align two radar satellite images of the mining area at different time phases, perform pixel-by-pixel conjugate multiplication on the registered composite image data to generate the original interferogram, and extract the complex phase angle components of the original interferogram, i.e., the initial interferometric phase map. .

[0158] S302: Based on the initial interferometric phase diagram The phase gradient magnitude is calculated to quantify the intensity of the abrupt change in the interference phase in space, i.e.:

[0159]

[0160] In the formula, The phase gradient magnitude is calculated by measuring the vector magnitude of the rate of phase change in two directions (x and y axes) to quantify the overall abrupt change in the interferometric phase in space. In a continuous deformation region without faults, the phase gradient magnitude is typically less than 0.1. The phase gradient magnitude changes gradually; the phase entanglement effect causes a sharp increase in the phase gradient magnitude at the fault boundary, and satisfies the condition that the phase gradient magnitude is greater than 1. When the phase difference between adjacent pixels exceeds When, it is determined that there is Integer multiple jump; The phase change rate is in the x-axis direction; This represents the rate of phase change along the y-axis.

[0161] S303: Identifies transition regions based on phase gradient magnitude, providing a spatial target map for subsequent compensation, i.e.:

[0162]

[0163] In the formula, The phase transition domain distribution map is generated by transforming the mathematical discontinuities in a continuous phase field into an operable binary spatial distribution, and the resulting binary spatial distribution map marks the pixel positions that need to be compensated.

[0164] S304: Calculate the estimated surface displacement based on mining plan data, for phase unwrapping. Jump compensation provides physical magnitude constraints to ensure that the compensation value k matches the actual deformation, that is:

[0165]

[0166] In the formula, The estimated value of vertical displacement on the ground surface is used to provide a priori displacement magnitude for phase compensation, avoiding erroneous compensation value k; the phase entanglement phenomenon is essentially a radar observation phase... For periodic mode repetition, the actual phase jump caused by the fault zone is precisely... Since the number of cycles is an integer multiple, k complete cycles are needed to compensate for the untangling phase to ensure continuity. The radar incident angle determines the projection weight of the vertical component; This is the satellite heading angle, used to decompose horizontal displacement onto the radar line-of-sight plane; This represents the east-west horizontal displacement, i.e., the horizontal movement of the Earth's surface along the x-axis. This represents the horizontal displacement in the north-south direction, i.e., the horizontal movement of the Earth's surface along the y-axis. This is the vertical pre-displacement field, reflecting the vertical subsidence of the land surface.

[0167] S3041: Calculate the radius of influence of mining based on the depth of the mining area and the fracture angle of the rock strata, and define the physical range of the mining impact, i.e.:

[0168]

[0169] In the formula, The radius of influence of mining is used to quantify the propagation distance of mining disturbances on the surface. The depth of the mining area determines the thickness of the overlying strata; The rock strata fracture angle is used to control the steepness of subsidence basins.

[0170] S3042: Calculate the vertical pre-displacement field based on the mining boundary and mining thickness of the mining area to predict the vertical subsidence of the ground surface, i.e.:

[0171]

[0172] In the formula, The vertical pre-displacement field represents the predicted vertical subsidence of the Earth's surface at location (x, y). A positive value indicates subsidence, while a negative value indicates uplift. q ​​is the subsidence coefficient, which is a lithological parameter based on experimental calibration. To determine the mining thickness, the maximum subsidence is directly controlled; This is the error function, used to accurately characterize the S-shaped distribution features of subsidence basins; These are the corner points of the mining boundary, used to define the geometric range of the mining area; The radius of influence of the mining operation is the x-axis radius component; The influence radius of the mining operation is the y-axis influence radius component, used to control the expansion range of subsidence along the x / y axis.

[0173] S3043: Based on the vertical pre-displacement field and the radius of influence of mining, calculate the east-west and north-south horizontal displacements, and convert the slope of the subsidence basin into horizontal displacement components, i.e.:

[0174]

[0175] In the formula, b is the horizontal movement coefficient based on experimental calibration, used to adjust the horizontal displacement ratio.

[0176] S305: Convert the estimated vertical displacement of the ground surface into the number of phase compensation cycles for the fault zone. Jump compensation provides a baseline value for the entire cycle, namely:

[0177]

[0178] In the formula, This represents the theoretical number of compensation cycles. It is a rounding function; This refers to the radar wavelength.

[0179] S306: Based on the constraint of the physical constraint field on the theoretical compensation cycle number, the geomechanical mechanism is integrated into the mathematical compensation, that is:

[0180]

[0181] The physical constraint field serves as an indicator of the fault's direction of motion. When its value is greater than or equal to 0, it indicates a normal fault, representing surface expansion. Therefore, the compensation value k is positively compensated to increase the phase value. When its value is less than 0, it indicates a reverse fault, representing surface compression. Therefore, the compensation value k is negatively compensated to decrease the phase value.

[0182] S307: Based on the compensation value, the initial interferometric phase map is improved to a physically solvable continuous phase, which is used to eliminate phase discontinuities in the fault zone, i.e.:

[0183]

[0184] In the formula, This is the reconstructed continuous phase field.

[0185] S4: Based on the reconstructed continuous phase field, an anisotropic diffusion tensor is constructed under physical constraints. An adaptive unwrapping algorithm is used to calculate along the steepest phase change path. Simultaneously, a fault edge protection mechanism is introduced to prevent the path from traversing discontinuous regions, ultimately obtaining a globally continuous unwrapped phase field. Specifically:

[0186] S401: Calculate the phase derivative field based on the reconstructed continuous phase field, generate a phase change gradient map, and identify regions of rapid phase change, i.e.:

[0187]

[0188] In the formula, The phase gradient is in the x-axis direction; The phase gradient is in the y-axis direction; The phase gradient field quantifies the phase change trend and intensity of each pixel. Let be the reconstructed continuous phase field value at position (x, y); Let be the reconstructed continuous phase field value at position (x+1, y). The reconstructed continuous phase field value at position (x, y+1) reflects the coordinate offset of adjacent pixels.

[0189] S402: Based on the phase derivative field and the principal gradient direction angle of the terrain, the directional weights are calculated to obtain the safest set of integration paths orthogonal to the fault strike, avoiding abrupt terrain changes, and constructing the most stable integration path along the direction of material continuity, i.e.:

[0190]

[0191] In the formula, The terrain-adaptive path safety weight coefficient is used to quantify the stability risk of each integral path during phase unwrapping; its value ranges from (0, 1], representing a continuous gradient from absolute danger to absolute safety; when... When the phase change direction is perpendicular to the principal gradient direction angle of the terrain, i.e., along the contour line direction, the corresponding surface material is a continuous region far from the fault; when When the phase change direction is parallel to the direction of the main gradient of the terrain, i.e., along the steep slope direction, it indicates that the current location is in a dangerous area and the phase change direction is parallel to the direction of the main gradient of the terrain, i.e., along the steep slope direction, corresponding to a fault activity zone or a terrain change zone. and This is the unit vector component of the terrain principal gradient direction angle, used to quantize terrain direction into computable mathematical components, i.e., the terrain principal gradient direction angle. Convert to the x and y components of the direction vector.

[0192] S403: Based on the value of the path safety weight coefficient, select the path perpendicular to the direction angle of the terrain principal gradient to construct the optimal path set, and force the path to follow the terrain contour line.

[0193] S404: Calculate the undistorted phase difference between adjacent pixels in the optimal path set, accurately eliminating the phase difference caused by... The jump error caused by periodic folding ensures that no phase ambiguity is introduced during path integration, i.e.:

[0194]

[0195] In the formula, The unbiased phase difference between the i-th and j-th points in the optimal path set is used as the atomic input of the integration path, thus completely eliminating... Risk of sudden propagation; The unbiased phase difference sinusoidal component between point i and point j in the optimal path set. The cosine component of the unbiased phase difference between point i and point j in the optimal path set is used to eliminate the effect of periodic folding and retain the directional information of the phase difference. The cosine reflects the projection length and the sine reflects the rotation direction. It is the arctangent function in the four quadrants; The i-th pixel in the optimal path set; Let j be the pixel at point j in the optimal path set.

[0196] S405: Based on the undistorted phase difference, the phase difference value is accumulated and integrated along the safe path with a preset reference point as the reference, generating a continuous unwrapped phase field inside the path in the fault zone, eliminating... The risk of transition error propagating along the integration path, namely:

[0197]

[0198] In the formula, The absolute phase value after unwrapping the target point represents the location of the surface point. The actual deformability at the location; The reference point phase value of the reconstructed continuous phase field, i.e., the absolute phase reference point of the integration process; It is a directional accumulation operator along the optimal path set.

[0199] S406: Calculate the fault barrier function using the physical constraint field and the jump boundary distance to obtain a real-time path interruption command, blocking the integral path from crossing the fault line, i.e.:

[0200]

[0201] In the formula, It is a geological fault activity intensity field, used to quantitatively describe the fault rupture risk at any point on the Earth's surface; The Euclidean distance from the current point (x, y) to the fault boundary; The boundary of the jump domain is the topological defect boundary in the phase field, indicating the occurrence of entanglement phase. The set of continuous curves with integer multiple jumps is obtained as follows:

[0202] S4061: Curl field for calculating phase gradient based on phase derivative field, accurately capturing the curl field in the phase field. Integer multiple jump defects provide an unambiguous criterion for jump boundary extraction, namely:

[0203]

[0204] In the formula, The curl field of the phase gradient reflects the loop integral value of the closed path; The starting gradient is in the x-axis direction; The gradient is the upper right gradient along the y-axis. For gradient compensation in the x-axis direction; Gradient compensation for the starting point in the y-axis direction.

[0205] S4062: Set the boundary of the jump domain based on the curl field of the phase gradient, that is:

[0206] .

[0207] S407: When the value of the geological fault activity intensity field is less than the preset judgment threshold, the current path is forcibly terminated to prevent the integral path from crossing the fault line and maintain the physical authenticity of the discontinuous area; based on the geological fault activity intensity field, the entire domain is divided into independent continuous sub-regions, that is, multi-sub-region unwrapping phase fields, thereby ensuring that the unwrapping results on both sides of the fault maintain physical independence and that the discontinuity of the fault zone deformation is not smoothed out.

[0208] S408: Global continuity correction based on multi-sub-region unwrapped phase field, calculating period offset.

[0209]

[0210] In the formula, To determine the optimal cycle offset; To unwrap the phase field of the Ath sub-region in the multi-sub-region unwrapping phase field; To unwrap the phase field of the Bth sub-region in the multi-sub-region unwrapping phase field; The absolute phase offset between two sub-regions A and B This is the operator for obtaining the median value of a sample.

[0211] S409: Calculates a continuous phase field with global reference alignment based on the unwrapped phase field of multiple sub-regions and the statistically optimal periodic offset, providing a distortion-free data foundation for deformation inversion, namely:

[0212]

[0213] In the formula, To achieve the final transformation from relative phase observation to absolute deformation field by establishing a unified global benchmark for the continuous phase field; The phase field is used to untangle the multi-sub-region.

[0214] S5: This method combines a unified, continuous phase field transformation across the entire domain with a point cloud dataset. Through 3D projection transformation, phase information is converted into spatial displacement vectors. Multi-period monitoring data is then fused to generate a time-series deformation field, and accuracy is cross-validated using ground monitoring stations. Ultimately, a 3D spatiotemporal evolution deformation field with millimeter-level accuracy is obtained, providing quantifiable and verifiable decision-making support for fault activity assessment and disaster early warning. Specifically:

[0215] S501: Calculates the radar line-of-sight displacement field based on a unified global benchmark continuous phase field, expressing the change in the Earth's surface along the radar observation direction, i.e.:

[0216]

[0217] In the formula, The radar line-of-sight displacement vector represents the linear displacement of a point on the ground along the direction of the satellite radar beam. A gradient operator for a continuous phase field with a unified global reference, quantizing the phase change rate between adjacent pixels.

[0218] S502: Based on point cloud datasets, the radar line-of-sight displacement vector is decomposed into vertical, east-west, and north-south displacement components to construct a three-dimensional deformation base, providing spatial omnidirectional diagnostic evidence for fault activity patterns and disaster risks.

[0219]

[0220] In the formula, This represents the vertical displacement component; This represents the east-west displacement component; This represents the displacement component in the north-south direction; This is the inverse of the satellite geometric projection matrix, used to map one-dimensional LOS displacement to a three-dimensional coordinate system; It is a three-dimensional displacement field.

[0221] S503: Collect the three-dimensional displacement fields of the mining area at different observation times and merge them into a three-dimensional displacement field set. It describes the evolution of surface deformation over time; where n is the number of three-dimensional displacement fields.

[0222] S504: Calculates the displacement rate field based on a three-dimensional displacement field ensemble, used for real-time monitoring of the accelerating trend of surface subsidence / uplift, i.e.:

[0223]

[0224] In the formula, This is a displacement rate field used to quantify the trend of changes in surface motion; and Represents the three-dimensional displacement field set in adjacent time phases and The displacement observation values.

[0225] S505: Calculate the cumulative displacement field based on the displacement rate field to reveal the spatiotemporal evolution law of deformation, namely:

[0226]

[0227] In the formula, To accumulate the displacement field, the displacement from the initial time was quantized. The cumulative displacement up to the current time t.

[0228] S506: By combining the measured displacement from GNSS monitoring stations, residual correction is applied to the cumulative displacement field to generate a calibration displacement field with millimeter-level accuracy, i.e.:

[0229]

[0230] In the formula, The magnitude of the residual vector is used to quantify the error between the cumulative displacement field and the measured ground data. This is the measured displacement field; This is the corrected final displacement field.

[0231] S507: Based on the physical constraint field and the corrected final displacement field, calculate the differential movement of the fault, quantitatively diagnose the fault activity mode and the relative displacement intensity of the hanging wall and footwall, and provide key mechanical basis for the risk of fault activation in mines, namely:

[0232]

[0233] In the formula, The vertical displacement on both sides of the fault is used to quantify the relative displacement between the hanging wall and footwall of the fault. This is a sign function used to determine whether a fault is normal or reverse.

[0234] S508: Calculate the activation risk index based on the vertical displacement on both sides of the fault to achieve quantitative disaster early warning, i.e.:

[0235]

[0236] In the formula, It is a fault activation risk index used to comprehensively assess the probability of disasters.

[0237] Example 2:

[0238] Building upon Example 1, this paper proposes an InSAR surface deformation information detection system based on unmanned aerial vehicles (UAVs), comprising a terrain gradient extraction module, a geomechanical constraint field construction module, a phase jump compensation module, a direction-adaptive phase unwrapping module, and a three-dimensional displacement field generation and verification module. Among these:

[0239] The terrain gradient extraction module collects point cloud data using UAV LiDAR, preprocesses it to construct a continuous terrain surface, and calculates the terrain gradient vector and principal gradient direction angle. This quantifies the surface slope and steepness, reveals the potential slip direction of the rock mass, and provides a topographic basis for fault analysis.

[0240] The geomechanical constraint field construction module integrates geological exploration data and mining direction data, calculates deformation sensitivity and directional synergy factors, and finally generates a physical constraint field. By combining the influence of geological stability and human mining activities, it predicts the risk of fault types such as collapse or rockburst.

[0241] The phase jump compensation module generates an initial interferometric phase map from radar imagery and calculates compensation values ​​based on physical constraint fields to eliminate phase jumps in the fault region. This solves the phase ambiguity problem caused by discontinuous deformation and lays a continuous and differentiable mathematical foundation for subsequent unwrapping.

[0242] The orientation-adaptive phase unwrapping module selects a safe integration path based on the terrain principal gradient direction, then unwrappes the phase along the contour line direction, and finally blocks dangerous paths through a fault barrier function. This achieves a globally continuous and accurately unwrapped phase field while maintaining the physical discontinuity of the fault.

[0243] The 3D displacement field generation and verification module decomposes phase information into 3D displacement components and fuses multi-period data to generate a spatiotemporal evolution deformation field. Finally, it verifies the accuracy using ground station data. By outputting 3D deformation fields and fault differential motion quantities with millimeter-level accuracy, it provides a quantitative basis for disaster early warning decisions.

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

Claims

1. A method for detecting InSAR surface deformation information based on unmanned aerial vehicles (UAVs), characterized in that, Includes the following steps: S1: Collect point cloud data of the mining area using the lidar system carried by the drone, construct a continuous terrain surface, calculate the terrain gradient vector based on the surface, and extract the main gradient direction angle of the terrain. S2: Based on the topographic gradient vector, the rock mass integrity coefficient of the mining area, and the mining direction angle, a physical constraint field is constructed to quantify the risk of terrain abrupt change and the interaction effect between mining and terrain; Step S2 includes: S201: Obtain the distribution map of rock mass integrity coefficient and dynamic matrix of mining direction in the mining area, and perform data fusion; S202: Based on the fused data, construct a deformation sensitivity function, namely: ； In the formula, Deformation sensitivity; This represents the terrain gradient vector; The normalized rock series number; S203: Based on the dynamic matrix of the mining direction and the principal gradient direction angle of the terrain, calculate the directional synergy factor, i.e.: ； In the formula, For directional coordination factors; This indicates the direction of mining operations; The direction angle of the principal gradient of the terrain; S204: Based on the deformation sensitivity and the orientation coordination factor, calculate the physical constraint field to indicate the fault type and activity level, i.e.: ； In the formula, For physical constraint fields; For directional coordination factors; Deformation sensitivity; S3: Based on the initial interferometric phase diagram and the physical constraint field, calculate the compensation value and compensate the initial interferometric phase diagram to obtain the reconstructed continuous phase field; S4: Based on the reconstructed continuous phase field and the physical constraint field, the phase is unwrapped using a direction-adaptive unwrapping algorithm to obtain a continuous phase field with a unified global reference. S5: Transform the global benchmark unified continuous phase field into a three-dimensional spatiotemporal evolution deformation field, and perform accuracy verification and disaster risk assessment.

2. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 1, characterized in that, Step S1 further includes: S101: Obtain the original 3D point cloud dataset of the mining area; S102: Preprocess the original 3D point cloud dataset to obtain a clean surface point cloud dataset; S103: Construct a continuous terrain surface based on the aforementioned pure surface point cluster. ; S104: Based on the continuous terrain surface, calculate the terrain gradient vector at each location, i.e.: ； In the formula, S is the topographic gradient vector; S is the elevation field, i.e., the continuous topographic surface. ; Let be the partial derivative of the elevation field in the x-direction; Let be the partial derivative of the elevation field in the y-direction; S105: Based on the terrain gradient vector, calculate the principal gradient direction angle of the terrain, i.e.: ； In the formula, The direction angle of the principal gradient of the terrain; It is the arctangent function.

3. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 2, characterized in that, Step S3 includes: S301: Generate initial interferometric phase map ; S302: Calculate the phase gradient magnitude of the initial interferometric phase pattern, i.e.: ； In the formula, The magnitude of the phase gradient; S303: Identify the phase transition region based on the phase gradient magnitude, that is: ； In the formula, This is a distribution diagram of the phase jump domain; S304: Based on mining plan data, calculate the estimated surface displacement, i.e.: ； In the formula, Estimated value of vertical displacement of the Earth's surface; The radar incident angle; This refers to the satellite's heading angle; This represents a horizontal displacement in the east-west direction. This represents a north-south horizontal displacement. For vertical pre-displacement field; S305: Convert the estimated surface displacement into the theoretical number of compensation cycles, i.e.: ； In the formula, This represents the theoretical number of compensation cycles. It is a rounding function; The radar wavelength; S306: Based on the physical constraint field, the theoretical compensation cycle number is constrained to obtain the final compensation value, i.e.: ； In the formula, k is the compensation value; S307: The initial interferometric phase map is compensated using the final compensation value to obtain the reconstructed continuous phase field, i.e.: ； In the formula, This is the reconstructed continuous phase field.

4. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 3, characterized in that, In step S304, the estimated surface displacement is calculated as follows: S3041: Calculate the radius of influence of mining based on the depth of the mining area and the fracture angle of the rock strata, i.e.: ； In the formula, The radius of influence of mining activities; The depth of the mining area; Rock strata fracture angle; S3042: Calculate the vertical pre-displacement field based on the mining boundary and mining thickness of the mining area to predict the vertical subsidence of the ground surface, i.e.: ； In the formula, The vertical pre-displacement field is represented by q, where q is the subsidence coefficient. For the thickness of the ore; It is the error function; This refers to the corner point of the mining boundary; The radius of influence of the mining operation is the x-axis radius component; The radius of influence of the mining operation is the radius of influence component along the y-axis; S3043: Calculate the east-west and north-south horizontal displacements based on the vertical pre-displacement field and the radius of influence of mining, i.e.: ； In the formula, b is the horizontal movement coefficient based on experimental calibration.

5. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 3, characterized in that, Step S4 includes: S401: Calculate the phase derivative field of the reconstructed continuous phase field, i.e.: ； In the formula, The phase gradient is in the x-axis direction; The phase gradient is in the y-axis direction; For phase gradient field; Let be the reconstructed continuous phase field value at position (x, y); S402: Based on the phase derivative field and the terrain principal gradient direction angle, calculate the path safety weight coefficient to determine the optimal integration path and merge them into the optimal path set, i.e.: ； In the formula, For terrain-adaptive path safety weighting coefficients; S403: Based on the optimal path set, calculate the unbiased phase difference between adjacent pixels, i.e.: ； In the formula, The unambiguous phase difference between the i-th and j-th points in the optimal path set; The i-th pixel in the optimal path set; The j-th pixel in the optimal path set; S404: Perform phase integration along the optimal integration path to obtain the preliminary unwrapped phase field, i.e.: ； In the formula, The absolute phase value after untangling at the target point; The reference point phase value for the reconstructed continuous phase field; For directional accumulation operators along the optimal path set; S405: Based on the aforementioned physical constraint field and phase jump boundary, construct the fault barrier function, namely: ； In the formula, This represents the intensity field of geological fault activity. The Euclidean distance from the current point (x, y) to the fault boundary; For the boundary of the jump domain; For physical constraint fields; S406: Based on the fault barrier function, the entire domain is divided into multiple continuous sub-regions for independent unwrapping to obtain the multi-sub-region unwrapped phase field; S407: Perform global continuity correction on the unwrapped phase field of the multi-sub-region to obtain a global reference-unified continuous phase field, i.e.: ； In the formula, A continuous phase field with a unified global reference. Unwrap the phase field for the multi-sub-region; To determine the optimal cycle offset.

6. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 5, characterized in that, The boundary of the transition domain is calculated as follows: S4051: The curl field for calculating the phase gradient based on the phase derivative field, i.e.: ； In the formula, The curl field represents the phase gradient; The starting gradient is in the x-axis direction; The gradient is the upper right gradient along the y-axis. For gradient compensation in the x-axis direction; Gradient compensation for the starting point in the y-axis direction; S4052: Set the boundary of the jump domain based on the curl field of the phase gradient, that is: ； In the formula, This is the boundary of the jump domain.

7. The method for detecting InSAR surface deformation information based on unmanned aerial vehicles according to claim 1, characterized in that, Step S5 includes: S501: The continuous phase field of the global reference is transformed into a radar line-of-sight displacement field, that is: ； In the formula, This is the displacement vector along the radar line of sight; The gradient operator for a continuous phase field with a unified global reference. The radar wavelength; S502: Decompose the radar line-of-sight displacement field into three-dimensional displacement components in the vertical, east-west, and north-south directions to construct a three-dimensional displacement field, namely: ； In the formula, This represents the vertical displacement component; This represents the east-west displacement component; This represents the displacement component in the north-south direction; It is the inverse of the satellite's geometric projection matrix; It is a three-dimensional displacement field; S503: Construct a time-series three-dimensional displacement field set based on multi-period monitoring data; S504: Calculate the displacement rate field and cumulative displacement field based on the aforementioned three-dimensional displacement field set, i.e.: ； In the formula, For displacement velocity field; and Represents the three-dimensional displacement field set in adjacent time phases and Displacement observations; ； In the formula, To accumulate the displacement field, the displacement from the initial time was quantized. The cumulative displacement up to the current time t; S505: The cumulative displacement field is corrected using data from ground monitoring stations to obtain the final displacement field, i.e.: ； In the formula, The magnitude of the residual vector; This is the measured displacement field; This is the corrected final displacement field; S506: Based on the final displacement field and the physical constraint field, calculate the fault differential motion and activation risk index, i.e.: ； In the formula, This represents the vertical displacement on both sides of the fault; This is a sign function used to determine normal / reverse faults; ； In the formula, This is the fault activation risk index.

8. A system for implementing the UAV-based InSAR surface deformation information detection method according to any one of claims 1-7, characterized in that, include: The terrain gradient extraction module is used to collect and process point cloud data of the mining area through the lidar system carried by the UAV, and to calculate the terrain gradient vector and the terrain principal gradient direction angle. The geomechanical constraint field construction module is used to construct physical constraint fields based on topographic data, rock mass data, and mining data. The phase jump compensation module is used to perform phase compensation based on the initial interference phase and the physical constraint field to obtain a reconstructed continuous phase field. An orientation-adaptive phase unwrapping module is used to perform orientation-adaptive unwrapping on the reconstructed continuous phase field to obtain a globally unified continuous phase field. The three-dimensional displacement field generation and verification module is used to convert the phase field into a three-dimensional spatiotemporal deformation field and perform accuracy verification and disaster risk assessment.

Images