Mining area earth surface three-dimensional displacement resolving method based on Farneback optical flow method and ascending and descending rail SAR (Synthetic Aperture Radar) collaborative observation

Abstract

The invention discloses a mining area earth surface three-dimensional displacement resolving method based on Farneback optical flow method and lifting rail SAR collaborative observation, and belongs to the technical field of synthetic aperture radar remote sensing monitoring. The method comprises the following steps: firstly, performing Farneback optical flow displacement field modeling on a time sequence SAR intensity image based on gray conservation hypothesis and second-order polynomial parameterization expansion, and extracting a dense two-dimensional displacement field with sub-pixel precision through a multi-scale pyramid layered structure and iterative weighted least square optimization; secondly, constructing a projection equation of earth surface three-dimensional displacement and radar sight line direction and azimuth direction displacement based on strict side-view imaging geometry in combination with lifting rail multi-view-angle SAR observation data, and integrating double-view-angle data to form an overdetermined observation equation set; and finally, solving the overdetermined equation set by adopting a least square criterion, and introducing a robust estimation strategy based on singular value decomposition to carry out robust inversion when the equation set is ill-conditioned or has a gross error, thereby solving high-precision three-dimensional displacement fields in the vertical, east-west and north-south directions of the earth surface.

Description

Technical Field

[0001] This invention belongs to the field of synthetic aperture radar remote sensing monitoring and deformation inversion technology, and in particular relates to a method for calculating the three-dimensional displacement of the surface of a mining area based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation. Background Technology

[0002] Mining-induced surface deformation is characterized by its wide spatial range, rapid evolution, and large deformation gradient. High-precision three-dimensional displacement monitoring is crucial for understanding mining patterns and mitigating disaster risks. In recent years, Interferometric Synthetic Aperture Radar (InSAR) technology, with its high spatiotemporal resolution and wide coverage, has become one of the core methods for monitoring deformation in mining areas. However, this technology heavily relies on the stability of the radar phase. In areas with intense mining and large deformation gradients, radar signals are prone to severe degradation of interferometric phase quality due to spatiotemporal decoherence, greatly limiting the effective application of traditional InSAR in such scenarios.

[0003] To overcome the phase incoherence problem, Pixel Offset Tracking (POT) technology, based on image registration principles, inverts deformation by calculating sub-pixel-level spatial displacement between SAR images, demonstrating unique advantages in capturing large-gradient deformations ranging from centimeters to meters. The further developed Adaptive Cross-Correlation Window (ACC-POT) algorithm optimizes monitoring performance through dynamic window adjustment. Nevertheless, existing POT-like methods still suffer from two significant drawbacks: first, their displacement estimation accuracy based on intensity cross-correlation is generally lower than that of InSAR technology based on phase measurement, especially when image quality is poor; second, the computational complexity of cross-correlation search increases dramatically with image resolution and window size, limiting its efficiency for engineering applications in high-resolution, large-scale monitoring scenarios.

[0004] To accurately reconstruct the three-dimensional displacement field of the Earth's surface, various technical approaches have been proposed, such as combining multi-track / multi-platform SAR imagery and fusing InSAR and POT observations. Among these, the method of obtaining multi-line deformation by combining ascending / descending orbit or different-view SAR data, and then inverting the three-dimensional displacement, is theoretically feasible. However, in practice, it is often limited by data fusion challenges caused by inconsistencies in resolution and revisit periods between different SAR platforms. Furthermore, the inherent accuracy limitations of POT technology directly affect the reliability of the final three-dimensional inversion results. In addition, observations from a single SAR orbit can only provide two displacement components: the line-of-sight and the azimuth. This results in inherent limitations and ill-posed problems in the three-dimensional deformation inversion process, which urgently need to be overcome through multi-view collaborative observations.

