DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification

Abstract

The invention discloses a DInSAR and leveling data fusion earth surface movement parameter inversion method based on earth surface subsidence velocity identification. The method comprises the steps that the subsidence velocity of any earth surface point at any moment in the mining process of a working face is calculated; determining a maximum vertical deformation speed theoretical value between two adjacent scene DInSAR images; the sinking speed of a certain point on the earth surface in the whole mining process of the working face is compared with the maximum vertical deformation speed theoretical value which can be monitored by the DInSAR, and whether the point is subjected to decoherence or not is judged; combining the DInSAR monitoring data of the non-decoherence area and the ground leveling data of the decoherence area to invert and solve the predicted parameters of the probability integral method; and according to a mine mining plan, surface subsidence deformation prediction is carried out by adopting a probability integral method. When the method is used for monitoring the surface subsidence deformation, the ground leveling workload can be reduced, effective data monitored by the DInSAR can be fully utilized, the loss of monitoring precision is not caused, and the accurate prediction of the surface subsidence is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of ground subsidence monitoring, and more particularly to a method for inverting ground movement parameters based on DInSAR and leveling data fusion based on ground subsidence velocity identification. Background Art

[0002] In recent years, with the increasing amount of coal resources mined, the large gradient deformation of the surface after coal mining has endangered people's daily life. The surface monitoring of the whole process of coal mining has become the basis for the reasonable control of surface subsidence. At present, the means of monitoring surface subsidence in various mining areas include leveling, RTK measurement, InSAR monitoring, etc. Although the data accuracy of the leveling and RTK monitoring of surface deformation is relatively reliable, the point observation, small amount of observation data, time-consuming and laborious observation work, and high observation cost are time-consuming and labor-intensive, and the amount of observation data is small. InSAR monitoring technology, which has emerged in recent years, has many advantages such as high precision, wide range, continuous monitoring, all-day and all-weather operation, convenient implementation, low cost, and can realize area monitoring. Although it makes up for the shortcomings of traditional monitoring technology to a certain extent, when the deformation speed of the mining area is fast, InSAR monitoring means will appear incoherence in the large gradient deformation area, resulting in low accuracy of monitoring results. Obviously, for the current mining surface subsidence monitoring technology, the above methods all have limitations and cannot obtain the full process surface subsidence information caused by working face mining in the observation area.

[0003] Due to the limitations of existing monitoring technologies, the monitoring process cannot accurately obtain surface subsidence information throughout the entire mining process of the working face. Therefore, it is urgent to propose a method for monitoring mining subsidence deformation to obtain more complete and high-precision measurement results. Summary of the invention

[0004] To this end, the technical problem to be solved by the present invention is to provide a surface movement parameter inversion method based on the fusion of DInSAR and leveling data based on surface subsidence velocity identification, which can avoid the time-consuming and labor-intensive leveling and the high cost, and avoid the problem of decoherence of DInSAR monitoring data caused by rapid deformation of the mining area.

[0005] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0006] The inversion method of surface movement parameters based on DInSAR and leveling data fusion based on surface subsidence velocity identification includes the following steps:

[0007] (1) Based on the sinking velocity of any point on the surface at any time in the subsidence area and the leveling data of the observation station in the target monitoring area, the sinking velocity V of any point on the surface at any time during the mining process of the working face is calculated. (x,y,t) ;

[0008] (2) Collect SAR data within the target monitoring area, perform DInSAR data processing, obtain the vertical deformation of the surface points, and determine the maximum vertical deformation velocity theoretical value V between two adjacent images. SAR ;

[0009] (3) The sinking velocity V of a point on the ground during the whole mining process is (x,y,t) The theoretical value V of the maximum vertical deformation velocity that can be monitored between two adjacent images of the DInSAR monitoring data at this point SAR By comparison, it is determined whether the DInSAR monitoring data of the point is incoherent during the whole process of mining at the working face, and then the incoherent area and non-incoherent area of ​​the DInSAR monitoring data of the mining subsidence area of ​​the working face are determined;

