Phase unwrapping method for differential interferogram of subsidence area based on fused random medium model

Abstract

The application provides a subsidence area differential interference image phase unwrapping method based on a fusion random medium model and belongs to the technical field of image processing. In view of the characteristics of mining subsidence, the method is derived from the mechanism and law of mining subsidence, deduces the surface deformation estimation caused by working face mining based on a random medium model, and then corrects the phase unwrapping result by using the surface deformation estimation, so as to accurately obtain the unwrapping phase of any point on the surface. The method fully considers the mining subsidence mechanism, uses the surface deformation estimation caused by underground coal mining to guide the phase unwrapping result and improve the phase unwrapping precision, and makes up for the defect that the existing unwrapping algorithm is helpless in the fast deformation area of the mine surface. The method can complete the unwrapping of the interferogram phase in the fast deformation area of the mine surface, and compared with the popular minimum cost flow method, has more reasonable unwrapping result.

Description

Technical Field

[0001] This invention relates to a phase unwrapping method for differential interferometric images of sinkhole regions based on a fused random medium model, belonging to the field of image processing technology. Background Technology

[0002] Synthetic Aperture Radar Interferometry (InSAR) technology, with its advantages of high ground resolution, wide coverage, and unique surface observation-based approach, has been widely used in constructing digital elevation models (DEMs) and monitoring surface deformation. However, the phase information in the interferogram is only the phase value modulo 2π. To obtain the absolute phase, phase unwrapping is necessary. The accuracy of phase unwrapping directly affects the accuracy of the constructed DEM and the monitored surface deformation; therefore, it is an extremely important step in SAR interferometric processing.

[0003] Existing phase unwrapping algorithms mainly fall into two categories: path-following methods and minimum norm methods. The former integrates the phase difference between adjacent pixels by selecting an integration path to detect the integer number of cycles in the 2π phase. The most classic method in this category is the Goldstein branching method, and other methods include minimum spanning tree methods and minimum cost flow methods. The latter uses minimizing the difference between the wrapped and unwrapped phase gradients as the criterion, transforming phase unwrapping into a global optimization problem of finding the minimum norm in mathematics. A typical example is the least squares method. Later, many scholars proposed phase unwrapping algorithms with different approaches for specific problems or to achieve specific goals. Xu Junyi et al. proposed a region-based phase unwrapping approach for large-scale, ultra-wideband P-band SAR interferometric phase maps; Huang Haifeng et al. proposed a fast multi-baseline, multi-frequency domain phase unwrapping algorithm for areas with drastic terrain changes; Liu Wanli et al. proposed a Kalman filter-based phase unwrapping method for interferometric images under high-noise conditions in mining areas. Furthermore, Danudirdjo et al. proposed an anisotropic phase unwrapping approach for synthetic aperture radar based on the projection shortening effect of interferometric images.

[0004] InSAR technology is widely used for monitoring mining subsidence due to its advantages of large-scale, high-precision, and continuous dynamic monitoring of surface deformation. However, surface subsidence caused by coal mining differs from other surface deformations, characterized by rapid deformation rates and large subsidence magnitudes. In interferometric phase maps, this manifests as large phase gradients and high fringe density between adjacent pixels. The discontinuity of phase differences between adjacent pixels greatly increases the difficulty of phase unwrapping, making it difficult for existing phase unwrapping algorithms to achieve phase unwrapping in interferometric phase maps of rapidly deforming areas of the mining surface. Summary of the Invention

[0005] To address the problems existing in the prior art, this invention provides a phase unwrapping method for fusion subsidence prediction estimation based on a stochastic medium model. This method starts from the mechanism and laws of mining subsidence, firstly deduces the estimated surface deformation caused by mining at the working face based on the stochastic medium model, and then uses the estimated surface deformation to correct the phase unwrapping result, so as to achieve the purpose of accurately obtaining the unwrapped phase.

[0006] To achieve the above objectives, this invention proposes a phase unwrapping method for differential interferometric images of subsidence areas that integrates a stochastic medium model. Targeting the characteristics of mining subsidence, this method starts from the mechanism and laws of mining subsidence, deduces the estimated surface deformation caused by mining on the working face based on a stochastic medium model, and then uses the estimated surface deformation to correct the phase unwrapping results, thereby accurately obtaining the unwrapped phase at any point on the surface.

[0007] The steps are as follows:

[0008] S1. Perform differential interferometric processing on the SAR image of the target mining area to obtain a differential interferometric fringe pattern containing deformation information on the surface of the target mining area. The differential interferometric processing steps include: registration, resampling, interferometry, and topographic phase removal.

[0009] S2. Using existing conventional phase unwrapping methods: branch cutting method and minimum cost flow method, the differential interferometric fringe pattern of the target mining area surface is unwrapped to obtain the initial unwrapped phase pattern of the target mining area surface.