[0005] China's independently developed LuTan-1 dual-satellite system, with its L-band anti-decoherence characteristics, high spatial resolution, and quasi-synchronous ascent-descent orbit collaborative observation capabilities, provides an ideal data platform for solving the aforementioned challenges of multi-view data fusion and 3D inversion. Meanwhile, the Farneback optical flow method, originating from the field of computer vision, based on the gray-level conservation assumption and parameterized motion field modeling, possesses the ability to estimate sub-pixel-level dense displacement fields with high accuracy and efficiency, providing a new algorithmic option for improving the accuracy and computational efficiency of extracting 2D displacement fields from SAR intensity images. However, how to combine the efficient displacement extraction capability of the Farneback optical flow algorithm with the unique collaborative observation advantages of the LuTan-1 satellite to develop a method for solving large-gradient 3D displacement fields in mining areas that balances accuracy, efficiency, and robustness remains a technical challenge that needs to be overcome. Currently, no research has systematically integrated the Farneback optical flow algorithm with LuTan-1 quasi-synchronous ascent-descent orbit SAR images to achieve a complete solution for high-precision, high-efficiency 3D displacement monitoring of the mining area surface. Summary of the Invention

[0006] This invention aims to overcome the aforementioned defects in the prior art. It addresses the problems of traditional InSAR being prone to incoherence in monitoring large gradients and three-dimensional deformation in mining areas, the difficulty in balancing accuracy and efficiency of traditional POT methods, and the ill-conditioned nature of three-dimensional inversion from a single SAR perspective. The invention provides a solution method that can achieve fast and robust acquisition of three-dimensional displacement fields with full vector and high spatiotemporal resolution.

[0007] To solve the above-mentioned technical problems, the present invention is achieved through the following technical solution:

[0008] This invention relates to a method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow and combined observations of rising and falling orbit SAR, comprising the following steps:

[0009] S1: Based on the gray-level conservation assumption and second-order polynomial parameterization expansion, Farneback optical flow displacement field modeling is performed on time-series SAR intensity images, and the surface dense two-dimensional displacement field with sub-pixel accuracy is extracted through multi-scale pyramid hierarchical structure and iterative weighted least squares optimization.

[0010] S2: Combine multi-view SAR observation data from both ascending and descending orbits to construct a rigorous geometric model for calculating the three-dimensional deformation of the surface in the mining area: Based on the geometric projection relationship between the line-of-sight and azimuth observation displacements of a single-track SAR and the actual three-dimensional displacement of the surface, establish mathematical projection equations; integrate dual-view observation data from ascending and descending orbits to form an overdetermined set of observation equations.

[0011] S3: Perform high-precision robust inversion on the overdetermined observation equations constructed in step S2: use the least squares criterion to solve the three-dimensional displacement components; when the equations have ill-conditioned properties or observation gross errors, further use a robust estimation strategy based on singular value decomposition to solve them, and finally obtain the three-dimensional displacement fields of the Earth's surface in the vertical, east-west and north-south directions.

[0012] As a preferred embodiment of the present invention, step S1 specifically includes:

[0013] S1.1: Perform second-order polynomial parameterization modeling on the local neighborhood grayscale distribution of each pixel in the SAR intensity image, and represent the pixel grayscale value as a quadratic function of its spatial coordinates;

[0014] S1.2: Constructing local polynomial coefficients representing the same surface region between adjacent temporal images. , , A set of equations relating the changes to the pixel displacement vector;

[0015] S1.3: Using the multi-scale pyramid hierarchical structure and iterative weighted least squares optimization, the dense two-dimensional displacement field with sub-pixel accuracy is obtained.

[0016] As a preferred embodiment of the present invention, in step S1.1, a set of matrices containing symmetric matrices is used. gradient vector and constant term The parameters fully characterize the structure and intensity features of the local image. The specific second-order polynomial parameterization model is as follows:

[0017]

[0018] In the formula, Represents the pixel grayscale value. Represents the two-dimensional coordinates of a pixel The matrix formed, symbol Represents transpose, symmetric matrix Characterizing local curvature, gradient vector Reflecting the spatial gradient distribution, constant term Indicates the baseline grayscale value.

[0019] Where, symmetric matrix for:

[0020]

[0021] gradient vector for: .

[0022] As a preferred technical solution of the present invention, step S1.3 specifically involves: introducing a multi-scale pyramid processing framework from coarse to fine; at each scale level, constructing and minimizing a weighted least squares objective function that considers neighborhood consistency; iteratively optimizing the displacement field of that layer; and finally accumulating the results of each layer to obtain a dense two-dimensional displacement field with sub-pixel accuracy at the original image scale. Furthermore, when solving the displacement field, a filtering window is set in conjunction with neighborhood information for optimization, and Gaussian weights or uniform weights are used to smooth the displacement field.