[0010] (4) Combine the DInSAR monitoring data of the non-decoherent area and the ground leveling data of the decoherent area to invert and obtain the estimated parameters using the probability integral method;

[0011] (5) According to the mine mining plan, the probability integral method is used to predict surface subsidence and deformation.

[0012] The above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification includes the following steps in step (1):

[0013] (1-1) Using the formula for calculating the subsidence velocity of a surface point in the subsidence area at any time and the formula for the spatial distribution function of the subsidence velocity of any surface point in the subsidence area, the formula for the subsidence velocity of any surface point in the subsidence area at any time is obtained;

[0014] (1-2) Collecting leveling data within the target monitoring area and obtaining the sinking velocity of a certain leveling point;

[0015] (1-3) Using the sinking velocity value of a certain leveling measurement point, the parameters in the sinking velocity formula at any point on the surface of the subsidence area at any time are inverted to obtain the sinking velocity at any point on the surface of the subsidence area at any time.

[0016] The above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification is

[0017] In step (1-1), the formula for calculating the sinking velocity of a surface point in the subsidence area at any time is:

[0018] v(t)=W m CK -ct (1-e -ct ) k-1 (1);

[0019] In formula (1), W m is the maximum subsidence value after the subsidence of a certain point on the surface stabilizes, t represents any time after coal seam mining, c is the lithology parameter, and k is the parameter to be fitted;

[0020] The formula for the spatial distribution function of the sinking velocity at any point on the surface of the subsidence area is:

[0021]

[0022] In formula (2), v max is the maximum sinking velocity of the surface point in the subsidence area during the subsidence process; x is the normal distance from the boundary of the working face in the strike direction, y is the normal distance from the boundary of the working face in the dip direction, and l is the maximum sinking velocity lag distance;

[0023] a is the morphological parameter of the main section of the strike, v x is the sinking velocity at any point on the strike main section;

[0024] b is the morphological parameter of the main section, according to the empirical formula Sure, is the average mining depth;

[0025] Substituting formula (1) into formula (2), we can obtain the formula for the sinking velocity at any point on the surface of the subsidence area at any time:

[0026]

[0027] The above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification is

[0028] In step (1-2), the leveling data in the target monitoring area is collected to calculate the sinking speed of a certain leveling point n at time t:

[0029]

[0030] In formula (4): W n (t) and W n (t-1) represents the sinking value of point n at the tth moment and the t-1th moment respectively; Δt is the number of days between two observations; W n (t) = H t -H0, where H t and H0 are the elevations of point n at time t and the first observation respectively.

[0031] In the above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification, in step (1-3),

[0032] According to the measured data of ground leveling, the sinking velocity value of a certain leveling point during the mining process is used to invert the parameters c and k by formula (1); then, according to the position of the surface point (x, y) that needs to be calculated, the sinking velocity of any surface point at any time in the subsidence area during the whole mining process of the working face is calculated by formula (3).

[0033] In the above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification, in step (2), SAR data in the target monitoring area is collected, and the DInSAR data processing method is as follows:

[0034] The dual-track method is used for data processing. The interference phase diagram is obtained by interfering two SAR images: the orbit error phase, atmospheric delay phase and random noise phase in the interference phase are ignored; the flat ground phase is eliminated by removing the flat ground effect and calculating the interference baseline; the terrain phase is eliminated by simulating with the digital elevation model DEM, and the surface deformation phase is finally retained;

[0035] According to the surface deformation phase, the surface vertical deformation R0 is inverted using equation (5):

[0036]

[0037] In formula (5), Δφ is the phase difference of the image pair, λ is the wavelength of the SAR radar signal, and θ is the angle between the radar line of sight and the vertical direction;

[0038] The maximum vertical deformation velocity theoretical value of the satellite used is determined according to the maximum vertical deformation velocity theoretical value that can be monitored between two adjacent images (6):

[0039]

[0040] In formula (6), Δt represents the time interval between two adjacent image data, λ / 2 represents half wavelength, and θ is the incident angle of the SAR radar.