[0010] S3. Based on the underground mining information and the geological and mining conditions of the target mining area, calculate the estimated value of surface deformation and subsidence caused by underground resource mining in the mining area; using the resolution of the surface SAR image of the target mining area as a benchmark, resample the surface deformation estimates that are inconsistent with the SAR image resolution to make the resolution consistent with the SAR image resolution. The resampling adopts the bilinear interpolation method; the coordinate system of the surface deformation estimate is transferred to the SAR image coordinate system, and a surface deformation map is generated. Then, the gradient between adjacent surface deformation estimate points in the surface deformation map is calculated.

[0011] S4. Based on the basic principle of phase interferometry, the regions in the surface deformation map are classified with λ*u / 2 as the reference. Regions with deformation gradients less than λ*u / 2 are marked as quasi-true deformation regions, and the rest are marked as deformation regions to be revised. Here, λ is the radar wavelength of the SAR image sensor, and u is the image resolution of the SAR image sensor.

[0012] S5. Using SAR radar parameters, the surface deformation information of the two types of regions that have been classified is converted into phase information, and marked as quasi-true phase region and phase region to be revised respectively.

[0013] S6. Define the region with the same coordinate range between the quasi-true phase region obtained in step S5 and the initial unwrapped phase map obtained in step S2 as the common part. Use the least squares method to perform phase matching processing on the common part, correct the quasi-true phase to the initial unwrapped phase, and obtain the correction polynomial f(i,j).

[0014] S7. The phase information of the phase region to be revised is superimposed using the modified polynomial f(i,j) to calculate the number of integer phase cycles between adjacent pixels after the revision.

[0015] S8. By fusing the integer number of phase cycles obtained in step S7 with the initial unwrapped phase obtained in step S2, the final unwrapped phase of the differential interferometric image of the depression area can be obtained.

[0016] Furthermore, in phase unwrapping, to obtain the actual surface deformation, phase ambiguity that is an integer multiple of 2π must be removed to restore the true phase. Therefore, the true phase and the wrapped phase satisfy the following relationship:

[0017]

[0018] That is, from the entanglement phase The optimal integer multiple of the 2π period n(i,j) is sought to achieve phase unwrapping.

[0019] In the formula, φ(i,j) is the wrapped phase value; φ(i,j) is the true phase value; n(i,j) is an integer; the number of pixels in the interferometric phase map is M×N, and the phase value directly observed in the pixel is between (-π,π).

[0020] Furthermore, the specific calculation process of the surface deformation value in step 3 is as follows: using the random medium model, by setting the expected inflection point, the surface deformation value within the inflection point is directly calculated using the random medium model, i.e., formula (2). The surface subsidence value outside the inflection point cannot be directly calculated using the result of the random medium model. It needs to be used after the model is corrected. The corrected result is formula (3). Using the piecewise function method, the surface deformation value after the working face is mined is expected, and then the phase integer number between adjacent pixels is detected based on the expected deformation value:

[0021] Surface subsidence value before the inflection point:

[0022] The stochastic medium model satisfies the following: a) the movement in each direction caused by mining is independent of direction; b) the surface deformation caused by mining of a large underground working face is a linear superposition of the surface deformation caused by mining of multiple small working faces; and the subsidence value W of surface point x caused by unit mining in a two-dimensional plane. e (x) is represented as:

[0023]

[0024] In the formula, r = H0 / tanβ is a constant called the main influence radius, H0 is the average mining depth, tanβ is the tangent of the main influence angle, and x represents the horizontal distance between any point on the ground and the inflection point.

[0025] Surface subsidence values ​​outside the inflection point:

[0026] When using a stochastic medium model to predict mining subsidence, there is a defect of excessively rapid convergence outside the inflection point. The main factor affecting the convergence of the subsidence curve is the constant r. By changing the value of r, the subsidence curve can be controlled. The ratio of the value of r outside the inflection point to the value of r inside the inflection point is approximately a constant k. Therefore, outside the inflection point, the subsidence value of the surface point caused by unit mining in the two-dimensional plane is expressed as:

[0027]

[0028] In the formula, k is a constant, defined as the convergence coefficient of the subsidence curve;

[0029] Expression of ground subsidence at any point in a three-dimensional plane:

[0030] The subsidence of surface point A(x,y) is the result of mining the coal seam in both the strike and dip directions. Under three-dimensional conditions, the subsidence value W(x,y) of surface point A(x,y) caused by mining in a certain mining area is... A Represented as:

[0031]

[0032] In the formula, f(x,y) is the spatial probability density function, and F represents the region, that is, the region where x∈[0,l] and y∈[0,L].