[0023] As a preferred embodiment of the present invention, step S2 specifically includes:

[0024] S2.1: Based on the geometric principles of SAR satellite side-looking imaging, strictly establish the three-dimensional displacement vector of the surface point in the northeast-central geographic coordinate system. The three-dimensional displacement vector is the mathematical projection equation between its two observable components in the single-track SAR image coordinate system. Including vertical direction East-west North-South Direction The two observable components are the radar line-of-sight displacement and the line-of-sight displacement. and satellite flight azimuth displacement The equation uses the satellite incident angle θ and heading angle α as the core geometric parameters;

[0025] S2.2: By combining the ascending and descending SAR observation data, the projection equations from two independent geometric perspectives are combined to construct a system based on the vertical direction. East-west North-South Direction The overdetermined observation equations are for unknowns, thus providing sufficient mathematical constraints for the effective separation of three-dimensional displacement components;

[0026] S2.3: Specific imaging geometric parameters based on ascending and descending orbit satellite imagery, i.e., incident angle. , and heading angle , Determine the coefficient matrix of the overdetermined system of equations. Furthermore, the geometric observability of the three-dimensional displacement solution is evaluated by analyzing the condition number or singular value distribution of the coefficient matrix B.

[0027] As a preferred embodiment of the present invention, in step S2.2, the overdetermined observation equations are expressed as follows:

[0028]

[0029] Among them, the observation vector , including the observed displacement components of ascending and descending orbits;

[0030] Vector to be determined ;

[0031] coefficient matrix It is determined by the incident angle θ and heading angle α of the satellite in ascending or descending orbit.

[0032] As a preferred embodiment of the present invention, the coefficient matrix The expression is:

[0033] .

[0034] As a preferred embodiment of the present invention, in step S3:

[0035] If the matrix is ​​nonsingular, the estimated value of the three-dimensional displacement vector D obtained by using the least squares criterion is:

[0036] ;

[0037] When the coefficient matrix of the overdetermined observation equation system is ill-conditioned or the observations contain gross errors, a robust estimation strategy based on singular value decomposition is adopted. By truncating small singular values, a generalized inverse matrix is ​​constructed, and a stable least squares solution is obtained under the minimum norm principle.

[0038] As a preferred embodiment of the present invention, the robust estimation strategy based on singular value decomposition in step S3 is specifically as follows:

[0039] Perform singular value decomposition on coefficient matrix B ,

[0040] In the formula, U represents an M×N orthogonal matrix, and its column vectors are... eigenvectors;

[0041] S represents an M×N diagonal matrix whose diagonal elements are The square root of the non-zero eigenvalue;

[0042] V represents an N×N orthogonal matrix whose column vectors are eigenvectors;

[0043] The generalized inverse matrix is ​​obtained by reasonably truncating the small singular values ​​in the singular value matrix S. And then through A robust estimate of the three-dimensional displacement vector is obtained.

[0044] The present invention has the following beneficial effects:

[0045] This invention innovatively introduces the Farneback optical flow method, derived from computer vision, into the monitoring of large gradient deformation in mining areas. It overcomes the problems of low computational efficiency and limited accuracy of traditional cross-correlation-based POT algorithms, and realizes the extraction of two-dimensional displacement fields from SAR images with high precision and high efficiency.

[0046] This invention fully utilizes the unique quasi-synchronous ascent and descent orbit collaborative observation advantage of the domestically produced LuTan-1 dual-satellite system. By constructing a multi-view geometric model, the three-dimensional deformation inversion is transformed from an underdetermined problem into an overdetermined problem, effectively solving the inherent problem that single-orbit SAR observations cannot separate three-dimensional displacement components.

[0047] The method proposed in this invention, while ensuring monitoring accuracy at the centimeter to sub-meter level, significantly reduces the calculation time from several hours in traditional methods to seconds, thereby significantly improving the timeliness and engineering application potential of three-dimensional deformation monitoring in mining areas.