[0041] The above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification, in step (3), during the whole process of mining, the subsidence velocity of any point on the surface at any time is obtained by formula (3), and the maximum vertical deformation velocity theoretical value V that can be monitored by DInSAR at this point is calculated. SAR Compare and determine whether decoherence occurs.

[0042] The above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification is compared as follows:

[0043] If the sinking velocity V of a point on the surface (x,y,t)During the whole process of mining subsidence, the maximum vertical deformation velocity V that can be monitored by DInSAR is exceeded once or multiple times. SAR When , the DInSAR monitoring data of this point loses coherence;

[0044] If the sinking velocity V of a point on the surface (x,y,t) During the whole process of mining, the maximum vertical deformation velocity V that can be monitored by DInSAR was not exceeded. SAR When , the InSAR data at this point does not lose coherence;

[0045] According to the above method, the surface subsidence basin is divided into two parts based on the DInSAR monitoring data, including the peripheral non-decoherent area and the internal decoherent area.

[0046] In the above-mentioned method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification, in step (4), the vertical deformation variables monitored by DInSAR in the non-decoherent area are selected as the measured data, and the observed surface subsidence data of the leveling point in the decoherent area are selected as the measured data, and the probability integral method prediction parameters are inverted respectively;

[0047] The formula for predicting sinking at any point in the probability integral model is:

[0048]

[0049] In formula (7), (x, y) is the coordinate of any point on the surface, W (x,y) is the sinking value at any point;

[0050] W0 is the maximum surface subsidence value, W0 = m × q × cos α, m is the coal thickness, q is the subsidence coefficient, and α is the coal seam inclination;

[0051] erf is the probability integral function, which is in the form of u is the integral parameter, r is the main influence radius, is the average mining depth; tanβ is the tangent value of the main influencing angle; r1=H1 / tanβ; r2=H2 / tanβ; r1 is the main influencing radius in the uphill direction; r2 is the main influencing radius in the downhill direction; H1 represents the mining depth of the working face in the uphill direction, and H2 represents the mining depth of the working face in the downhill direction;

[0052] L is the calculated length of the inclined working surface, D1 is the dip length, s1 and s2 are the offset distances of the inflection points up and down the hill, respectively, and θ0 is the propagation angle of mining influence;

[0053] l1 is the calculated length of limited mining in the strike direction, l1 = D3-s3-s4, D3 is the strike length, s3 and s4 are the offset distances of the left and right turning points of the strike direction respectively;

[0054] The initial probability integral method prediction parameters are set and substituted into formula (7) for simulation calculation. The curve fitting method is used to continuously invert until accurate probability integral method prediction parameters are obtained. The initial probability integral method prediction parameters include the subsidence coefficient q, the horizontal movement coefficient b, the uphill and downhill turning point offset distances s1 and s2, the left and right turning point offset distances s3 and s4, the main influence angle tangent value tanβ, and the mining influence propagation angle θ0.

[0055] In the above-mentioned DInSAR and leveling data fusion surface movement parameter inversion method based on surface subsidence velocity identification, in step (5), the mining size, mining depth, and coal thickness information of the planned working face of the mine are collected, and the probability integral method prediction parameters obtained in step (4) are used to predict the surface subsidence after the working face is mined;

[0056] The final subsidence curve predicted by the probability integral method should fit the DInSAR monitoring measurement results in the non-decoherent area; in the decoherent area, the maximum subsidence value should be close to the maximum subsidence value calculated by the level feature point observation data.

[0057] The technical solution of the present invention achieves the following beneficial technical effects:

[0058] 1. The present invention discloses a method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification. The invention uses the data of the non-decoherent area monitored by DInSAR as a benchmark, and the decoherent area is combined with the ground leveling data, and the probability integral method is used to invert and obtain the probability integral method prediction parameters reflecting the surface subsidence characteristics of the area, and then predict the surface subsidence deformation according to the underground mining plan. This method can set up only a small number of leveling observation points in the large deformation area of ​​the subsidence basin during surface subsidence monitoring, and then with the help of DInSAR non-decoherent monitoring data, it can invert and obtain more accurate probability integral method prediction parameters, which can greatly reduce the workload of ground leveling measurement without causing loss of monitoring accuracy.

