A method for calculating three-dimensional deformation field of coal mine goaf based on Beidou constraint InSAR

Abstract

The application discloses a kind of InSAR coal mine goaf three-dimensional deformation field resolving methods based on Beidou constraint, comprising: obtaining the LOS direction deformation rate of coal mine goaf and pre-processing ascending track and descending track;Based on the geometric relationship equation of the LOS direction deformation rate of ascending track and descending track and the deformation rate of three directions of east-west, north-south and vertical is constructed to obtain the three-dimensional deformation rate of InSAR;According to the observation data of each Beidou monitoring station in the coal mine goaf, the three-dimensional deformation rate of each Beidou station is obtained;Based on the data fidelity term of InSAR observation value, the constraint term based on Beidou observation value and the edge maintaining regularization term for stable solution, the objective function of three-dimensional deformation field resolving is constructed;Based on the three-dimensional deformation rate of InSAR and the three-dimensional deformation rate of each Beidou station, the objective function is solved using alternating direction multiplier method, and the three-dimensional deformation rate of each geographic grid point is obtained.The application realizes the high-precision three-dimensional deformation field modeling of coal mine goaf.

Description

Technical Field

[0001] This invention relates to the field of surveying and mapping science and technology, specifically to an InSAR method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints. Background Technology

[0002] Underground goaf areas formed by coal mining are highly susceptible to geological disasters such as surface subsidence and collapse, posing a serious threat to mining infrastructure, the ecological environment, and the safety of people's lives and property. Therefore, continuous and accurate monitoring of surface deformation in coal mine goaf areas is crucial.

[0003] Currently, InSAR technology has become the mainstream technology for deformation monitoring in mining areas due to its all-weather, all-day, wide-area, and high spatial resolution monitoring capabilities. However, traditional InSAR technologies (such as D-InSAR, PS-InSAR, and SBAS-InSAR) mainly monitor one-dimensional deformation along the radar line-of-sight (LOS). To obtain the true three-dimensional deformation vector (east-west, north-south, and vertical directions), at least three SAR datasets with different geometric perspectives (such as ascending orbit, descending orbit, and near-polar orbit) are usually required, and decomposed through mathematical models.

[0004] This method has obvious limitations:

[0005] ① High cost: It requires the acquisition of multi-source SAR data, which greatly increases data cost and processing complexity.

[0006] ② Geometric weakening: North-south deformation is least sensitive to SAR geometry. The accuracy of north-south deformation calculated solely from multi-angle InSAR data is extremely low, or even unreliable.

[0007] ③ Unstable solution: When the available SAR data has poor geometric configuration or the deformation signal is weak, the ill-conditioned problem of the three-dimensional decomposition equation becomes prominent, and the solution result has a large error.

[0008] The BeiDou Navigation Satellite System (BDS) can provide three-dimensional surface deformation information with millimeter-level to centimeter-level accuracy, but its drawback is low spatial resolution (only point observations). Summary of the Invention

[0009] To overcome the shortcomings of existing technologies that rely solely on InSAR technology, such as high cost, geometric weakening, and unstable solution, as well as the low spatial resolution of those that rely solely on BeiDou navigation technology, this invention provides a BeiDou-constrained InSAR method for solving the three-dimensional deformation field of coal mine goaf. By combining the advantages of InSAR's high-density planar observation with the advantages of BeiDou's high-precision point-based three-dimensional observation, the invention overcomes the shortcomings of each method and achieves high-precision three-dimensional deformation field modeling of coal mine goaf.

[0010] According to one aspect of the present invention, a method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints using InSAR is provided, comprising: acquiring and preprocessing the LOS deformation rates of the ascending and descending orbits of the coal mine goaf; constructing geometric relationship equations between the LOS deformation rates of the ascending and descending orbits and the deformation rates in the east-west, north-south, and vertical directions based on the preprocessing results to obtain the three-dimensional deformation rate of InSAR; calculating the three-dimensional deformation rate of each BeiDou station based on the observation data of each BeiDou monitoring station in the coal mine goaf; constructing an objective function for calculating the three-dimensional deformation field based on the data fidelity term of the InSAR observations, the constraint term based on the BeiDou observations, and the edge-preserving regularization term for stable calculation; and solving the objective function using the alternating direction multiplier method based on the three-dimensional deformation rate of InSAR and the three-dimensional deformation rate of each BeiDou station to obtain the three-dimensional deformation rate of each geographic grid point.

[0011] Furthermore, the expression for the geometric relation equation is:

[0012] ,

[0013] in, For LOS deformation rate, Let be the unit geometric vector in the direction of the satellite's line of sight. , , Let be the deformation rates in the east-west, north-south, and vertical directions, respectively.

[0014] Furthermore, the expression for the objective function is:

[0015] ,

[0016] in, Let be the objective function. For InSAR observation vectors, The geometric coefficient matrix, Let be the three-dimensional deformation vector to be determined. , , , Let be the deformation rates in the east-west, north-south, and vertical directions, respectively. For BeiDou observation vectors, For the sampling matrix of Beidou stations, This represents the overall weighting coefficient for the BeiDou constraint term. To preserve the regularization term for the edges, The weight coefficients for the edge-preserving regularization term.