[0033] Let the length of the rectangular mining area along the strike be l, and the width along the dip be L. Let the mining dimensions be x∈[0,l] and y∈[0,L]. Considering the independence of the probabilities in the x and y directions, when the lower left corner of the mining area is taken as the origin, the stochastic model can be expressed as a piecewise function. The surface deformation value can be further expressed as:

[0034]

[0035] In the formula, r is the radius of influence of the coal seam strike. 1(2) These represent the main radii of influence of the coal seam in the uphill and downhill directions, respectively; W max =m·q·cosα is the maximum subsidence value of the ground surface, m is the coal seam mining thickness, q is the subsidence coefficient, and α is the coal seam dip angle. Then, the subsidence gradient function of the ground point A(x,y) in the azimuth direction β is the first derivative of the subsidence. The gradient in a certain direction is denoted as W'(x,y). β .

[0036] Furthermore, the principle of region classification is to compare the gradient W'(x,y) in a certain direction. β The relationship between the size of λ*u / 2 and the actual deformation zone is determined as follows: if the former is smaller, it is designated as the quasi-true deformation zone; otherwise, it is designated as the deformation zone to be revised.

[0037] Using formula (7), the surface deformation information of the two types of regions that have been classified is converted into phase information, and marked as quasi-true phase region and phase region to be revised respectively;

[0038] Using the formula f(i,j)=a*i 2 +b*j 2 The formula +c*ij+d*i+e*j+f calculates the correction polynomial, where a, b, c, d, e, and f are the parameters to be determined, and (i,j) are the plane coordinates. To determine the expression for f(i,j), we need to find the six parameters a, b, c, d, e, and f. The method for determining these six parameters is to select the same points in the common part and solve them using the least squares method.

[0039] Calculate the corrected phase φ′ using the formula. defo (i,j)=φ defo (i,j)+f(i,j), and the number of integer cycles of the phase.

[0040] Furthermore, based on the unwrapped phase expression after phase correction:

[0041] According to the basic principles of D-InSAR technology, the vertical surface subsidence W is related to the deformation phase along the radar line of sight. The relationship is as follows:

[0042]

[0043] In the formula, λ is the radar satellite wavelength; θ is the line-of-sight angle; The deformation phase along the radar line of sight is the remaining part after removing the reference ellipsoid phase, terrain phase, atmospheric delay phase, noise phase, etc. from the interferometric phase diagram. Therefore, the formula for converting surface deformation into phase in step S5 is:

[0044]

[0045] Therefore, according to step S7, the number of integer cycles n(x,y) of the predicted phase of the table sinking is:

[0046] n(x,y)=Int[φ′ defo / (2π)] (8)

[0047] In the formula, Int[·] represents the floor function, and φ′ defo Indicates the corrected deformation phase;

[0048] Substituting equation (8) into (1) yields the unwrapping phase at any point on the ground in step S8:

[0049]

[0050] The initial phase unwrapping result is corrected using the phase integer number n(x, y).

[0051] Beneficial effects:

[0052] This invention addresses the problem that existing conventional phase unwrapping methods are unable to accurately achieve differential interferometric phase unwrapping in mining areas. Considering the characteristics of mining subsidence, the proposed method starts from the mechanism and laws of mining subsidence. First, it uses a stochastic medium model to estimate the surface deformation caused by mining at the working face. Then, it uses the estimated surface deformation to correct the phase unwrapping result, so as to achieve the goal of accurately obtaining the unwrapped phase.

[0053] This invention fully considers the mechanism of mining subsidence, using the estimated surface deformation caused by underground coal mining to guide the phase unwrapping results and improve the accuracy of phase unwrapping. The proposed method overcomes the shortcomings of existing unwrapping algorithms in areas of rapid surface deformation in mining areas. The proposed phase unwrapping algorithm based on subsidence prediction can complete the unwrapping of the interferogram phase in areas of rapid surface deformation in mining areas, and has more reasonable unwrapping results compared with the popular minimum cost flow method. Attached Figure Description

[0054] Figure 1 This is a schematic diagram of phase unwrapping according to the present invention. The upper diagram shows the absolute phase, and the lower diagram shows the wrapped phase.

[0055] Figure 2 This is a schematic diagram of the phase unwrapping method for differential interferometric image maps of the depression region based on the random medium model of the present invention:

[0056] Figure 3 A two-dimensional coordinate system diagram;

[0057] Figure 4 A three-dimensional coordinate system diagram;

[0058] Figure 5 This is a schematic diagram of Example 1. In the figure, (a): the simulated mining face and the surface subsidence contour lines predicted by the random medium model; (b): the true phase of the simulated surface deformation; (c): the entanglement phase of the simulated surface deformation; (d): the integer number of phases detected by this method; (e): the unentangled phase after unentanglement by the method in this paper; (f): a histogram of the difference between the true phase and the unentangled phase.