[0059] 2. When using this method to analyze the law of surface subsidence, it avoids the time-consuming and labor-intensive leveling measurement, which is too costly, and avoids the problem of DInSAR monitoring data losing coherence due to rapid deformation of the mining area. The probability integral method is used to invert the optimal probability integral method prediction parameters based on the DInSAR monitoring data of the non-losing coherence area and the measurement data of the leveling feature points in the losing coherence area, and predict the settlement information of any point on the surface during the entire mining process of the working face.

[0060] 3. Through the spatial distribution function model V (x,y,t)The sinking velocity of a certain point on the surface during the entire mining process of the working face is calculated and compared with the maximum theoretical deformation rate monitored by DInSAR, and the DInSAR incoherence and non-incoherence areas are divided; the incoherence area is controlled by the leveling feature point data, and the non-incoherence area is controlled by the DInSAR monitoring data. The probability integral method is used to predict the subsidence of the entire process in order to achieve the control process of mining subsidence in the entire mining area, providing theoretical guidance for in-depth research on the law of surface subsidence. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 The fusion flow chart of InSAR data and leveling data of the present invention;

[0062] Figure 2 Relationship diagram between maximum sinking speed and working face distribution position;

[0063] Figure 3 Schematic diagram of the division of decoherent regions and non-decoherent regions;

[0064] Figure 4 The probability integral method is used to fit the subsidence value diagram on the dip profile at a certain moment. DETAILED DESCRIPTION

[0065] InSAR monitoring has many advantages such as high precision, wide range, continuous monitoring, all-day and all-weather operation, easy implementation, low cost, and area monitoring. However, when the surface deformation speed of the monitored mining area is too fast, decoherence will occur, resulting in unreliable monitoring data accuracy. Although the data accuracy of leveling observation is more reliable, this method is based on point observation, the amount of observation data is small, the observation work is time-consuming and laborious, and the observation cost is too high.

[0066] The present application discloses a method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification. The method uses the data of the non-decoherent area monitored by DInSAR as a benchmark, and the decoherent area is combined with ground leveling data, and the probability integral method is used to invert and obtain the probability integral method prediction parameters reflecting the surface subsidence characteristics of the area, and then the surface subsidence deformation is predicted according to the underground mining plan.

[0067] The specific probability integral prediction method flow chart is as follows: Figure 1 shown.

[0068] The first step is to find the formula for calculating the sinking velocity at any point on the surface at any time, collect leveling data in the target monitoring area, obtain the sinking velocity of a certain point on the surface, and calculate the sinking velocity of any point on the surface at any time during the working face mining process.

[0069] The formula for calculating the sinking speed of a surface point in the subsidence area at any time is:

[0070] v(t)=W m CK -ct (1-e -ct ) k-1 (1);

[0071] In formula (1), W m is the maximum subsidence value after the subsidence of a certain point on the surface stabilizes, t represents any time after coal seam mining, c is the lithology parameter, and k is the parameter to be fitted;

[0072] The spatial distribution function formula of the sinking velocity at any point on the surface of the subsidence area (such as Figure 2 shown) is:

[0073]

[0074] In formula (2), v max is the maximum sinking speed of the surface point in the subsidence area during the subsidence process; x is the normal distance from the boundary of the working face (opening cut or stop mining line) in the strike direction, y is the normal distance from the boundary of the working face (transport tunnel or return air tunnel) in the dip direction, and l is the maximum sinking speed lag distance;

[0075] a is the morphological parameter of the main section of the strike, v x is the sinking velocity at any point on the strike main section;

[0076] b is the morphological parameter of the main section, according to the empirical formula Sure, is the average mining depth;

[0077] Substituting formula (1) into formula (2), we can obtain the formula for the sinking velocity at any point on the surface of the subsidence area at any time:

[0078]

[0079] Collect the leveling data in the target monitoring area and calculate the sinking speed of a certain leveling point n at time t:

[0080]