[0017] Furthermore, the expression corresponding to the edge-preserving regularization term is:

[0018] ,

[0019] in, To preserve the regularization term for the edges, For deformation component index, This represents the total number of geographic grid points. Indices representing geographic grid points For gradient operators, It is an L1 norm. Three-dimensional deformation vector The c-th component, and These are first-order difference operators along the horizontal and vertical directions, respectively.

[0020] Furthermore, the objective function is solved using the alternating direction multiplier method to obtain the three-dimensional deformation rate of each geographic grid point, including: splitting the objective function into a function containing a data fidelity term and a BeiDou constraint term. Norm subproblem, and a problem involving margin-preserving regularization. The norm subproblem can be solved using the conjugate gradient method. Norm subproblem and its solution using the soft-threshold shrinkage operator The norm subproblem involves iteratively updating the 3D deformation vector, auxiliary variables, and dual variables in the objective function until convergence, to obtain the 3D deformation rate for each geographic grid point.

[0021] Furthermore, the acquisition of the LOS-direction deformation rate of the ascending and descending orbits includes: using the observation data from each BeiDou monitoring station to perform tropospheric atmospheric delay correction on the acquired SAR image data of the ascending and descending orbits; and using time-series InSAR technology to process the corrected SAR image data of the ascending and descending orbits to obtain the LOS-direction deformation rate of the ascending and descending orbits.

[0022] Furthermore, the solution method also includes: step S5, outputting and verifying the three-dimensional deformation rate of each geographic grid point.

[0023] According to one aspect of this invention, a BeiDou-constrained InSAR three-dimensional deformation field calculation system for coal mine goaf areas is provided. This system is implemented using the aforementioned BeiDou-constrained InSAR three-dimensional deformation field calculation method for coal mine goaf areas. The system includes: a LOS-direction deformation rate acquisition module, used to acquire and preprocess the LOS-direction deformation rates of the ascending and descending orbits of the coal mine goaf area; and an InSAR three-dimensional deformation rate acquisition module, used to construct geometric relationship equations between the LOS-direction deformation rates of the ascending and descending orbits and the deformation rates in the east-west, north-south, and vertical directions based on the preprocessing results, to obtain the InSAR three-dimensional deformation rate. The system includes: a BeiDou 3D deformation rate acquisition module, used to calculate the 3D deformation rate of each BeiDou monitoring station based on the observation data of each BeiDou monitoring station in the coal mine goaf; an objective function construction module, used to construct the objective function for 3D deformation field calculation based on the data fidelity term of InSAR observations, the constraint term based on BeiDou observations, and the edge-preserving regularization term for stable calculation; and a comprehensive 3D deformation rate calculation module, used to solve the objective function using the alternating direction multiplier method based on the 3D deformation rate of InSAR and the 3D deformation rate of each BeiDou station to obtain the 3D deformation rate of each geographic grid point.

[0024] According to one aspect of the present invention, an electronic device is provided, including a memory and a processor, wherein the memory stores program instructions that are executed by the processor, and the processor invokes the program instructions to execute the method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints using InSAR.

[0025] According to one aspect of the present invention, a non-transitory computer-readable storage medium is provided, the non-transitory computer-readable storage medium storing computer instructions that cause the computer to execute the aforementioned method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints using InSAR.

[0026] The above technical solution first acquires the ascending and descending orbit InSAR time-series LOS deformation rate and sparse BeiDou three-dimensional deformation rate data; then, geocoding and resampling the ascending and descending orbit InSAR time-series LOS deformation rate to establish a geometric model of InSAR LOS observation and three-dimensional deformation; finally, constructing an objective function that includes an InSAR data fidelity term, a BeiDou constraint term, and a spatial regularization term, and limiting the spatial regularization term to a total variational (norm) regularization term that preserves the edges; and accordingly, using the alternating direction multiplier method (ADMM) and a soft threshold shrinkage operator, iteratively solving the objective function to obtain the three-dimensional deformation field.

[0027] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0028] (1) The present invention limits the spatial regularization term to a total variational (norm) regularization term that preserves the edge; and accordingly, it uses the alternating direction multiplier method (ADMM) and the soft threshold shrinkage operator to iteratively solve the objective function to obtain the three-dimensional deformation field, so as to solve the technical contradiction of blurred edges of goaf subsidence caused by traditional smoothing regularization (such as Laplace smoothing).

[0029] (2) The present invention removes the main noise source (atmospheric delay) in the InSAR observations by using the observation data of the Beidou monitoring station, which significantly improves the input accuracy of the data fidelity term, thereby making the final calculated three-dimensional deformation field more reliable and accurate.

[0030] (3) This invention combines the advantages of InSAR planar high-density observation with the advantages of Beidou point-based high-precision three-dimensional observation, overcoming the shortcomings of existing technologies that only use InSAR technology, such as high cost, geometric weakening and unstable solution, or the shortcomings of Beidou navigation technology, such as low spatial resolution, and realizes high-precision three-dimensional deformation field modeling of coal mine goaf. Attached Figure Description