[0048] This method provides a reliable technical approach for obtaining full-vector, high spatiotemporal resolution three-dimensional deformation fields of the surface in mining areas. It not only helps to deepen the understanding of mining subsidence mechanisms, but also has important value for the accurate identification and early warning of geological disasters such as landslides and ground fissures. At the same time, it demonstrates the huge application potential of the domestically produced LT-1 satellite in the field of high-precision engineering deformation monitoring.

[0049] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

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

[0051] Figure 1 This is a flowchart illustrating the method in this invention;

[0052] Figure 2 This is a flowchart illustrating step S1 in this invention;

[0053] Figure 3 This is a flowchart illustrating step S2 in this invention;

[0054] Figure 4 This is a flowchart illustrating step S3 in this invention. Detailed Implementation

[0055] 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, and 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.

[0056] Please see Figure 1-4 As shown, the present invention is a method for calculating the three-dimensional displacement of the surface of a mining area based on the Farneback optical flow method and the joint observation of the rising and falling orbit SAR. S1: Based on the gray level conservation assumption and the second-order polynomial parameterization expansion, the Farneback optical flow displacement field is modeled on the time-series SAR intensity image, and the surface dense two-dimensional displacement field with sub-pixel accuracy is extracted through the multi-scale pyramid hierarchical structure and iterative weighted least squares optimization.

[0057] S2: Combine multi-view SAR observation data from both ascending and descending orbits to construct a rigorous geometric model for calculating the three-dimensional deformation of the surface in the mining area: Based on the geometric projection relationship between the line-of-sight and azimuth observation displacements of a single-track SAR and the actual three-dimensional displacement of the surface, establish mathematical projection equations; integrate dual-view observation data from ascending and descending orbits to form an overdetermined set of observation equations.

[0058] S3: Perform high-precision robust inversion on the overdetermined observation equations constructed in step S2: use the least squares criterion to solve the three-dimensional displacement components; when the equations have ill-conditioned properties or observation gross errors, further use a robust estimation strategy based on singular value decomposition to solve them, and finally obtain the three-dimensional displacement fields of the Earth's surface in the vertical, east-west and north-south directions.

[0059] One specific application of this embodiment is:

[0060] S1: Based on the Farneback optical flow method, this algorithm extracts the two-dimensional displacement field of the Earth's surface from temporal SAR intensity images with high precision and efficiency. Its core principle is to mathematically model the local gray-level distribution of the image using a second-order polynomial. By analyzing the changes in polynomial coefficients between adjacent temporal images, a linear equation is established between the coefficient changes and pixel displacements, thus transforming the complex image matching problem into an efficient mathematical solution problem. To handle large displacements and ensure stability, the algorithm employs a multi-scale pyramid framework from coarse to fine for iterative optimization. Ultimately, the output is a dense two-dimensional displacement field (including range and azimuth displacements) with sub-pixel precision, laying a reliable data foundation for subsequent three-dimensional calculations.

[0061] The specific implementation steps are as follows:

[0062] S1.1: Establish a second-order polynomial parameterized representation model of the local gray-level distribution of SAR intensity images.

[0063] The Farneback optical flow model is based on a second-order polynomial expansion model, which parameterizes the local neighborhood gray-level distribution of each pixel in the image.

[0064]

[0065] In the formula, Represents the pixel grayscale value; Represents the two-dimensional coordinates of a pixel The matrix formed; symbols Indicates transpose; symmetric matrix Characterizing local curvature; gradient vector Reflects spatial gradient distribution; constant term This represents the baseline grayscale value.

[0066] Where, symmetric matrix for:

[0067]

[0068] gradient vector for:

[0069]

[0070] This step primarily involves using the fundamental assumption of gray-level conservation to perform second-order polynomial parameterization modeling of the local neighborhood gray-level distribution of each pixel in the SAR intensity image. This represents the pixel gray-level value as a quadratic function of its spatial coordinates, using a set of symmetric matrices. gradient vector and constant term The parameters fully characterize the structure and intensity features of the local image.

[0071] S1.2: Construct a set of equations relating polynomial coefficient transformation and displacement vector between adjacent temporal images.

[0072] If we assume the original pixel position is as follows:

[0073]

[0074] The position of the pixel after movement is:

[0075]

[0076] in ; ; Then, based on the polynomial coefficient transformation relationship, a system of linear equations about the displacement vector is constructed. If If it is not singular, then we can obtain:

[0077]

[0078] According to theoretical derivation, there must be one. However, this requirement may not be met in practice, so an average approximation is used. Let:

[0079]

[0080] but:

[0081]

[0082] This step mainly involves theoretically deriving the variation law of local polynomial coefficients representing the same surface area between adjacent temporal images, and establishing the relationship between displacement vectors and polynomial coefficients. The precise mathematical relationship between the changes transforms the original displacement problem based on pixel grayscale matching into a problem of solving a system of linear equations constructed by parameter changes.

[0083] S1.3: A robust solution to the displacement field is achieved by using a multi-scale pyramid hierarchical structure and iterative weighted least squares optimization.

[0084] Considering the consistency within the neighborhood, a weighted least squares method is used to construct the objective function for optimization to obtain the displacement:

[0085]

[0086] The pyramid layering process is as follows: a Gaussian pyramid is built for each image scene. The upper pyramid is obtained by downsampling the lower pyramid, and the bottom pyramid is the original image. Let the target displacement in the original image be... Then the first The displacement of the layer is:

[0087]

[0088] in, This represents the number of layers in the pyramid image. The optical flow result of the top layer is reflected in the next lower layer and used as the optical flow estimate for that layer. Then the optical flow of the second-to-last layer can be expressed as:

[0089]

[0090] In each image layer, by minimizing the objective function By obtaining the optical flow of this layer and iterating sequentially, the target displacement of the bottom layer of the pyramid can be obtained. :

[0091]

[0092] By accumulating the optical flow results at different levels of the pyramid in a coarse-to-fine manner, a stable estimate of the target pixel displacement is achieved. Furthermore, to obtain a smooth and reliable optical flow field, the Farneback optical flow algorithm incorporates neighborhood information to set a certain filtering window for optimization and uses Gaussian or uniform weights for smoothing.

[0093] This step introduces a multi-scale pyramid processing framework from coarse to fine. At each scale level, the displacement field of that layer is iteratively optimized by constructing and minimizing a weighted least squares objective function that considers neighborhood consistency. Finally, the results of each layer are accumulated to obtain a dense two-dimensional displacement field with sub-pixel accuracy at the original image scale.

[0094] S2: Construct a rigorous geometric model to convert the two-dimensional observation displacements (line of sight and azimuth) extracted in S1 into actual three-dimensional surface deformation. The specific process is as follows: First (S2.1), establish the observation geometric equations for a single satellite, clarifying how the vertical, east-west, and north-south displacements of the Earth's surface are projected onto the satellite's line of sight and flight direction. Then (S2.2), addressing the "ill-conditioned" problem of solving three-dimensional vectors from a single perspective, combine the ascending and descending orbit data acquired quasi-synchronously by the LuTan-1 satellite. Simultaneously establish the two sets of observation equations to construct an overdetermined system of equations with more equations than unknowns, providing sufficient mathematical constraints for the accurate separation of the three-dimensional components. Finally (S2.3), determine the coefficient matrix of the equation system based on specific satellite parameters and evaluate its solvability.

[0095] The specific implementation steps are as follows:

[0096] S2.1: Derive the rigorous geometric projection relationship between the observed displacements in the line of sight and azimuth of a single-track SAR and the actual three-dimensional displacements of the Earth's surface.

[0097] In a simplified model that only considers vertical settlement:

[0098]

[0099] In the formula Indicates vertical deformation. Indicates deformation along the line of sight. This represents the incident angle of the corresponding pixel. However, surface deformation caused by coal mining usually includes significant vertical subsidence and horizontal movement, so this simplified model cannot accurately reflect the complex three-dimensional characteristics of actual deformation in the mining area in a strict sense.

[0100] To reveal the projection mechanism of surface deformation in the mining area along the satellite line of sight, it is necessary to establish a spatial projection model between satellite imaging geometry and three-dimensional surface displacement.

[0101]