[0081] In formula (4): W n (t) and W n (t-1) represents the sinking value of point n at the tth moment and the t-1th moment respectively; Δt is the number of days between two observations; W n (t) = H t -H0, where H t and H0 are the elevations of point n at time t and the first observation respectively.

[0082] According to the measured data of ground leveling, the sinking velocity value of a certain leveling point during the mining process is used to invert the parameters c and k by formula (1); then, according to the position of the surface point (x, y) that needs to be calculated, the sinking velocity of any surface point at any time in the subsidence area during the whole mining process of the working face is calculated by formula (3).

[0083] The second step is to collect SAR data within the target monitoring area, perform DInSAR data processing, obtain the vertical deformation of the surface points, and determine the theoretical value of the maximum vertical deformation velocity between two adjacent images.

[0084] Collect SAR data within the target monitoring area. The DInSAR data processing method is as follows:

[0085] The dual-track method is used for data processing. The interference phase diagram is obtained by interfering two SAR images: the orbit error phase, atmospheric delay phase and random noise phase in the interference phase are ignored; the flat ground phase is eliminated by removing the flat ground effect and calculating the interference baseline; the terrain phase is eliminated by simulating with the digital elevation model DEM, and the surface deformation phase is finally retained;

[0086] According to the surface deformation phase, the surface vertical deformation R0 is inverted using equation (5):

[0087]

[0088] In formula (5), Δφ is the phase difference of the image pair, λ is the wavelength of the SAR radar signal, and θ is the angle between the radar line of sight and the vertical direction;

[0089] The maximum vertical deformation velocity theoretical value of the satellite used is determined according to the maximum vertical deformation velocity theoretical value that can be monitored between two adjacent images (6):

[0090]

[0091] In formula (6), Δt represents the time interval between two adjacent image data, λ / 2 represents half wavelength, and θ is the incident angle of the SAR radar.

[0092] The third step is to calculate the sinking speed V of a certain point on the surface during the whole mining process of the working face. (x,y,t) The theoretical value V of the maximum vertical deformation velocity that can be monitored between two adjacent images of the DInSAR monitoring data at this point SAR By comparison, it is determined whether the DInSAR monitoring data of the point is incoherent during the whole process of working face mining, and then the incoherent area and non-incoherent area of ​​the DInSAR monitoring data of the working face mining subsidence are determined (such as Figure 3 shown).

[0093] During the whole process of mining subsidence of the working face, the sinking velocity of a certain point on the surface is obtained by formula (3) and compared with the theoretical value V of the maximum vertical deformation velocity that can be monitored by DInSAR (obtained by formula (6)).

[0094] If the sinking velocity V of a point on the surface (x,y,t) During the whole mining process of the working face, the maximum vertical deformation velocity V that can be monitored by DInSAR is exceeded once or multiple times. SAR When , the DInSAR monitoring data of this point loses coherence;

[0095] If the sinking velocity V of a point on the surface (x,y,t) During the whole process of mining, the maximum vertical deformation velocity V that can be monitored by DInSAR was not exceeded. SAR There is no decoherence in the InSAR data at this point.

[0096] According to the above method, the surface subsidence basin can be divided into two parts based on DInSAR monitoring data, including the peripheral non-decoherent area and the internal decoherent area (such as Figure 3 shown).

[0097] The fourth step is to invert the DInSAR monitoring data in the non-decoherent area and the ground leveling data in the decoherent area to obtain the predicted parameters using the probability integral method.

[0098] In the non-decoherent area, the vertical deformation monitored by DInSAR is selected as the measured data, and in the decoherent area, the observation data of the leveling feature points (such as the maximum subsidence value and the inflection point) are selected as the measured data. The formula for predicting the subsidence of any point in the probability integral model is:

[0099]

[0100] In formula (7), (x, y) is the coordinate of any point on the surface, W (x,y) is the subsidence value of any point; W0 is the maximum subsidence value of the surface, W0 = m × q × cos α (m is the coal thickness, q is the subsidence coefficient, α is the coal seam inclination), erf is the probability integral function, its form is u is the integral parameter, r is the main influence radius, is the average mining depth; tanβ is the tangent value of the main influencing angle; r1=H1 / tanβ; r2=H2 / tanβ; r1 is the main influencing radius in the uphill direction; r2 is the main influencing radius in the downhill direction; H1 represents the mining depth of the working face in the uphill direction, and H2 represents the mining depth of the working face in the downhill direction;