[0031] 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 accompanying drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0032] Figure 1 A flowchart illustrating a method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints, provided as an embodiment of the present invention.

[0033] Figure 2 This is a schematic diagram illustrating the principle of a method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints, provided in an embodiment of the present invention. Detailed Implementation

[0034] It should be noted that:

[0035] The terms “comprising” and “having”, and any variations thereof, in the specification, claims, and accompanying drawings of this invention are intended to cover a non-exclusive inclusion, such as a process, method, system, product, or apparatus that includes a series of steps or units, not necessarily limited to those explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0036] The block diagrams shown in the accompanying drawings are merely functional entities and do not necessarily correspond to physically independent entities. That is, these functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor devices and / or microcontroller devices. The flowcharts shown in the accompanying drawings are merely illustrative and do not necessarily include all content and operations / steps, nor do they necessarily have to be performed in the described order. For example, some operations / steps can be decomposed, while others can be combined or partially combined; therefore, the actual execution order may change depending on the specific circumstances.

[0037] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined to form new technical solutions. Such combinations are not bound by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0038] Example 1:

[0039] Please refer to the appendix for details. Figure 1 and 2 This embodiment provides a method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints using InSAR. This embodiment uses a goaf in a large coal mine as the monitoring area. The following will demonstrate the existing techniques employed... The specific steps of regularization methods:

[0040] Step S1: Acquire SAR image data of the coal mine goaf area for both ascent and descent, and process the data using time-series InSAR technology to obtain the LOS deformation rate of the ascent and descent orbits; simultaneously acquire observation data from each BeiDou monitoring station within the coal mine goaf area, and calculate the three-dimensional deformation rate of each BeiDou station.

[0041] Step S1 aims to acquire and preprocess data. Specifically, 50 Sentinel-1A ascending-orbit SAR images and 50 descending-orbit SAR images of the coal mine goaf from 2019 to 2021 are collected. Temporal InSAR technology is used to process the data to obtain the spatiotemporally continuous annual average LOS deformation rate for both ascending and descending orbits. and Data from 15 BeiDou monitoring stations within and around the coal mine's goaf were collected, and the three-dimensional deformation rate of each monitoring station (including the east-west deformation rate calculated by BeiDou) was calculated. North-South Horizontal Deformation Rate Vertical deformation rate ).

[0042] To obtain high-density surface LOS-oriented deformation fields, time-series InSAR technology is required. The specific technology chosen (SBAS-InSAR or PS-InSAR) depends on the surface cover characteristics of the monitoring area.

[0043] PS-InSAR (Permanent Scatterer InSAR): This technique prioritizes point targets (such as buildings, exposed rocks, artificial corner reflectors, and infrastructure) that maintain high coherence over long time series. If there are many such stable structures on the surface of a mining area, PS-InSAR can provide millimeter-level, high-precision time series of point deformation.

[0044] SBAS-InSAR (Short Baseline Set InSAR): This technique is more suitable for areas with surface cover dominated by distributed scatterers (such as bare soil, sparsely vegetated areas, and farmland), which often lack PS points. By combining interferogram pairs with short spatiotemporal baselines, SBAS can maintain high coherence, thereby obtaining deformation fields with higher density and wider spatial coverage in low coherence areas.

[0045] Considering that the surface of coal mine goaf areas is mostly bare soil, farmland, or vegetation cover, stable PS point sources are relatively sparse. To obtain a more spatially continuous and comprehensive deformation field, this embodiment preferably employs SBAS-InSAR technology. By performing SBAS time-series processing on the collected ascending and descending orbit Sentinel-1A images (50 scenes each), two high-spatial-density InSAR ascending orbit LOS-direction deformation rates covering the entire monitoring area were finally obtained. and InSAR de-orbiting LOS deformation rate .

[0046] Step S2: Convert the LOS values ​​of the ascending and descending orbits to the deformation rate in a unified manner to the geodetic coordinate system and perform geographic grid resampling.

[0047] In this embodiment, the LOS deformation rates of the ascending and descending orbits are unified to the UTM-WGS84 coordinate system and resampled to a 50m×50m grid.

[0048] Step S3: For each geographic grid point, establish the geometric relationship equation between the LOS deformation rate of the ascending and descending orbits and the deformation rates in the east-west, north-south, and vertical directions to obtain the three-dimensional deformation rate of InSAR.

[0049] In step S3, for each geographic grid point Establish the SAR observation geometry equation:

[0050]

[0051] in, For geographic grid points InSAR ascending orbit LOS deformation rate, For geographic grid points The deformation rate of the LOS direction of the InSAR descent orbit. , , The geographic grid points to be found Deformation rates in the east-west, north-south, and vertical directions; This is the unit geometric vector in the direction of the satellite's line of sight.