[0102] In this model, , , These represent the vertical, east-west, and north-south directions in the geographic coordinate system, respectively. , These represent the deformation in the east-west direction and the deformation in the north-south direction, respectively. The heading angle is defined as the angle between the direction of true north and the direction of the satellite's flight, rotated clockwise. This represents the projection of the deformation of a ground target point onto the azimuth direction. However, traditional single-satellite observations are constrained within a two-dimensional plane formed by the line of sight and the azimuth direction, resulting in inherent ill-conditioning problems in three-dimensional deformation field inversion.

[0103] This step mainly involves rigorously establishing the three-dimensional displacement vector (vertical direction) of the surface point in the northeast-central geographic coordinate system based on the geometric principles of SAR satellite side-looking imaging. East-west North-South Direction ) and its two observable components in the single-track SAR image coordinate system (radar line-of-sight displacement) Satellite flight azimuth displacement The mathematical projection equation between the satellite and the incident angle explicitly includes the satellite's angle of incidence. With heading angle As a core geometric parameter.

[0104] S2.2: Integrate dual-view observation data from the ascending and descending orbits to construct an overdetermined set of observation equations to overcome the ill-conditioned nature of the inversion.

[0105] To overcome this dimensional limitation, a multi-platform collaborative observation system needs to be constructed. This involves fusing observational data acquired by two or more satellites with heterogeneous orbits (i.e., different imaging geometries) during quasi-synchronous periods, thereby providing sufficient observational constraints for the accurate separation of three-dimensional deformation components. The LuTan-1 satellite, with its unprecedented ascending and descending orbit collaborative observation capabilities, provides an ideal solution to the aforementioned three-dimensional deformation inversion problem. Operating in the more penetrating L-band, the system possesses strong anti-decoherence capabilities, achieving a maximum spatial resolution of 1.6 meters (line-of-sight) × 1.9 meters (azimuth) in strip 1 mode, enabling clear identification of subtle deformation features of targets such as mining structures and landslide cracks. Crucially, the LuTan-1 dual-satellite operation in collaborative observation mode compresses the revisit interval of the same region to less than 8 days, significantly reducing deformation signal confusion caused by time asynchrony and providing a reliable data foundation for quasi-synchronous, high-precision three-dimensional deformation field reconstruction.

[0106] Therefore, based on the provided functional relationship between the line of sight and azimuth and the three-dimensional deformation in the geographic coordinate system, the functional model of the three-dimensional deformation observation equation is constructed as follows:

[0107]

[0108] This transforms the 3-D deformation solution from solving an underdetermined system of equations to solving an overdetermined system of equations. and The deformations in the LOS and AZI directions, respectively, are obtained from the ascending SAR satellite. and These are the LOS and AZI deformations obtained from a descending SAR satellite, respectively. , These are the incident angles for satellites in ascending and descending orbits, respectively. , These are the heading angles of the satellites in ascending and descending orbits, respectively.

[0109] This step primarily addresses the inherent limitation of single SAR orbit observations in directly separating three-dimensional displacements. By combining the ascending and descending orbit observation data acquired quasi-synchronously by the LuTan-1 system, the projection equations from the two independent geometric perspectives are combined to form a system containing four (or preferably three) independent equations, solving for only three unknowns. , , The overdetermined equations provide sufficient mathematical constraints for the effective separation of three-dimensional displacement components.

[0110] S2.3: Determine the coefficient matrix of the overdetermined system of equations and analyze the geometric observability of the three-dimensional displacement solution.

[0111] The simplified expression of the observation equations is as follows:

[0112]

[0113] in, The geometric transformation coefficient matrix can be represented as:

[0114]

[0115] This step is primarily based on the specific imaging geometry parameters (incident angle) of the ascending and descending satellite imagery. , and heading angle , The coefficient matrix of the overdetermined equation system is calculated precisely. The behavior of this matrix directly determines the stability and accuracy of the three-dimensional displacement inversion. Analyzing its condition number or singular value distribution can evaluate the separability of different displacement components and the rationality of the solution scheme.