[0101] L is the calculated length of the inclined working surface, D1 is the dip length, s1 and s2 are the offset distances of the inflection points of the upper and lower hills, and θ0 is the propagation angle of mining influence. l1 is the calculated length of strike-limited mining, l1=D3-s3-s4, D3 is the strike length, s3 and s4 are the offset distances of the left and right inflection points of the strike, respectively.

[0102] The initial probability integral method prediction parameters (including the subsidence coefficient q, the offset distances s1 and s2 of the uphill and downhill turning points, the offset distances s3 and s4 of the left and right turning points of the strike direction, the main influence angle tangent tanβ, and the mining influence propagation angle θ0) are set and substituted into formula (6) for simulation calculation. The accuracy error between the calculation results and the measured data is compared, and the inversion is continuously performed until the accurate probability integral method prediction parameters are obtained.

[0103] Step 5: Use the probability integral method to predict surface subsidence and deformation according to the mine mining plan

[0104] Collect information such as mining size, mining depth, coal thickness, etc. of the planned working face of the mine, use the probability integral method obtained in step 4 to estimate the parameters, and predict the subsidence of the surface after the working face is mined.

[0105] The final subsidence curve predicted by the probability integral method should fit the DInSAR monitoring measurement results in the non-decoherent area; in the decoherent area, the maximum subsidence value should be close to the maximum subsidence value calculated by the level feature point observation data.

[0106] When using this method to analyze the law of surface subsidence, it avoids the time-consuming and labor-intensive leveling measurement, and the high cost, while avoiding the problem of DInSAR monitoring data losing coherence due to rapid deformation of the mining area. The probability integral method is used to invert the optimal probability integral method prediction parameters based on the DInSAR monitoring data of the non-losing coherence area and the measurement data of the leveling feature points in the losing coherence area, and the settlement information of any point on the surface during the entire mining process of the working face is predicted.

[0107] Specific application examples:

[0108] Calculate the sinking velocity of a working face on the main dip section (with InSAR data of the working face during mining and leveling feature point data on the main dip section).

[0109] Select the sinking velocity values ​​at several moments at the maximum sinking value, and use formula (1) to invert parameters c and k. After obtaining parameters c and k, select the 214th day from the mining of the working face (the maximum sinking velocity measured by leveling) as the time parameter t. Since this point is located on the main dip profile, the x value is 0. The surface sinking velocity V at different y distances is calculated by formula (3): (x,y,t) .

[0110] The observation data of the Sentinel-1 satellite are used as the fitting data of the probability integration method in the non-decoherence region. The observation revisit period of the Sentinel-1 satellite is 12 days, the wavelength of the C band is 5.6 cm, and the radar incident angle is about 41.6°. Substituting it into formula (6), the critical decoherence deformation rate is calculated to be 1.72 mm / d.

[0111] V (x,y,t) The obtained sinking velocity was compared with the critical incoherence deformation rate, and finally the critical distance (distance from the goaf) y for the incoherence of InSAR data in the dip direction was determined to be 300m.

[0112] Using formula (7), the probability integral method predicts that the parameter q is 0.33, the outer turning deviation is 0.12, the mining impact propagation angle coefficient is 0.8, the horizontal movement coefficient is 0.3, and the main impact angle tangent tanβ is 1.3. The probability integral method is used to fit the InSAR monitoring data curve. Figure 4 ,from Figure 4 It can be seen that in the decoherent area, the subsidence value obtained by using the InSAR data is very different from the subsidence value data measured by leveling. The subsidence value data obtained by using the probability integral fitting algorithm of the present application is highly consistent with the subsidence value data measured by leveling.

[0113] Obviously, the above embodiments are merely examples for the purpose of clear explanation, and are not intended to limit the implementation methods. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the implementation methods here. The obvious changes or modifications derived therefrom are still within the scope of protection of the claims of this patent application.