[0052] In this step, the goal is to provide each geographic grid point with... Establish its InSAR observations (LOS to deformation rate) ) and the actual three-dimensional deformation rate to be determined ( Mathematical transformation model between ( ).

[0053] The geometric equation is a linear projection model:

[0054] ,

[0055] Or write:

[0056] ,

[0057] in, It is the deformation rate in the radar line-of-sight direction (i.e., the LOS deformation rate from step S1). It is the geographic grid point Three unknowns to be determined. It is a unit geometric vector describing the geometry of satellite observation in the line-of-sight direction.

[0058] This unit geometric vector is for each geographic grid point. And each orbit (ascending / descending) is different; it is calculated based on the metadata of the SAR image (such as orbital parameters and imaging geometry):

[0059] ,

[0060] ,

[0061] ,

[0062] in, , , These are geographic grid points Unit vectors in the east-west, north-south, and vertical directions. Geographic grid points The radar beam incident angle at that location, It is the satellite flight azimuth angle (or orbital inclination-related azimuth angle) at grid point i.

[0063] In step S3, it is necessary to calculate the respective geometric vectors for the ascending orbit (asc) and descending orbit (desc) datasets. and Thus, for geographic grid points that simultaneously have both ascending and descending orbit observations... This yields the SAR observation geometry equations, which are used to solve for the three unknowns. It provides the foundation.

[0064] Step S4: Construct a three-dimensional deformation field solution model that integrates BeiDou constraints. The objective function of the three-dimensional deformation field solution model includes a data fidelity term based on InSAR observations, a constraint term based on BeiDou observations, and an edge-preserving regularization term for stable solution. Based on the three-dimensional deformation rate of InSAR and the three-dimensional deformation rate of each BeiDou station, the objective function is solved using the alternating direction multiplier method to obtain the three-dimensional deformation rate of each geographic grid point.

[0065] Step S41, construct the penalized least squares objective function, the corresponding expression is:

[0066]

[0067] in, The objective function (i.e., the cost function) to be optimized. For InSAR observation vectors (i.e., all) (a set) The geometric coefficient matrix (by) constitute), Let be the three-dimensional deformation vector to be determined. , For BeiDou observation vectors (i.e., all) The set, (This refers to the three-dimensional deformation rate observation value of the k-th Beidou monitoring station). For the sampling matrix of Beidou stations, This represents the overall weighting coefficient for the BeiDou constraint term. For regularization functions, represents the weighting coefficients of the space regularization term.