[0116] S3 aims to accurately and robustly derive the complete three-dimensional displacement field of the Earth's surface from the overdetermined equations constructed in S2. First (S3.1), the overdetermined equations are solved using the classical least squares method, ideally obtaining the best estimate of the three-dimensional displacements (vertical, east-west, and north-south), effectively suppressing random observation errors. However, when poor observation geometry (e.g., close viewing angles from ascending and descending orbit satellites) leads to an "ill-conditioned" coefficient matrix, or when gross errors exist in the data, the least squares solution may be unstable. Therefore, step (S3.2) introduces a robust estimation strategy based on singular value decomposition (SVD). By analyzing and reasonably truncating small singular values, a generalized inverse matrix is ​​constructed, thereby obtaining a stable and reliable solution insensitive to errors under the minimum norm principle. Finally, the output is a three-dimensional displacement field with centimeter to sub-meter accuracy, providing direct evidence for deformation analysis and disaster early warning in mining areas.

[0117] The specific implementation steps are as follows:

[0118] S3.1: The best linear unbiased estimate of the three-dimensional displacement components is obtained by using the least squares criterion.

[0119] If matrix If the surface is non-singular, then the best estimate of 3-D surface deformation in the geographic coordinate system can be obtained using the least squares method:

[0120]

[0121] This step involves applying the classical least squares estimation criterion to the overdetermined observation equations constructed in step S2, under the ideal condition that the observation noise satisfies the Gaussian white noise assumption, and directly obtaining the three-dimensional displacement vector by solving the normal equations. , , The optimal estimate in the least squares sense is a solution that can effectively suppress the influence of random observation errors.

[0122] S3.2: Introduce a robust solution strategy based on singular value decomposition for ill-conditioned problems.

[0123] When matrix If it is strange, then Since it is a singular matrix, the values ​​obtained using least squares are not unique. To solve this problem, its generalized inverse matrix can be derived by using the SVD method, and the least squares solution can be obtained under the minimum norm principle.

[0124]

[0125] In the formula, U represents an M×N orthogonal matrix. The eigenvectors constitute the elements of this matrix; S represents an M×N diagonal matrix. The eigenvalues ​​form the diagonal elements of the matrix; V represents an N×N orthogonal matrix. The eigenvectors constitute the elements of this matrix. Assuming the generalized inverse matrix of B is B+, the optimal estimate of 3-D surface deformation in the geographic coordinate system can be obtained using the SVD method:

[0126]

[0127] This step addresses the issue that when the observation geometry leads to ill-conditioned coefficient matrices (e.g., close incident angles on ascending and descending orbits) or gross errors exist in the observations, the least squares solution may be unstable or unreliable. To address this, singular value decomposition (SVD) is employed. By analyzing and appropriately truncating small singular values, a generalized inverse of the coefficient matrix is ​​constructed, thereby obtaining a stable least squares solution that is insensitive to observation errors under the minimum norm principle.

[0128] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0129] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to the specific implementations described. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation, characterized in that, Includes the following steps: S1: Based on the gray-level conservation assumption and second-order polynomial parameterization expansion, Farneback optical flow displacement field modeling is performed on time-series SAR intensity images, and the surface dense two-dimensional displacement field with sub-pixel accuracy is extracted through multi-scale pyramid hierarchical structure and iterative weighted least squares optimization. S2: Combine multi-view SAR observation data from both ascending and descending orbits to construct a rigorous geometric model for calculating the three-dimensional deformation of the surface in the mining area: Based on the geometric projection relationship between the line-of-sight and azimuth observation displacements of a single-track SAR and the actual three-dimensional displacement of the surface, establish mathematical projection equations; integrate dual-view observation data from ascending and descending orbits to form an overdetermined set of observation equations. S3: Perform high-precision robust inversion on the overdetermined observation equations constructed in step S2: use the least squares criterion to solve the three-dimensional displacement components; when the equations have ill-conditioned properties or observation gross errors, further use a robust estimation strategy based on singular value decomposition to solve them, and finally obtain the three-dimensional displacement fields of the Earth's surface in the vertical, east-west and north-south directions.

2. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 1, characterized in that, Step S1 specifically includes: S1.1: Perform second-order polynomial parameterization modeling on the local neighborhood grayscale distribution of each pixel in the SAR intensity image, and represent the pixel grayscale value as a quadratic function of its spatial coordinates; S1.2: Constructing local polynomial coefficients representing the same surface region between adjacent temporal images. , , A set of equations relating the changes to the pixel displacement vector; S1.3: Using the multi-scale pyramid hierarchical structure and iterative weighted least squares optimization, the dense two-dimensional displacement field with sub-pixel accuracy is obtained.

3. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 2, characterized in that, In step S1.1, a set of symmetric matrices is used. gradient vector and constant term The parameters fully characterize the structure and intensity features of the local image. The specific second-order polynomial parameterization model is as follows: In the formula, Represents the pixel grayscale value. Represents the two-dimensional coordinates of a pixel The matrix formed, symbol Represents transpose, symmetric matrix Characterizing local curvature, gradient vector Reflecting the spatial gradient distribution, the constant term Indicates the baseline grayscale value. Where, symmetric matrix for: gradient vector for: .

4. A method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation, as described in claim 2 or 3, is characterized in that... Step S1.3 specifically involves: introducing a multi-scale pyramid processing framework from coarse to fine; at each scale level, constructing and minimizing a weighted least squares objective function that considers neighborhood consistency; iteratively optimizing the displacement field of that layer; and finally accumulating the results of each layer to obtain a dense two-dimensional displacement field with sub-pixel precision at the original image scale. When solving the displacement field, a filtering window is set in combination with neighborhood information for optimization, and Gaussian weights or uniform weights are used to smooth the displacement field.

5. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 1, characterized in that, Step S2 specifically includes: S2.1: Based on the geometric principles of SAR satellite side-looking imaging, strictly establish the three-dimensional displacement vector of the surface point in the northeast-central geographic coordinate system. The three-dimensional displacement vector is the mathematical projection equation between its two observable components in the single-track SAR image coordinate system. Including vertical direction East-west North-South Direction The two observable components are the radar line-of-sight displacement and the line-of-sight displacement. and satellite flight azimuth displacement The equation uses the satellite incident angle θ and heading angle α as the core geometric parameters; S2.2: By combining the ascending and descending SAR observation data, the projection equations from two independent geometric perspectives are combined to construct a system based on the vertical direction. East-west North-South Direction The overdetermined observation equations are for unknowns, thus providing sufficient mathematical constraints for the effective separation of three-dimensional displacement components; S2.3: Specific imaging geometric parameters based on ascending and descending orbit satellite imagery, i.e., incident angle. , and heading angle , Determine the coefficient matrix of the overdetermined system of equations. Furthermore, the geometric observability of the three-dimensional displacement solution is evaluated by analyzing the condition number or singular value distribution of the coefficient matrix B.

6. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 5, characterized in that, In step S2.2, the overdetermined observation equation set is expressed as: Among them, observation vector , including the observed displacement components of ascending and descending orbits; Vector to be determined ; coefficient matrix It is determined by the incident angle θ and heading angle α of the satellite in ascending or descending orbit.

7. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 5, characterized in that, The coefficient matrix The expression is: 。 8. A method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation, as described in claim 6 or 7, characterized in that... In step S3: If the matrix is ​​nonsingular, the estimated value of the three-dimensional displacement vector D obtained by using the least squares criterion is: ； When the coefficient matrix of the overdetermined observation equation system is ill-conditioned or the observations contain gross errors, a robust estimation strategy based on singular value decomposition is adopted. By truncating small singular values, a generalized inverse matrix is ​​constructed, and a stable least squares solution is obtained under the minimum norm principle.

9. The method for calculating three-dimensional surface displacement in mining areas based on Farneback optical flow method and ascending-descending orbit SAR collaborative observation as described in claim 8, characterized in that, In step S3, the robust estimation strategy based on singular value decomposition is specifically as follows: Perform singular value decomposition on coefficient matrix B , In the formula, U represents an M×N orthogonal matrix, and its column vectors are... eigenvectors; S represents an M×N diagonal matrix whose diagonal elements are The square root of the non-zero eigenvalue; V represents an N×N orthogonal matrix whose column vectors are eigenvectors; The generalized inverse matrix is ​​obtained by reasonably truncating the small singular values ​​in the singular value matrix S. And then through A robust estimate of the three-dimensional displacement vector is obtained.

Images