[0059] Figure 6This is a schematic diagram of Example 2. In the figure, (a): contour lines of surface subsidence predicted by the mining face and random medium model; (b): surface deformation caused by underground mining obtained by a ground 3D laser scanner; (c): the true phase of surface deformation caused by underground mining; (d): the entanglement phase of surface deformation caused by underground mining; (e): the integer number of phases detected by this method; (f): the unentangled phase after untangling by this method.

[0060] Figure 7 This is a schematic diagram of Example 3. In the figure, (a): the surface subsidence contour lines predicted by the mining face and the random medium model; (b): the entanglement phase of surface deformation caused by underground mining; (cb): the integer number of phase cycles detected by the minimum cost flow method phase unwrapping method and the method proposed in this paper; (d): the unwrapped phase after unwrapping by this method.

[0061] Figure 8 This is a statistical histogram of the difference between the unwound phase and the original wrapped phase after rewinding. Detailed Implementation

[0062] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0063] Untangling methods

[0064] In an M×N interferometric phase diagram, the phase values ​​directly observed in a pixel are between (-π, π), and are called the entangled phase. To obtain the actual surface deformation, phase ambiguity that is an integer multiple of 2π must be eliminated to restore the true phase, ensuring that the true phase and the entangled phase satisfy the following relationship:

[0065]

[0066] In the above formula, Let φ(i,j) be the entangled phase value; φ(i,j) be the true phase value; and n(i,j) be an integer. Phase unwrapping is the process of unwrapping the entangled phase. The process of seeking the optimal integer multiple of 2π period n(i,j) involves the following one-dimensional entanglement and unentanglement phase relationship: Figure 1 As shown.

[0067] Depend on Figure 1Theoretically, if the absolute value of the phase difference between adjacent pixels is less than π, phase unwrapping can be achieved by integrating the phase difference. However, in real SAR image interferograms, factors such as radar shadows, phase noise, and rapid surface deformation cause discontinuities in the phase data, making the integration of the phase difference extremely difficult. In particular, when using D-InSAR technology to monitor surface deformation in mining areas, the large-scale subsidence that occurs shortly after underground coal mining leads to significant discontinuities between adjacent pixels in the interferometric phase map. This results in inaccurate integration of the phase difference during phase unwrapping, affecting the unwrapping effect.

[0068] Untangling Strategy

[0069] Compared to the unpredictability of surface deformation caused by earthquakes, landslides, and glacial movement, surface deformation caused by underground mining follows certain patterns, and the deformation value and spatial distribution at any point on the surface can be predicted. Therefore, we can use the difference between the predicted results of two adjacent pixels to detect the integer number of phase jumps, thereby achieving phase unwrapping.

[0070] Figure 2 A flowchart based on this approach is provided, with the specific process explained below:

[0071] 1) Perform differential interferometric processing on SAR images (including registration, resampling, interferometry, and terrain phase removal) to obtain a differential interferometric fringe pattern containing deformation information of the target area.

[0072] 2) Use existing conventional phase unwrapping methods (such as branch cutting method, minimum cost flow method, etc.) to unwrap the interferometric phase map to obtain the initial unwrapped phase map of the target area.

[0073] 3) Based on the underground mining information and the geological and mining conditions of the target area, calculate the estimated value of surface deformation caused by underground resource mining; resample the calculated surface deformation values ​​using SAR image resolution as a benchmark; unify them to the radar coordinate system and calculate the deformation gradient between adjacent points.

[0074] 4) Based on the fundamental principle of phase interferometry, and using λ*u / 2 as the reference (λ and u are the radar wavelength and image resolution of the SAR image sensor, respectively), the deformation map obtained in step 3 is classified. The classification principle is: regions with deformation gradients less than λ*u / 2 are designated as quasi-true deformation regions, and the rest are deformation regions to be revised.

[0075] 5) Using radar sensor parameters, based on formula (7), the surface deformation of the two types of areas classified in step 4 is converted into phase information, and recorded as the quasi-true phase area and the phase area to be revised respectively.

[0076] 6) The common part of the quasi-true phase region obtained in step 5 and the initial unwrapped phase map obtained in step 2 is subjected to phase matching processing using the least squares method. The quasi-true phase is corrected to the initial unwrapped phase, and the corrected polynomial f(i,j) is obtained.

[0077] 7) The phase information of the phase region to be revised is corrected using the modified polynomial f(i,j), and the number of integer cycles of the corrected phase is calculated using formula (8).

[0078] 8) The final unwound phase is obtained by fusing the integer number of phases obtained in step 7 with the initial unwound phase obtained in step 2.