[0068] In this embodiment, the three-dimensional velocity of 15 BeiDou stations is used as the constraint point ( Set constraint weights The closer a grid point is to a BeiDou station, the stronger its binding force. (The larger the value). The Laplacian operator is chosen as the smoothing constraint. Smoothing factor Determined using the L-curve method. The traditional method employed... Smoothing (such as the Laplace operator) Its goal is to minimize the square of the gradient (the rate of change of deformation). When encountering a real, large gradient value (such as at the edge of a goaf), The norm will impose a severe penalty, forcing the model to smooth it out, resulting in the edges of the goaf in the solution being severely blurred, making it impossible to accurately delineate the disaster boundary.

[0069] When constructing the objective function in step S4, the BeiDou constraint term Weighting coefficients in It is a key parameter that determines the extent to which BeiDou high-precision data can boost InSAR data.

[0070] In this embodiment, It is not a single, global scalar, but rather an adjustable weight function related to data quality and spatial location. Its value should simultaneously consider the accuracy of the BeiDou observations and the spatial influence (i.e., distance) of this constraint on neighboring grid points.

[0071] 1. Weights based on observation accuracy

[0072] Each BeiDou monitoring station Both provide a three-dimensional rate and its accuracy (standard deviation) In the model, It should be implemented as a BeiDou constrained weight matrix. Its diagonal elements (i.e., the first) The weight of each BeiDou station is inversely proportional to the square of its accuracy:

[0073] ,

[0074] in, Let k be the observation standard deviation of the BeiDou station. This refers to the weighting based on the accuracy of BeiDou station k. This means that the higher the accuracy calculated from BeiDou data (…), the higher the weighting. The smaller the value, the higher the weight the site receives in the objective function. The larger it is, the stronger its binding force will be.

[0075] 2. Spatial distance-based weighting (soft constraint)

[0076] The BeiDou constraint should not only apply to the grid point where the BeiDou site is located, but also influence its neighboring grid points in the form of a soft constraint, and this influence should decrease with distance. In the implementation (of Example 1), this is explicitly stated as follows: the closer a grid point is to the BeiDou site, the stronger its constraint. This can be achieved through a spatial interpolation function or a distance-weighted function. For example, for any grid point that is not a BeiDou site... The constraint weights it is subject to It could be based on the nearest BeiDou site. distance Functions, such as the Gaussian decay function:

[0077] ,

[0078] in, Grid points of non-BeiDou stations Constraint weights based on spatial distance, It is the maximum weight value of the BeiDou constraint term. Grid points of non-BeiDou stations With the latest BeiDou site distance, This is the characteristic distance of spatial weight decay, used to control the range of influence. For indexes of grid points or BeiDou sites.

[0079] In a better implementation, these two factors are combined to construct a spatially continuously varying BeiDou weight field. At the BeiDou site Location, weight Due to its high precision ( The weight is determined by the location, while in areas far from the site, the weight is determined by the location. Then based on distance It rapidly decays to a very small value (or zero). This ensures that BeiDou data, while providing high-precision anchor points, does not unduly constrain areas far from the stations that should be dominated by InSAR data.

[0080] Step S42: Construct a large system of equations for all geographic grid points (approximately 100,000) in the entire coal mine goaf area, and solve iteratively using the conjugate gradient method to finally output the three-dimensional deformation rate at each geographic grid point. .

[0081] The defect in the above embodiment 1, and its technical contradiction, lies in: traditional Laplace smoothing ( The mathematical essence of norm regularization is to penalize large gradients in the model, forcing the generation of spatially smooth solutions. However, the specific application scenario of this invention—coal mine goaf—is not entirely smooth in its physical deformation pattern. Subsidence basins caused by coal mine goafs have clear, sharp physical boundaries where the deformation gradient is abrupt (i.e., discontinuous). Example 1 employs smoothing constraints (…). Regularization is used to solve a physical phenomenon that is inherently non-smooth / piecewise smooth. The regularization term over-penalizes true gradient abrupt changes, resulting in the goaf edges in the solution being blurred and smoothed, failing to accurately reflect the true situation of surface fractures and gradient change zones, thus affecting the accuracy of disaster assessment.

[0082] Example 2:

[0083] refer to Figure 1 To solve the problem in Example 1 To address the technical contradictions caused by edge blurring resulting from regularization, this embodiment modifies steps S41 and S42 as follows.

[0084] Modification of step S41:

[0085] The objective function in Example 1 Regularization term Set as an edge-preserving Normative total variational regularization term. Modified objective function. for:

[0086] ,

[0087] in, The L1 norm (Manhattan distance) is used. It is the L2 norm (Euclidean distance). For the calculation of L1 norm (). Calculate the square of the L2 norm. For the objective function Data items (InSAR data fidelity item + BeiDou constraint item):

[0088] ,

[0089] This is the total variation term. This embodiment preferably uses an anisotropic total variation, whose discrete form is defined as:

[0090] ,

[0091] in, For the deformation component index (E, N, U). This represents the total number of geographic grid points. For gradient operators, ; Traverse all grid points; and These are the first-order finite difference operators in the horizontal and vertical directions, respectively. Let c be the c-th component of the three-dimensional deformation vector X.

[0092] In order to solve the problem in (Example 1) The edge blurring issue caused by smoothing regularization was addressed by adopting... The total variation of the norm as a space regularization term The core theoretical basis lies in the fact that the physical phenomenon of coal mine goaf subsidence has the characteristic of being segmented and smooth: the deformation is smooth (small gradient) inside the subsidence basin and in the stable zone far away from the goaf; however, at the edge of the subsidence basin, there are abrupt changes in the deformation gradient (discontinuity) caused by surface fracturing or severe bending.

[0093] This embodiment uses -TV Regularization Its goal is to minimize the absolute value of the gradient (the rate of change of deformation). The key property of the norm is that it promotes sparsity. When applied to gradients, it tends to produce a gradient-sparse solution: it allows for large gradient values ​​at a few locations (i.e., the actual edges of the goaf) (because...). The penalty for large values ​​is much smaller. Simultaneously, it strongly compresses the gradients of other regions (smooth regions) to zero, thus effectively suppressing noise. Therefore, The -TV regularization term, as an edge-preserving regularization term, perfectly matches the piecewise smoothing physical characteristics of goaf areas. It can effectively smooth noise while preserving the gradient discontinuities at the edges of subsidence basins with high fidelity, significantly improving the ability to reconstruct goaf boundaries.

[0094] Modification of step S42 (new solver):

[0095] Due to the objective function Includes non-differentiable Due to the norm term, the conjugate gradient method in Example 1 fails. This example introduces the Alternating Direction Multiplier Method (ADMM) for iterative solution. The implementation details of the ADMM solver are as follows:

[0096] ① Problem Refactoring and Variable Splitting:

[0097] Introducing auxiliary variables and order The objective function is reconstructed into an equivalent constrained optimization problem:

[0098] ,

[0099] in, For L1 regularization terms related to Z in ADMM, .

[0100] ② Construct the augmented Lagrangian function:

[0101] Introducing dual variables and penalty parameters :

[0102] ,

[0103] in, To construct the augmented Lagrange function, dual variables transpose, The parameter is the augmented Lagrange Grange penalty parameter for the ADMM algorithm.

[0104] ③ADMM Iteration Steps:

[0105] In the In each iteration, updates are performed alternately. , , :

[0106] (a) Solving the X-subproblem Least squares):

[0107] ,

[0108] in, This is the index of the number of iterations in the ADMM algorithm. For the first The three-dimensional deformation vector of the next iteration. For the first Auxiliary variables for the next iteration ( ), For the first The dual variable (multiplier) of the next iteration is a smooth convex problem, which can be solved efficiently by the conjugate gradient method (CG) mentioned in Example 1.

[0109] (b) Solving the Z-subproblem Proximal operators):

[0110] ,

[0111] in, For the first The auxiliary variable for the next iteration gives this subproblem a closed-form solution, namely the soft-threshold shrinkage operator. The closed-form solution to the Z-subproblem is:

[0112] ,

[0113] in, The input vector for the soft thresholding operator is... Soft threshold operator Defined as:

[0114]

[0115] in, The threshold of the soft threshold operator ( ), The j-th element of vector v ).