Claims

1. A method for inverting surface movement parameters based on DInSAR and leveling data fusion based on surface subsidence velocity identification, characterized in that: The steps include: (1) Based on the sinking velocity of any point on the surface at any time in the subsidence area and the leveling data of the observation station in the target monitoring area, the sinking velocity V of any point on the surface at any time during the mining process of the working face is calculated. (x,y,t) ; (2) Collect SAR data within the target monitoring area, perform DInSAR data processing, obtain the vertical deformation of the surface points, and determine the maximum vertical deformation velocity theoretical value V between two adjacent images. SAR ; (3) The sinking velocity V of a point on the ground during the whole mining process is (x,y,t) The theoretical value V of the maximum vertical deformation velocity that can be monitored between two adjacent images of the DInSAR monitoring data at this point SAR By comparison, it is determined whether the DInSAR monitoring data of the point is incoherent during the whole process of mining at the working face, and then the incoherent area and non-incoherent area of ​​the DInSAR monitoring data of the mining subsidence area of ​​the working face are determined; (4) Combine the DInSAR monitoring data of the non-decoherent area and the ground leveling data of the decoherent area to invert and obtain the estimated parameters using the probability integral method; (5) According to the mine mining plan, the probability integral method is used to predict surface subsidence and deformation.

2. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 1 is characterized in that: In step (1), the following steps are included: (1-1) Using the formula for calculating the subsidence velocity of a surface point in the subsidence area at any time and the formula for the spatial distribution function of the subsidence velocity of any surface point in the subsidence area, the formula for the subsidence velocity of any surface point in the subsidence area at any time is obtained; (1-2) Collecting leveling data within the target monitoring area and obtaining the sinking velocity of a certain leveling point; (1-3) Using the sinking velocity value of a certain leveling measurement point, the parameters in the sinking velocity formula at any point on the surface of the subsidence area at any time are inverted to obtain the sinking velocity at any point on the surface of the subsidence area at any time.

3. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 2 is characterized in that: In step (1-1), the formula for calculating the sinking velocity of a surface point in the subsidence area at any time is: v(t)=W m cke -ct (1-e -ct ) k-1 (1); In formula (1), W m is the maximum subsidence value after the subsidence of a certain point on the surface stabilizes, t represents any time after coal seam mining, c is the lithology parameter, and k is the parameter to be fitted; The formula for the spatial distribution function of the sinking velocity at any point on the surface of the subsidence area is: In formula (2), v max is the maximum sinking velocity of the surface point in the subsidence area during the subsidence process; x is the normal distance from the boundary of the working face in the strike direction, y is the normal distance from the boundary of the working face in the dip direction, and l is the maximum sinking velocity lag distance; a is the morphological parameter of the main section of the strike, v x is the sinking velocity at any point on the strike main section; b is the morphological parameter of the main section, according to the empirical formula Sure, is the average mining depth; Substituting formula (1) into formula (2), we can obtain the formula for the sinking velocity at any point on the surface of the subsidence area at any time:

4. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 3 is characterized in that: In step (1-2), the leveling data in the target monitoring area is collected to calculate the sinking speed of a certain leveling point n at time t: In formula (4): W n (t) and W n (t-1) represents the sinking value of point n at the tth moment and the t-1th moment respectively; Δt is the number of days between two observations; W n (t) = H t -H0, where H t and H0 are the elevations of point n at time t and the first observation respectively.

5. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 4 is characterized in that: In steps (1-3), According to the measured data of ground leveling, the sinking velocity value of a certain leveling point during the mining process is used to invert the parameters c and k by formula (1); then, according to the position of the surface point (x, y) that needs to be calculated, the sinking velocity of any surface point at any time in the subsidence area during the whole mining process of the working face is calculated by formula (3).

6. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 1, characterized in that: In step (2), SAR data within the target monitoring area is collected, and the DInSAR data processing method is as follows: The dual-track method is used for data processing. The interference phase diagram is obtained by interfering two SAR images. The orbit error phase, atmospheric delay phase and random noise phase in the interference phase are ignored. The flat ground phase is eliminated by removing the flat ground effect and calculating the interference baseline. The terrain phase is eliminated through simulation with the digital elevation model DEM, and the surface deformation phase is finally retained; According to the surface deformation phase, the surface vertical deformation R0 is inverted using equation (5): In formula (5), Δφ is the phase difference of the image pair, λ is the wavelength of the SAR radar signal, and θ is the angle between the radar line of sight and the vertical direction; The maximum vertical deformation velocity theoretical value of the satellite used is determined according to the maximum vertical deformation velocity theoretical value that can be monitored between two adjacent images (6): In formula (6), Δt represents the time interval between two adjacent image data, λ / 2 represents half wavelength, and θ is the incident angle of the SAR radar.

7. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 6 is characterized in that: In step (3), during the whole mining process of the working face, the subsidence velocity at any point on the surface at any time is calculated by formula (3), which is equal to the theoretical value V of the maximum vertical deformation velocity that can be monitored by DInSAR at that point. SAR Compare and determine whether decoherence occurs.

8. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 7 is characterized in that: The comparison method is: If the sinking velocity V of a point on the surface (x,y,t) During the whole process of mining subsidence, the maximum vertical deformation velocity V that can be monitored by DInSAR is exceeded once or multiple times. SAR When , the DInSAR monitoring data of this point loses coherence; If the sinking velocity V of a point on the surface (x,y,t) During the whole process of mining, the maximum vertical deformation velocity V that can be monitored by DInSAR was not exceeded. SAR When , the InSAR data at this point does not lose coherence; According to the above method, the surface subsidence basin is divided into two parts based on the DInSAR monitoring data, including the peripheral non-decoherent area and the internal decoherent area.

9. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 8, characterized in that: In step (4), the vertical deformation monitored by DInSAR in the non-decoherent area is selected as the measured data, and the surface subsidence data observed by the leveling point in the decoherent area is selected as the measured data, and the probability integral method prediction parameters are inverted respectively; The formula for predicting sinking at any point in the probability integral model is: In formula (7), (x, y) is the coordinate of any point on the surface, W (x,y) is the sinking value of any point; W0 is the maximum surface subsidence value, W0 = m × q × cos α, m is the coal thickness, q is the subsidence coefficient, and α is the coal seam inclination; erf is the probability integral function, which is in the form of u is the integral parameter, r is the main influence radius, is the average mining depth; tanβ is the tangent value of the main influencing angle; r1=H1 / tanβ; r2=H2 / tanβ; r1 is the main influencing radius in the uphill direction; r2 is the main influencing radius in the downhill direction; H1 represents the mining depth of the working face in the uphill direction, and H2 represents the mining depth of the working face in the downhill direction; L is the calculated length of the inclined working surface, D1 is the dip length, s1 and s2 are the offset distances of the inflection points up and down the hill, respectively, and θ0 is the propagation angle of mining influence; l1 is the calculated length of limited mining in the strike direction, l1 = D3-s3-s4, D3 is the strike length, s3 and s4 are the offset distances of the left and right turning points of the strike direction respectively; The initial probability integral method prediction parameters are set and substituted into formula (7) for simulation calculation. The curve fitting method is used to continuously invert until accurate probability integral method prediction parameters are obtained. The initial probability integral method prediction parameters include the subsidence coefficient q, the horizontal movement coefficient b, the uphill and downhill turning point offset distances s1 and s2, the left and right turning point offset distances s3 and s4, the main influence angle tangent value tanβ, and the mining influence propagation angle θ0.

10. The method for inverting surface movement parameters by fusing DInSAR and leveling data based on surface subsidence velocity identification according to claim 9, characterized in that: In step (5), the mining size, mining depth, and coal thickness information of the planned working face of the mine are collected, and the probability integral method prediction parameters obtained in step (4) are used to predict the subsidence of the surface after the working face is mined; The final subsidence curve predicted by the probability integral method should fit the DInSAR monitoring measurement results in the non-decoherent area; in the decoherent area, the maximum subsidence value should be close to the maximum subsidence value calculated by the level feature point observation data.

Images