[0079] To implement the above method, two problems need to be solved: the first is the prediction of surface deformation values, and the second is the expression of the unwrapping phase.

[0080] Surface deformation caused by mining is expected

[0081] Mining subsidence prediction is a core issue in mining subsidence science. Among numerous subsidence prediction methods, the stochastic medium model proposed by Polish scholar J. Litwinisyn is one of the most mature methods. This method has been widely used in China, and extensive production practice has confirmed that the accuracy of the predicted deformation values ​​within the inflection point can reach over 90% (the location of the inflection point is as follows). Figure 3 (As shown); for the portion outside the inflection point, after model correction, the predicted value of the surface deformation result can also reach 90%. Therefore, it is theoretically reliable to use a stochastic medium model and a piecewise function to predict the surface deformation value after mining, and then to detect the phase integer number between adjacent pixels based on the predicted deformation value.

[0082] 1) Expression of surface subsidence value before the inflection point

[0083] The stochastic medium model has two fundamental assumptions: a) the movement in each direction caused by mining is independent of direction; b) the surface deformation caused by mining of a large underground working face is a linear superposition of the surface deformation caused by mining of multiple small working faces. Under these two assumptions, the subsidence W of surface point x caused by unit mining in a two-dimensional plane is... e (x) can be represented as:

[0084]

[0085] In the formula, r = H0 / tanβ is a constant, called the primary influence radius. Here, H0 is the average mining depth; tanβ is the tangent of the primary influence angle.

[0086] 2) Expression of surface subsidence values ​​outside the inflection point

[0087] When using a stochastic medium model to predict mining subsidence, there is a defect of excessively rapid convergence outside the inflection point. Equation (2) shows that the main factor affecting the convergence of the subsidence curve is r; as long as r changes, the boundary convergence will change. Therefore, by changing the value of r, the subsidence curve can be controlled. Studies have found that the ratio of the value of r outside the inflection point to the value of r inside the inflection point is approximately a constant k. Therefore, outside the inflection point, the subsidence value of the surface point caused by unit mining in the two-dimensional plane can be expressed as:

[0088]

[0089] In the formula, k is a constant, which can be defined as the convergence coefficient of the subsidence curve.

[0090] 3) Expression of ground subsidence value at any point in a three-dimensional plane

[0091] Surface movement is a three-dimensional problem, meaning the subsidence of surface point A(x,y) is a result of mining the coal seam along both the strike direction (x-axis) and the dip direction (y-axis). Under three-dimensional conditions, what is the subsidence value W(x,y) of surface point A(x,y) caused by mining in a certain mining area? A It can be represented as:

[0092]

[0093] In the formula, f(x,y) is the spatial probability density function.

[0094] like Figure 4 As shown, assuming a rectangular mining area has a mining length of l along the strike and a mining width of L along the dip, with mining dimensions x∈[0,l] and y∈[0,L] respectively, considering the independence of probabilities in the x and y directions, when the lower left corner of the mining area is taken as the origin, equation (4) can be expressed as:

[0095]

[0096] In the formula, r is the radius of influence of the coal seam strike. 1(2) These represent the main radii of influence of the coal seam in the uphill and downhill directions, respectively; W max =m·q·cosα is the maximum subsidence value of the surface, m is the coal seam mining thickness, q is the subsidence coefficient, and α is the coal seam dip angle.

[0097] The settlement gradient function of ground point A(x,y) in the azimuth direction β is the first derivative of the settlement, denoted as W'(x,y). β .

[0098] Untangled phase representation based on predicted results

[0099] According to the basic principles of D-InSAR technology, the vertical surface subsidence W is related to the deformation phase along the radar line of sight. The relationship between them can be represented as:

[0100]

[0101] In the formula, λ is the radar satellite wavelength; θ is the line-of-sight angle; The deformed phase along the radar line of sight is the remaining portion after removing the reference ellipsoid phase, terrain phase, atmospheric delay phase, noise phase, etc., from the interferometric phase diagram. Therefore:

[0102]

[0103] Therefore, the number of integer phase cycles n(x,y) detected based on the predicted surface subsidence can be expressed as:

[0104]

[0105] In the formula, Int[·] represents the floor operation, and the meanings of the other parameters are the same as those described above.

[0106] Substituting equation (8) into equation (1) yields the unwrapped phase at any point on the ground, as shown in equation (9) below:

[0107]

[0108] In the formula, W(x,y) is a piecewise function, and its specific expression is (5).

[0109] Example 1