[0116] (c) Solving the Y-subproblem (updating dual variables):

[0117]

[0118] in, For the first The dual variable of the next iteration.

[0119] Repeat the iteration until convergence.

[0120] The superiority of Example 2 over Example 1 stems from the fact that both use fundamentally different mathematical norms when constructing the regularization constraint terms, thus resolving the core technical contradiction between the mathematical model assumptions and the actual physical deformation present in Example 1. Example 1 uses... Norm regularization (such as Laplace smoothing) is mathematically essentially a penalty for the square of the model gradient. This penalty mechanism imposes extremely high weights on large gradient values ​​(i.e., locations of abrupt deformation changes), therefore Regularization strongly tends to produce a globally smooth solution, blurring real, sharp gradient abrupt changes (such as goaf edges) at all costs to minimize the objective function. Regularization term. However, surface subsidence in coal mine goaf areas is physically a typical piecewise smooth phenomenon: the interior and exterior regions of the basin are relatively smooth, but there are physical gradient abrupt changes at the subsidence edges. Therefore, Example 1 The smoothing assumption contradicts this physical reality, leading to inevitable edge blurring and boundary distortion. In contrast, Example 1 employs... Norm total variation (regularization) Its mathematical essence is to penalize the absolute value of the model gradient. The norm has a well-known sparsity-enhancing property; when applied to the gradient (i.e., TV), it tends to produce a gradient-sparse solution. This means the model allows for large gradient values ​​in a few locations (i.e., the sharp edges of the goaf), while forcing the gradient close to zero in most areas (the flat interior or exterior of the basin). This piecewise smooth mathematical property perfectly matches the physical characteristics of goaf subsidence, therefore... -TV regularization can effectively distinguish between real edge abrupt changes and random noise, suppressing noise while maintaining the sharp edges of the subsidence basin, thus solving the problem of... The technical contradiction of edge blurring caused by regularization. To achieve this theoretically superior... The model, Example 1, also features corresponding theoretical innovations in the solver. Because... The norm is not differentiable at zero, causing the efficient conjugate gradient (CG) method in Example 1 to fail. The alternating direction multiplier method (ADMM) introduced in Example 1 is a solution to this problem. - A standard and efficient theoretical framework for mixed norm optimization problems. ADMM cleverly decomposes this complex, non-differentiable problem into a single variable using the variable splitting technique. The least squares subproblem (which can still be solved efficiently using CG) and a The proximal operator subproblem. The key to the latter is that it has a closed analytical solution in the form of a soft-threshold contraction operator, which is precisely what is implemented. The ADMM algorithm is a core mathematical tool for studying norm sparsity. Therefore, it not only theoretically guarantees the sparsity of non-differentiable features... - Efficient solution to the TV objective function, and computationally efficient implementation of edge preservation (through soft thresholding shrinkage) and stable inversion (through... The approach to solving subproblems.

[0121] Example 3:

[0122] Based on Example 1 or Example 2, this embodiment adds an atmospheric delay correction step before step S11 (InSAR processing).

[0123] InSAR interferometric phase Contains multiple signals .in, The phase caused by deformation, Phase caused by terrain For noise phase, Phase delay caused by atmospheric delay is one of the main error sources in InSAR deformation monitoring. BeiDou / GNSS data can provide high-precision, high-temporal-resolution station atmospheric delay data, which can be used to accurately correct InSAR data.

[0124] Implementation steps:

[0125] ① Obtain BeiDou ZTD: Using the BeiDou / GNSS observation data obtained in step S12, calculate the total zenith delay (ZTD) for each station and each SAR imaging time. ZTD is divided into dry delay (ZHD) and wet delay (ZWD).

[0126] ②ZHD Calculation: ZHD changes slowly in space and time, and can be accurately calculated using the Saastamoinen model:

[0127]

[0128] in, The standard atmospheric pressure of the station. The geographical latitude of the station. This represents the elevation of the station.

[0129] ③ZWD Calculation: ZWD varies drastically in space and time. Using the BeiDou PPP calculation mode or the GAMIT network calculation mode, ZWD is estimated epoch-by-epoch as a random walk parameter.

[0130] (a) Spatial Interpolation and Correction: Sparse GNSS station ZWD data is interpolated onto a dense grid of InSAR imagery. The interpolated ZTD (ZHD+ZWD) is projected onto the LOS direction of InSAR and converted into atmospheric correction phase calculated by GNSS. and from the original interference phase Deducted from the middle.

[0131] (b) Post-processing: Perform the time-series InSAR solution in step S1 using the atmospherically corrected interferogram to obtain higher precision. and .

[0132] (c) Beneficial effects: This embodiment removes the main noise source (atmospheric delay) in InSAR observations by using BeiDou data, which significantly improves the input accuracy of the data fidelity item in step S41, thereby making the final calculated three-dimensional deformation field more reliable and accurate.

[0133] Example 4:

[0134] Based on Example 1 or Example 2, this embodiment adds step S5 after step S4, which outputs and verifies the output and verifies the three-dimensional deformation rate of each geographic grid point.

[0135] After calculation, the three-dimensional deformation field of the entire coal mine goaf was obtained. Verification showed that the calculated vertical deformation agreed well with the leveling measurement results, with a correlation coefficient of 0.98. The root mean square error (RMSE) of the east-west deformation, compared with the results from BeiDou stations, was ±1.2 mm / year. The accuracy of the north-south deformation was also significantly better than the traditional unconstrained method, with an RMSE of ±2.1 mm / year. The main drawback of this method is that the edges of the goaf are blurred.

[0136] In step S5, the output and verification of the three-dimensional deformation field, in order to objectively evaluate the accuracy of the solution results, it is necessary to use independent data that was not involved in the modeling for verification.

[0137] In practice, cross-validation can be used:

[0138] In step S1, data from 15 BeiDou monitoring stations were acquired. In step S4, when building the model, data from all 15 stations is not used for constraints. Instead, a portion of the stations (e.g., 3 stations) are reserved as independent verification points, and only the remaining 12 stations are used as constraint terms. The vector is used to solve the three-dimensional deformation field.

[0139] After the model is solved (step S4), the calculated three-dimensional deformation field is extracted. The deformation rate values ​​at these three independent verification points.

[0140] Will The values ​​of the three-dimensional deformation field calculated by the model and the actual observation values ​​from these three BeiDou stations are compared. (i.e., independent BeiDou data used for verification) are compared point-by-point and component-by-component (E, N, U), and accuracy statistics such as root mean square error (RMSE) are calculated to serve as an objective evaluation of the out-of-model conformity accuracy.

[0141] Leveling is the most classic and accurate method for obtaining elevation (vertical) deformation. If a leveling route is set up in the monitoring area and simultaneous observations are conducted, then these leveling data (providing only high-precision vertical deformation) will yield accurate data. (This is to verify the solution of the present invention) The best independent data source for components.

[0142] During verification, the vertical deformation rate calculated by this invention is used. Extraction is performed along the leveling route, and the vertical rate is compared with that obtained from leveling measurements (e.g., drawing profiles, calculating correlation coefficients and RMSE) to evaluate the accuracy of the vertical deformation component solution.

[0143] Based on the same technical concept as the aforementioned embodiments, this invention also provides a BeiDou-constrained InSAR three-dimensional deformation field calculation system for coal mine goaf areas. This system is implemented using the aforementioned BeiDou-constrained InSAR three-dimensional deformation field calculation method for coal mine goaf areas, comprising: a LOS-direction deformation rate acquisition module, used to acquire and preprocess the LOS-direction deformation rates of the ascending and descending orbits of the coal mine goaf area; and an InSAR three-dimensional deformation rate acquisition module, used to construct geometric relationship equations between the LOS-direction deformation rates of the ascending and descending orbits and the deformation rates in the east-west, north-south, and vertical directions based on the preprocessing results, to obtain the InSAR three-dimensional deformation field. The system includes: a 3D deformation rate acquisition module (using BeiDou data to calculate the 3D deformation rate of each BeiDou monitoring station within the coal mine goaf); an objective function construction module (using InSAR data fidelity terms, BeiDou data constraint terms, and edge-preserving regularization terms for stable calculations); and a comprehensive 3D deformation rate calculation module (using alternating direction multiplier method to solve the objective function based on the InSAR 3D deformation rate and the 3D deformation rate of each BeiDou station to obtain the 3D deformation rate of each geographic grid point).

[0144] Based on the same technical concept as the foregoing embodiments, the present invention also provides an electronic device, including a memory and a processor, wherein the memory stores program instructions that are executed by the processor, and the processor calls the program instructions to execute the aforementioned method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints in InSAR.

[0145] Based on the same technical concept as the foregoing embodiments, the present invention also provides a non-transitory computer-readable storage medium that stores computer instructions that cause the computer to execute the aforementioned method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints using InSAR.

[0146] In summary, the above embodiments first acquire the InSAR temporal LOS deformation rate and sparse BeiDou 3D deformation rate data; then, geocoding and resampling the InSAR temporal LOS deformation rate are performed to establish a geometric model of InSAR LOS observation and 3D deformation; finally, an objective function is constructed that includes an InSAR data fidelity term, a BeiDou constraint term, and a spatial regularization term, and the spatial regularization term is constrained to a total variational (norm) regularization term that preserves the edges; and accordingly, the objective function is iteratively solved using the Alternating Direction Multiplier Method (ADMM) and a soft threshold shrinkage operator to obtain the 3D deformation field.

[0147] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; 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 or all of the technical features therein; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints using InSAR, characterized in that, include: Obtain the LOS deformation rate of the lifting and lowering rails in the goaf of a coal mine and perform preprocessing. Based on the preprocessing results, geometric equations relating the LOS deformation rate of the ascending and descending orbits to the deformation rates in the east-west, north-south, and vertical directions are constructed to obtain the three-dimensional deformation rate of InSAR. Based on the observation data from each Beidou monitoring station in the coal mine goaf area, the three-dimensional deformation rate of each Beidou station was calculated. The objective function for solving the three-dimensional deformation field is constructed based on the data fidelity term of InSAR observations, the constraint term of BeiDou observations, and the edge-preserving regularization term for stable solution. The expression corresponding to the edge-preserving regularization term is: ， in, To preserve the regularization term for the edges, For deformation component index, This represents the total number of geographic grid points. Indices representing geographic grid points For gradient operators, It is an L1 norm. Three-dimensional deformation vector The c-th component, and These are first-order difference operators along the horizontal and vertical directions, respectively; Based on the three-dimensional deformation rate of InSAR and the three-dimensional deformation rate of each BeiDou station, the objective function is solved using the alternating direction multiplier method to obtain the three-dimensional deformation rate of each geographic grid point.

2. The method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints according to claim 1, characterized in that, The expression for the geometric relation equation is: ， in, For LOS deformation rate, Let be the unit geometric vector in the direction of the satellite's line of sight. , , Let be the deformation rates in the east-west, north-south, and vertical directions, respectively.

3. The method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints as described in claim 1, characterized in that, The expression for the objective function is: ， in, Let be the objective function. For InSAR observation vectors, The geometric coefficient matrix, Let be the three-dimensional deformation vector to be determined. , , , Let be the deformation rates in the east-west, north-south, and vertical directions, respectively. For BeiDou observation vectors, For the sampling matrix of Beidou stations, This represents the overall weighting coefficient for the BeiDou constraint term. To preserve the regularization term for the edges, The weight coefficients for the edge-preserving regularization term.

4. The method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints as described in claim 1, characterized in that, The objective function is solved using the alternating direction multiplier method to obtain the three-dimensional deformation rate of each geographic grid point, including: The objective function is split into a function that includes a data fidelity term and a BeiDou constraint term. Norm subproblem, and a problem involving margin-preserving regularization. The norm subproblem can be solved using the conjugate gradient method. Norm subproblem and its solution using the soft-threshold shrinkage operator Norm subproblem; The three-dimensional deformation vector, auxiliary variables, and dual variables in the objective function are alternately and iteratively updated until convergence is achieved to obtain the three-dimensional deformation rate of each geographic grid point.

5. The method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints according to claim 1, characterized in that, The acquisition of the LOS deformation rate of the ascending and descending orbits includes: Tropospheric atmospheric delay correction was performed on the acquired SAR image data for both ascending and descending orbits using observation data from various BeiDou monitoring stations. Based on the corrected ascending and descending SAR image data, the LOS deformation rate of ascending and descending orbits was obtained by time-series InSAR technology.

6. The method for calculating the three-dimensional deformation field of coal mine goaf based on BeiDou constraints according to claim 1, characterized in that, The solution method also includes: Step S5: Output and verify the three-dimensional deformation rate of each geographic grid point.

7. A BeiDou-constrained InSAR system for calculating the three-dimensional deformation field of coal mine goaf, characterized in that, The method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints, as described in any one of claims 1 to 6, includes: The LOS deformation rate acquisition module is used to acquire the LOS deformation rate of the rising and falling rails in the goaf of the coal mine and perform preprocessing. The InSAR three-dimensional deformation rate acquisition module is used to construct geometric relationship equations between the LOS deformation rate of the ascending and descending orbits and the deformation rates in the east-west, north-south, and vertical directions based on the preprocessing results, so as to obtain the three-dimensional deformation rate of InSAR. The BeiDou 3D deformation rate acquisition module is used to calculate the 3D deformation rate of each BeiDou monitoring station based on the observation data of each BeiDou monitoring station in the coal mine goaf area. The objective function construction module is used to construct the objective function for solving the three-dimensional deformation field based on the data fidelity term based on InSAR observations, the constraint term based on BeiDou observations, and the edge-preserving regularization term for stable solution. The integrated three-dimensional deformation rate calculation module is used to solve the objective function based on the three-dimensional deformation rate of the InSAR and the three-dimensional deformation rate of each BeiDou station, and to obtain the three-dimensional deformation rate of each geographic grid point.

8. An electronic device, characterized in that, The system includes a memory and a processor. The memory stores program instructions that are executed by the processor. The processor invokes the program instructions to execute the method for calculating the three-dimensional deformation field of a coal mine goaf based on BeiDou constraints, as described in any one of claims 1 to 6.

9. A non-transitory computer-readable storage medium, characterized in that, The non-transitory computer-readable storage medium stores computer instructions that cause the computer to execute the InSAR three-dimensional deformation field calculation method for coal mine goaf based on Beidou constraints as described in any one of claims 1 to 6.

Images