[0110] A mining face was simulated using the geological and mining conditions in Table 1. The mining depth and size of the face, as well as the predicted surface subsidence contour lines after mining using a stochastic medium model, are shown below. Figure 5 As shown in (a); Figure 5 (b) and Figure 5 (c) These are the true and entangled phases of surface deformation monitored by the RADARSAT-2 satellite using differential interferometry, simulated with the satellite parameters in Table 2. (i.e., the entangled phases that need to be unwrapped in differential interferometry processing.)

[0111] Table 1. Geological and mining conditions for simulated mining.

[0112]

[0113] Table 2 Basic parameters of simulated SAR satellites

[0114]

[0115] Entangled phase estimation

[0116] In reality, the parameters used in predicting mining subsidence cannot fully and accurately represent geological mining conditions. This is why the prediction of surface deformation caused by underground coal seam mining using stochastic medium models is not entirely accurate. Therefore, when applying the derived formulas for phase integer cycle detection, another set of parameters from Table 3 is used to reflect the reality that does not conform to actual mining conditions.

[0117] Table 3 Parameters used for phase unwrapping

[0118]

[0119] Using equations (8), (9), and the parameters in Table 3, the integer number of phase cycles was detected and the final unwrapping phase was calculated, and the results are as follows: Figure 5 As shown in (d) and (e) in the figure; (a): contour lines of surface subsidence predicted by the simulated mining face and random medium model; (b): true phase of simulated surface deformation; (c): entanglement phase of simulated surface deformation; (d): integer number of phases detected by this method; (e): unentangled phase after unentanglement by the method in this paper; (f): histogram of the difference between the true phase and the unentangled phase;

[0120] Results Evaluation

[0121] Will Figure 5 The untangling phase in (e) and Figure 5 (b) Perform a difference operation on the true phase. Figure 5 (f) is the statistical histogram of the phase difference. From Figure 5 As shown in (f), the number of pixels with negative phase differences is much greater than the number of pixels with positive phase differences. This is because the sinking coefficient q used when detecting the integer number of phase cycles is too small, which makes the detected integer number of phase cycles smaller than the true integer number, thus increasing the probability that the unwrapped phase is smaller than the true phase. The average difference and standard deviation of the phase difference are -0.0058 rad and 1.7405 rad, respectively. Both are less than π, indicating that the phase unwrapping method for differential interferometric images of mining areas proposed in this paper is feasible under noise-free conditions.

[0122] Example 2: Wannian Mine 13266 Working Face

[0123] Wannian Mine is located in Fengfeng Mining District, Handan City, Hebei Province. The mining area of ​​working face 13266 from December 4, 2009 to March 11, 2010, and the predicted surface subsidence contour lines after mining using a stochastic medium model are as follows: Figure 6 (a) is shown below (the expected parameters are shown in Table 4 below); Figure 6(b) The spatial distribution of surface deformation obtained by a 3D laser scanner during the mining operation; the true phase and entanglement phase of the surface deformation simulated using the satellite parameters in Table 2 (i.e., the entanglement phase that needs to be unwrapped in differential interferometry) are respectively as follows: Figure 6 (c) and Figure 6 As shown in (d).

[0124] Table 4. Parameters used for phase unwrapping

[0125]

[0126] Using equations (8), (9), and the parameters in Table 4, the integer number of phase cycles was detected and the final unwrapping phase was calculated, and the results are as follows: Figure 6 (e) and Figure 6 As shown in (f).

[0127] Subtracting the unwrapped phase from the true phase yields an average phase difference of -1.56 rad, which is less than π. If the phase discrepancy after unwrapping is considered as noise, the calculated signal-to-noise ratio is 58.56. Both indicators demonstrate the high effectiveness of the proposed method in unwrapping the entangled phase caused by surface deformation resulting from underground mining.

[0128] Example 3

[0129] The study area is the 15223 working face of Jiulong Mine, which also belongs to Fengfeng Mining Area, Handan City, Hebei Province. The mining scope of the working face during the period from January 3, 2014 to March 6, 2014, and the predicted surface subsidence contour lines after mining using the stochastic medium model are as follows: Figure 7 As shown in (a), the expected parameters are shown in Table 5:

[0130] Table 5. Parameters used for phase unwrapping

[0131]

[0132]

[0133] The image data used for differential interferometry in this experiment is from the RadarSat-2 radar satellite. The acquisition date and related parameters of the images are shown in Table 6 below.

[0134] Table 6 RadarSat-2 Image Parameter Information

[0135]

[0136] Wrapping phase estimation:

[0137] During the experiment, GAMMA software was used, and the conventional "two-track method" was employed to perform differential interferometric processing on SAR images. The resulting differential interferometric phase is as follows: Figure 7 As shown in (b).

[0138] To obtain the true phase caused by surface deformation, we used two methods to estimate the entanglement phase: the minimum cost flow method and the phase unwrapping method based on subsidence prediction proposed in this paper. The results are as follows: Figure 7 (c) and (d) in the middle.

[0139] Result evaluation:

[0140] from Figure 7 (c) It can be seen that due to the high fringe density and phase noise in the interferometric phase map, the minimum cost flow method cannot achieve phase unwrapping of all pixels in the study area; while the phase unwrapping method based on subsidence prediction proposed in this paper not only fully considers the mechanism of mining subsidence, but also overcomes the shortcomings of conventional phase unwrapping methods in areas with rapid surface deformation in mining areas.

[0141] Will Figure 7 (d) is rewound and connected with the untangled phase. Figure 7 The interference in (b) is subjected to difference calculation, and the statistical histogram of the phase difference is as follows: Figure 8 As shown. From Figure 8 It can be seen that the phase difference between the unwound phase and the original wrapped phase after rewinding conforms to a normal distribution. The average difference and standard deviation of the phase difference are -0.0020 rad and 0.9011 rad, respectively, both of which are less than π. This further indicates that the method can achieve phase unwinding well while maintaining a high degree of consistency with the original wrapped phase.

[0142] To address the difficulty in achieving phase unwrapping in differential interferometric phase maps for mining areas, this paper proposes a phase unwrapping method based on subsidence prediction, starting from the mechanism of mining subsidence and utilizing the fundamental theory of stochastic medium models. Based on the predicted values ​​of surface deformation caused by underground mining, the unwrapped phase of any point on the surface in the interferometric phase map is derived.

[0143] This method was validated for its feasibility, effectiveness, and applicability using noise-free simulated data, noisy measured 3D laser scanning data, and real InSAR imagery. Theoretical analysis and practical verification show that this method can not only achieve phase unwrapping of differential interferometric images in rapidly deformable areas of mining surfaces, but also has excellent unwrapping performance, thereby further improving the application of D-InSAR technology in mining areas.

Claims

1. A phase unwrapping method for differential interferometric images of sinkhole regions based on a fused random medium model, characterized in that: In view of the characteristics of mining subsidence, starting from the mechanism and law of mining subsidence, the estimated surface deformation caused by mining is derived based on the stochastic medium model. Then, the estimated surface deformation is used to correct the phase unwrapping result, so as to accurately obtain the unwrapped phase at any point on the surface. The steps are as follows: S1. Perform differential interferometric processing on the SAR image of the target mining area to obtain a differential interferometric fringe pattern containing deformation information on the surface of the target mining area. The differential interferometric processing steps include: registration, resampling, interferometry, and topographic phase removal. S2. Using existing conventional phase unwrapping methods: branch cutting method and minimum cost flow method, the differential interferometric fringe pattern of the target mining area surface is unwrapped to obtain the initial unwrapped phase pattern of the target mining area surface. S3. Based on the underground mining information and the geological and mining conditions of the target mining area, calculate the estimated value of surface deformation and subsidence caused by underground resource mining in the mining area; using the resolution of the surface SAR image of the target mining area as a benchmark, resample the surface deformation estimates that are inconsistent with the SAR image resolution to make the resolution consistent with the SAR image resolution. The resampling adopts the bilinear interpolation method; the coordinate system of the surface deformation estimate is transferred to the SAR image coordinate system, and a surface deformation map is generated. Then, the gradient between adjacent surface deformation estimate points in the surface deformation map is calculated. S4. Based on the basic principle of phase interferometry, the regions in the surface deformation map are classified with λ*u / 2 as the reference. Regions with deformation gradients less than λ*u / 2 are marked as quasi-true deformation regions, and the rest are marked as deformation regions to be revised. Here, λ is the radar wavelength of the SAR image sensor, and u is the image resolution of the SAR image sensor. S5. Using SAR radar parameters, the surface deformation information of the two types of regions that have been classified is converted into phase information, and marked as quasi-true phase region and phase region to be revised respectively. S6. Define the region with the same coordinate range between the quasi-true phase region obtained in step S5 and the initial unwrapped phase map obtained in step S2 as the common part. Use the least squares method to perform phase matching processing on the common part, correct the quasi-true phase to the initial unwrapped phase, and obtain the correction polynomial f(i,j). S7. The phase information of the phase region to be revised is superimposed using the modified polynomial f(i,j) to calculate the number of integer phase cycles between adjacent pixels after the revision. S8. By fusing the integer number of phase cycles obtained in step S7 with the initial unwrapped phase obtained in step S2, the final unwrapped phase of the differential interferometric image of the depression area can be obtained.

2. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 1, characterized in that, In phase unwrapping, to obtain the actual surface deformation, it is necessary to remove... For phase ambiguity that is an integer multiple of the actual phase, the true phase is restored. Therefore, the following relationship exists between the true phase and the entangled phase: (1)， That is, from the entangled phase Seeking the best Period integer multiple This achieves phase untangling; In the formula, This represents the entanglement phase value; This is the true phase value; The values ​​are integers; the interferometric phase map has M×N pixels, and the phase values ​​directly observed in each pixel are in... between.

3. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 2, characterized in that, The specific calculation process of the surface deformation value in step 3 is as follows: Using the random medium model, by setting the expected inflection point, the surface deformation value within the inflection point is directly calculated using the random medium model, i.e., formula (2). The surface subsidence value outside the inflection point cannot be directly calculated using the result of the random medium model. It needs to be used after the model is corrected. The corrected result is formula (3). Using the piecewise function method, the surface deformation value after the working face is mined is expected, and then the phase integer number between adjacent pixels is detected based on the expected deformation value: Surface subsidence value before the inflection point: The stochastic medium model satisfies the following: a) the movement in each direction caused by mining is independent of direction; b) the surface deformation caused by mining of a large underground working face is a linear superposition of the surface deformation caused by mining of multiple small working faces; and the surface deformation caused by mining of unit units in a two-dimensional plane is... sinking value Represented as: (2)， In the formula, The constant is called the principal influence radius. It is the average mining depth. It is the main influencing angle tangent, where x represents the horizontal distance between any point on the ground and the inflection point; Surface subsidence values ​​outside the inflection point: When using a stochastic medium model to predict mining subsidence, there is a defect of excessively rapid convergence outside the inflection point. The main factor affecting the convergence of the subsidence curve is the constant r. By changing the value of r, the subsidence curve can be controlled. The ratio of the value of r outside the inflection point to the value of r inside the inflection point is approximately a constant k. Therefore, outside the inflection point, the subsidence value of the surface point caused by unit mining in the two-dimensional plane is expressed as: (3)， In the formula, k is a constant, defined as the convergence coefficient of the subsidence curve; Expression of ground subsidence at any point in a three-dimensional plane: surface point The subsidence is a result of mining the coal seam in both the strike and dip directions. Under three-dimensional conditions, the subsidence of a certain mining area causes surface points to... sinking value Represented as: (4)， In the formula, Let F be the spatial probability density function, and let F represent the region. The area; Let the length of the rectangular mining area along the strike be... The width of the mining along the dip is The mining dimensions are respectively Considering The independence of directional probabilities, when the lower left corner of the mining area is taken as the origin of the coordinate system, allows for a piecewise function representation of the stochastic model. The surface deformation value can then be further expressed as: (5)， In the formula, The radius of influence is mainly determined by the strike of the coal seam. These represent the main influence radii of the coal seam in the uphill and downhill directions, respectively. This represents the maximum surface subsidence value. For coal seam mining thickness, This is the subsidence coefficient. Let β be the dip angle of the coal seam. Then, the settlement gradient function of the ground point A(x, y) in the direction of azimuth β is the first derivative of the settlement. The gradient in a certain direction is denoted as W'(x, y). β .

4. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 3, characterized in that, The principle of region classification is to compare the gradient W'(x, y) in a certain direction. β The relationship between the size of λ*u / 2 and the size of the former is such that if the former is smaller, it is defined as the quasi-real deformation region; Otherwise, it will be designated as a deformation area to be revised; Using formula Calculate the corrected polynomial, where a, b, c, d, e, and f are the parameters to be determined, and (i, j) are the plane coordinates; determine The expression is to find the six parameters a, b, c, d, e, and f. The method to determine these six parameters is to select the same points in the common part and solve them using the least squares method. Calculate the corrected phase using the formula and the number of integer cycles of the phase; This represents the deformation phase along the radar line of sight.

5. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 4, characterized in that, Based on the phase-corrected untangled phase representation: Based on the basic principles of D-InSAR technology, vertical subsidence of the earth's surface Deformation phase along the radar line of sight The relationship is as follows: (6) In the formula, For radar satellite wavelength; The viewing angle; The deformation phase along the radar line of sight is the remaining part after removing the reference ellipsoid phase, terrain phase, atmospheric delay phase, noise phase, etc. from the interferometric phase diagram. Therefore, the formula for converting surface deformation into phase in step S5 is: (7)。 6. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 5, characterized in that, According to step S7, the number of integer cycles of the phase of the predicted subsidence value is determined. : (8) In the formula This represents the integer operation. This indicates the corrected deformation phase.

7. The phase unwrapping method for differential interferometric images of the sinkhole region based on the fused random medium model according to claim 6, characterized in that, Substituting equation (8) into (1) yields the unwrapping phase at any point on the ground in step S8: (9) The initial phase unwrapping result is corrected using the phase integer number n(x, y).

Images