CN104062660A - Mining area earth surface time sequence deformation monitoring method based on time domain discrete InSAR interference pair - Google Patents

Abstract

The invention discloses a time-series surface deformation monitoring method in a mining area based on a time-domain discrete InSAR interference pair, which obtains a time-domain discrete InSAR interference pair that does not cover the entire time-series process; obtains high coherence points in the time-domain discrete InSAR interference pair; establishes The relationship equations between the dynamic subsidence model parameters of the mining area and the unwrapped phase of the high-coherence points of the time-domain discrete InSAR interferometric pair are used to estimate the parameters of the dynamic subsidence model of the high-coherence points in the mining area; finally, the parameters of the mining area can be estimated at any time The time-series deformation of the ground surface, thus realizing the time-series deformation monitoring of the mining area surface based on time-domain discrete InSAR interferometric pairs. The present invention realizes the use of discrete InSAR interferometric pairs that do not cover the time-series process to monitor the time-series deformation of the surface of the mining area. The concept is ingenious, the process is simple, the monitoring results are accurate and effective, the application prospect of InSAR technology is greatly broadened, and the cost and technology of time-series deformation monitoring in the mining area are reduced. limit.

Description

A kind ofly based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring method

Technical field

The present invention relates to a kind ofly based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring method.

Background technology

It is a kind ofly can obtain centimetre even remote sensing technology of grade deformation of earth's surface that interfering synthetic aperture radar is measured (Interferometric Synthetic Aperture Radar is called for short InSAR).Its ultimate principle is exactly by two width or above synthetic-aperture radar (Synthetic Aperture Radar is called for short SAR) satellite image being carried out to differential interferometry processing, extracting centimetre even millimetre-sized radar line of sight direction deformation values from phase differential.For traditional Geodetic Technique, the advantage such as that InSAR has is round-the-clock, continuous space covering, high precision, low cost.But, traditional difference interfering synthetic aperture radar is measured (differential interferometric SAR, D-InSAR) be subject to the factor impacts such as dephasing pass, space, time dephasing pass and atmospheric disturbance, and cannot obtain the InSAR sequential Deformation Field on earth's surface.In order to overcome these restrictions, some sequential InSAR deformation monitoring methods are suggested, as permanent scatterer interferometry (permanent scatter interferometry, PS-InSAR) and Small Baseline Subset technology (small baseline subset interferometry, SBAS-InSAR).At present sequential InSAR technology be widely used with the deformation monitoring such as urban surface, volcano, landslide, underground mining in.

Yet, because traditional sequential InSAR technology (PS-InSAR or SBAS-InSAR) needs a large amount of SAR data, and require InSAR to interfere covering program process when whole.For the front and back SAR image causing because of reasons such as Ground Deformation are excessive, space-time dephasing pass is serious, cannot interfere the time domain discrete InSAR forming to interfere right, traditional sequential InSAR technology is helpless.Yet this phenomenon is comparatively common in the deformation monitoring of mining area.Therefore, not only cost is higher, data demand is harsh to utilize traditional sequential InSAR technical monitoring mining area surface sequential deformation, and is difficult to the sequential deformation of the large magnitude in monitoring mining area, and it has restricted the range of application of InSAR technology in mining area greatly.Therefore; how to utilize time domain discrete InSAR to interfere the deformation of monitoring mining area surface high precision sequential; for widening InSAR application space; reduce mining area sequential deformation monitoring cost, and instruct mining area safety production, the geologic hazard of early warning mining area surface and ecological environmental protection all to play an important role with this.

Summary of the invention

The object of the present invention is to provide and a kind ofly based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring method, it has overcome, and traditional sequential InSAR method monitoring cost is high, data demand is harsh, and cannot monitor the defects such as large magnitude Deformation Field.

Based on time domain discrete InSAR, interfere a right mining area surface sequential deformation monitoring method, comprise following step:

Step 1: the time domain discrete InSAR interference of obtaining program process while not covering whole SAR image is right;

By all SAR images that cover mining area to be monitored according to time order and function, utilize difference interfering synthetic aperture radar to measure D-InSAR and carry out differential interferometry between two, obtain the reconciliation of coherence map group and twine figure group, from coherence map group, find out the InSAR interference that cannot interfere right, and reject corresponding with it coherence map and conciliate and twine phase diagram, thereby the InSAR that obtains the time domain discrete of program process while not covering interfere to and corresponding coherence map group conciliate and twine figure group;

Described coherence map, is the foundation of evaluating two width SAR image similarity degrees, in interfering processing, generates;

Described cannot the interference, refer to that the coherence of coherence map is less than coherence's threshold value of setting;

Step 2: obtain the high coherent point that time domain discrete InSAR interferes centering;

While never covering according to coherence's threshold value of setting, the InSAR of the time domain discrete of program process interferes the high coherent point of interfering centering to extracting each time domain discrete InSAR in corresponding coherence map group;

Coherence's threshold value that the coherence map that described high coherent point is all time domain discrete is all set the coherence of this point;

Step 3: utilize low-pass filtering to weaken atmosphere delay and the noise at high coherent point place, and ignore and move horizontally the contribution to deformation to radar line of sight, the solution of high coherent point (i, j) twines phase place δ φ and is:

$\mathrm{\δ\φ}(i,j)=\frac{4\mathrm{\π}}{\mathrm{\λ}\cos \mathrm{\θ}(i,j)}[W(t_B,i,j)-W(t_A,i,j)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}\mathrm{\Δh}(i,j)}{r\sin \mathrm{\θ}(i,j)}$

Wherein, λ is radar wavelength, the radar incident angle that θ (i, j) is high coherent point, and r is that radar satellite is apart from the distance of target, B _⊥be the vertical parallax length of two width SAR images, t _b, t _abe respectively the acquisition time of two width SAR images, above-mentioned parameter all directly obtains from the header file of corresponding SAR image;

W is the surface subsidence value of high coherent point, interferes right solution to twine phase diagram obtain from time domain discrete InSAR, and Δ h is the elevation residual error in high coherent point, is coefficient to be asked;

Step 4: choosing mining area dynamic settling model is W (Δ t)=f (Δ t, P);

In formula, Δ t is with respect to mining area initial settlement moment t ₀interval time; F is dynamic settling Model Mapping function; P is model solve for parameter, and number is numP;

Step 5: set up dynamic settling model solve for parameter P and solution and twine the relation equation of phase place;

$\mathrm{\δ\φ}(i,j)=\frac{4\mathrm{\π}}{\mathrm{\λ}\cos \mathrm{\θ}(i,j)}[f({\mathrm{\Δt}}_B,P,i,j)-f(\mathrm{\Δ}{t}_A,P,i,j)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}\mathrm{\Δh}(i,j)}{r\sin \mathrm{\θ}(i,j)}$

Wherein, Δ t _b=t _b-t ₀with Δ t _a=t _a-t ₀for the SAR image time after proofreading and correct;

Step 6: utilization is more than or equal to numP+1 time domain discrete InSAR and interferes right solution twine phase place and corresponding radar wavelength λ, incidence angle θ, oblique distance r and respectively interfere right vertical parallax length B _⊥respectively in the formula described in substitution step 5, by a plurality of equations simultaneousness system of equations that obtain, calculate the solve for parameter P of dynamic settling model and the elevation residual delta h of high coherent point, described time domain discrete InSAR interferes right solution to twine solution that phase place obtains from step 1 to twine phase diagram and obtain; The mining area dynamic settling model that the solve for parameter P substitution step 4 of dynamic settling model is chosen, calculates the earth's surface sequential sedimentation of any time, realizes and interferes right mining area sequential deformation monitoring based on time domain discrete InSAR.

The selected mining area of described step 4 dynamic model comprises any one in Knothe model, Gompertz model, Logistic model, Richards model or Weibull model.

In described step 6, during solving model parameter P to be assessed, utilize the solve for parameter globally optimal solution of Genetic algorithm searching dynamic settling model, thus the mining area dynamic settling model that substitution step 4 is chosen.

In described step 2, coherence's threshold value is 0.2-0.4.

Beneficial effect

The invention provides and a kind ofly based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring method, the SAR image that covers mining area to be monitored is generated to interferogram between two according to time order and function order, reject the interference that mining area surface to be monitored cannot interfere right, stay the InSAR that does not cover the time domain discrete of program process when whole interfere right; Obtain the high coherent point that time domain discrete InSAR interferes centering; Set up mining area dynamic settling model parameter and time domain discrete InSAR and interfere right high coherent point solution to twine the relation equation group of phase place, based on this system of equations, estimate the dynamic settling model parameter of the high coherent point in mining area; Finally use this parameter can estimate the sequential deformation of any time mining area surface, thereby realized based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring.The present invention has realized utilizing and has not covered the discrete InSAR of program process when whole and interfere the sequential deformation to monitoring mining area surface; be skillfully constructed; process is simple; monitoring result accurate and effective; broken through that traditional sequential InSAR technical requirement data volume is large, InSAR interferes covering whole time span and cannot monitor the limitations such as quick distortion; greatly widened InSAR technology application prospect, reduced mining area sequential deformation monitoring cost and technical limitation, for mining area surface geo-hazard early-warning assessment and ecological environmental protection provide important technical support.

Accompanying drawing explanation

Fig. 1 time domain discrete InSAR interferes forming schematic diagram;

The schematic flow sheet of Fig. 2 described method of the present invention;

Fig. 3 simulation utilizes time domain discrete InSAR to interfere the sinking field to high coherent point place, the earth's surface obtaining;

Fig. 4 utilizes time domain discrete InSAR to interfere the sequential deformation at the high coherent point place to estimating, in figure, all to take for the first scape image capturing time be benchmark the time;

Fig. 5 utilizes time domain discrete InSAR to interfere the statistic histogram with simulation deformation differences to the sequential deformation of monitoring.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described further.

As shown in Figure 2, be the schematic flow sheet of the method for the invention, a kind ofly based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring method, comprise following step:

Step 1: the time domain discrete InSAR interference of obtaining program process while not covering whole SAR image is right;

By all SAR images that cover mining area to be monitored according to time order and function, utilize difference interfering synthetic aperture radar to measure D-InSAR and carry out differential interferometry between two, obtain coherence map group and solution and twine figure group, from coherence map group, find out the InSAR interference that cannot interfere right, and from rejecting the coherence map group answer in contrast, conciliate and twine figure group, thereby the InSAR that obtains the time domain discrete of program process while not covering interfere to and corresponding coherence map group conciliate and twine figure group;

Described coherence map, is the foundation of evaluating two width SAR image similarity degrees, in interfering processing, generates;

Described cannot the interference, refer to that the coherence of coherence map is less than coherence's threshold value of setting;

Step 2: obtain the high coherent point that time domain discrete InSAR interferes centering;

While never covering according to coherence's threshold value of setting, the InSAR of the time domain discrete of program process interferes the high coherent point of interfering centering to extracting each time domain discrete InSAR in corresponding coherence map group;

Coherence's threshold value that the coherence map that described high coherent point is all time domain discrete is all set the coherence of this point;

Step 3: utilize low-pass filtering to weaken atmosphere delay and the noise at high coherent point place, and ignore and move horizontally the contribution to deformation to radar line of sight, the solution of high coherent point (i, j) twines phase place δ φ and is:

$\mathrm{\δ\φ}(i,j)=\frac{4\mathrm{\π}}{\mathrm{\λ}\cos \mathrm{\θ}(i,j)}[W(t_B,i,j)-W(t_A,i,j)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}\mathrm{\Δh}(i,j)}{r\sin \mathrm{\θ}(i,j)}$

Wherein, λ is radar wavelength, the radar incident angle that θ (i, j) is high coherent point, and r is that radar satellite is apart from the distance of target, B _⊥be the vertical parallax length of two width SAR images, t _b, t _abe respectively the acquisition time of two width SAR images, above-mentioned parameter all directly obtains from the header file of corresponding SAR image;

W is the surface subsidence value of high coherent point, interferes right solution to twine phase diagram obtain from time domain discrete InSAR, and Δ h is the elevation residual error in high coherent point, is coefficient to be asked;

Step 4: choosing mining area dynamic settling model is W (Δ t)=f (Δ t, P);

In formula, Δ t is with respect to mining area initial settlement moment t ₀interval time; F is dynamic settling Model Mapping function; P is model solve for parameter, and number is numP;

Step 5: set up dynamic settling model solve for parameter P and solution and twine the relation equation of phase place;

$\mathrm{\δ\φ}(i,j)=\frac{4\mathrm{\π}}{\mathrm{\λ}\cos \mathrm{\θ}(i,j)}[f({\mathrm{\Δt}}_B,P,i,j)-f(\mathrm{\Δ}{t}_A,P,i,j)]+\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}\mathrm{\Δh}(i,j)}{r\sin \mathrm{\θ}(i,j)}$

Wherein, Δ t _b=t _b-t ₀with Δ t _a=t _a-t ₀for the SAR image time after proofreading and correct;

Step 6: utilization is more than or equal to numP+1 time domain discrete InSAR and interferes right solution twine phase place and corresponding radar wavelength λ, incidence angle θ, oblique distance r and respectively interfere right vertical parallax length B _⊥respectively in the formula described in substitution step 5, by a plurality of equations simultaneousness system of equations that obtain, calculate the solve for parameter P of dynamic settling model and the elevation residual delta h of high coherent point, described time domain discrete InSAR interferes right solution to twine solution that phase place obtains from step 1 to twine phase diagram and obtain; The mining area dynamic settling model that the solve for parameter P substitution step 4 of dynamic settling model is chosen, calculates the earth's surface sequential sedimentation of any time, realizes and interferes right mining area sequential deformation monitoring based on time domain discrete InSAR.

As shown in Figure 1, suppose that the available SAR image that covers mining area to be monitored has n+1 width, its acquisition time is respectively (t ₀, t ₁..., t _n) ((t wherein ₀<t ₁< ... <t _n)).Because the factors such as speed of deformation is very fast, space-time dephasing is closed, atmosphere delay is serious cause part SAR image and follow-up image to interfere, as shown in the black fork in Fig. 1, thereby formed, do not cover when whole the discrete InSAR of program process interfere rightly on time, this phenomenon is comparatively common in mining area surface deformation.Yet for such data, except adopting the average velocity that precision is lower, traditional sequential InSAR technology cannot therefrom be obtained accurate earth's surface, mining area sequential Deformation Field.Therefore, the present invention is directed to this phenomenon has proposed to interfere right mining area surface sequential deformation monitoring method based on time domain discrete InSAR.

Now suppose that it is t by the time that k width solution twines phase diagram _aand t _btwo width SAR images obtain, the phase value that this solution twines any high coherent point (i, j) in phase diagram can be expressed as:

$\begin{array}{c}{\mathrm{\&delta;\&phi;}}_k(i,j)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}[\mathrm{LOS}(t_B,i,j)-\mathrm{LOS}(t_A,i,j)]+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;k}\mathrm{\&Delta;h}(i,j)}{r\sin \mathrm{\&theta;}(i,j)}\\ +[{\mathrm{\&phi;}}_{\mathrm{atm}}(t_B,i,j)-{\mathrm{\&phi;}}_{\mathrm{atm}}(t_{A}i,j)]+\mathrm{\&Delta;}{n}_k(i,j)\end{array}---\left(1\right)$

In formula, δ φ _kbe that the solution that k width solution twines phase diagram twines phase place, LOS for observation constantly earth's surface at the deformation values of radar line of sight direction, B _{⊥ k}be that k width solution twines that in phase diagram, to interfere right vertical parallax, Δ h be vertical error, λ is radar wavelength, the oblique distance distance that r is radar, and θ is radar incident angle, φ _atmfor atmosphere delay phase place, Δ n is noise.

Because mining area surface deformation mainly be take sinking W as main, and its be radar line of sight to the main contributions of deformation, therefore, the present invention suppose radar line of sight to move horizontally contribution be negligible, that is: W=LOS/cos θ.For atmosphere delay phase place and the noise in formula (1), can it be weakened by the mode of low-pass filtering.So far, formula (1) can be reduced to:

${\mathrm{\&delta;\&phi;}}_k(i,j)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}[W(t_B,i,j)-W(t_A,i,j)]+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;}\mathrm{\&Delta;h}(i,j)}{r\sin \mathrm{\&theta;}(i,j)}---\left(2\right)$

Mining area surface dynamic settling model is a lot, (these models are an available mapping function W (t)=f (t all more typically Knothe model, Gompertz model, Logistic model, Weibull model etc., P) represent, as previously mentioned, W (t) is t surface subsidence constantly, f is the mapping function of dynamic settling model, and P is model parameter).Because the time t of mining area dynamic settling model is with respect to initial settlement relative time constantly, therefore, dynamic settling model f (t, P) substitution formula (2) is being needed to utilize before the first scape SAR image time t ₀proofread and correct the acquisition time (t of all SAR images ₀, t ₁..., t _n), SAR image time (the Δ t after its correction ₀, Δ t ₁..., Δ t _n)=(0, t ₁-t ₀..., t _n-t ₀).Therefore, the SAR image time Δ t after correction _mtime, the dynamic settling model on earth's surface can be expressed as:

W(Δt _m)＝f(Δt _m,P)m＝0,1,…,n (3)

Formula (3) substitution formula (2) can be obtained:

$\mathrm{\&delta;\&phi;}(i,j)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}[f({\mathrm{\&Delta;t}}_B,P,i,j)-f(\mathrm{\&Delta;}{t}_A,P,i,j)]+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;}\mathrm{\&Delta;h}(i,j)}{r\sin \mathrm{\&theta;}(i,j)}---\left(4\right)$

For high coherent point (i, j), all time domain discrete InSAR interfere all setting up the equation suc as formula (4) at this point.The solve for parameter number (being the number of P) of supposing mining area dynamic settling model is numP, the solve for parameter number in formula (4) is numP+1 (in formula, Δ h is also unknown number), therefore, as long as it is right to available InSAR interference to be more than or equal in theory numP+1, no matter interfere whether discrete in time domain, all can estimate model parameter P.Afterwards, use this parameter and in conjunction with dynamic settling model, can predict the earth's surface sequential sedimentation of any time, thereby realize, based on time domain discrete InSAR, interfere right mining area sequential deformation monitoring.

In order more clearly to understand content of the present invention, it is example that the present invention will be take the conventional dynamic settling model---Logistic model---in a mining area, and verifies validity of the present invention in conjunction with simulated experiment.

First, suppose to have 14 scapes to cover the available SAR image in mining area to be studied, the time interval is between two 46 days.Afterwards, supposing within the SAR image cover time, to have exploited one in mining area to be studied, to adopt be 600 meters deeply, and adopting thick is 2.5m, and length and width are respectively 1150 and 150m, and average exploitation rate is the workplace of 2.5m/ days, and utilizes fast Lagrangian analysis software FLAC ^3Dthe Ground Deformation that the underground mining of having simulated causes.

It is right that the InSAR that utilizes random number generator simulation to generate cannot to interfere interferes, and produces 13 random numbers between 1 to 1800, if certain random number surpasses 1200, thinks that its corresponding InSAR interferes interfering, and the simulation of rejecting its correspondence is sunk.

Then utilize high coherent point position and elevation residual error in random number generator simulation deformation map, and interfere the solution to generating to twine phase transition for sinking the time domain discrete InSAR that utilizes of simulation, thereby obtain the sinking field at high coherent point place, region to be monitored;

Utilize the method for the invention to process SAR image data, concrete steps are as follows:

Step 1: the time domain discrete InSAR interference of obtaining program process while not covering whole SAR image is right;

By all SAR images that cover mining area to be monitored according to time order and function, utilizing difference interfering synthetic aperture radar to measure D-InSAR interferes between two, obtain the reconciliation of coherence map group and twine figure group, from coherence map group, find out the InSAR interference that cannot interfere right, and reject the coherence map answer in contrast, thereby the InSAR that obtains the time domain discrete of program process while not covering interfere to and corresponding coherence map group conciliate and twine figure group;

Find that the 2nd, 7 and 12 phases are right for the InSAR that cannot interfere interferes according to the time order and function coherence map group that differential interferometry obtains between two utilizing 14 scape SAR images, rejected, the InSAR that obtains the time domain discrete of program process while not covering interfere to and corresponding coherence map group conciliate and twine figure group;

Step 2: obtain the high coherent point that time domain discrete InSAR interferes centering;

While never covering according to coherence's threshold value 0.3 of setting, the InSAR of the time domain discrete of program process interferes the high coherent point of interfering centering to extracting each time domain discrete InSAR in corresponding coherence map group;

Interfere the solution of the high coherent point of centering to twine phase diagram each time domain discrete InSAR and be converted to field, as shown in Figure 3;

Step 3: utilize low-pass filtering to weaken atmosphere delay and the noise at high coherent point place, and ignore and move horizontally the contribution to deformation to radar line of sight, the solution of high coherent point (i, j) twines phase place δ φ and is:

$\mathrm{\&delta;\&phi;}(i,j)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}[W(t_B,i,j)-W(t_A,i,j)]+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;}\mathrm{\&Delta;h}(i,j)}{r\sin \mathrm{\&theta;}(i,j)}$

Wherein, λ is radar wavelength, the radar incident angle that θ (i, j) is high coherent point, and r is that radar satellite is apart from the distance of target, B _⊥be the vertical parallax length of two width SAR images, t _b, t _abe respectively the acquisition time of two width SAR images, above-mentioned parameter all directly obtains from the header file of corresponding SAR image;

W is the surface subsidence value of high coherent point, interferes right solution to twine phase diagram obtain from time domain discrete InSAR, and the elevation residual error that Δ h is high coherent point, is coefficient to be asked;

Step 4: choosing mining area dynamic settling model is W (Δ t)=f (Δ t, P); In the present embodiment, select Logistic model;

When utilizing Logistic model description mining area surface dynamic settling process, its expression formula is:

$W\left(\mathrm{\&Delta;t}\right)=\frac{W_0}{1+{\mathrm{ae}}^{-\mathrm{b\&Delta;t}}}---\left(5\right)$

In formula, W ₀for maximum sinking value, a, b is the form parameter of Logistic curve, three is model solve for parameter P, that is: P=[W ₀, a, b], Δ t is with respect to initial settlement correction time constantly.

Step 5: set up mining area dynamic settling model solve for parameter P and solution and twine the relation equation of phase place;

After being replaced with to Logistic model, dynamic settling model in formula (4) can obtain:

${\mathrm{\&delta;\&phi;}}_k(i,j)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}[\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_B}}-\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_A}}]+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;}\mathrm{\&Delta;h}(i,j)}{r\sin \mathrm{\&theta;}(i,j)}---\left(6\right)$

As previously mentioned, δ φ is that solution twines phase place, B _⊥for interfering right vertical parallax, λ is radar wavelength, the oblique distance distance that r is radar, and θ is radar incident angle.These parameters all can be conciliate to twine phase diagram from SAR camera file and be obtained.W ₀for maximum sinking value, a, b is the form parameter of Logistic curve, and Δ h is vertical error, and rear four are solve for parameter.

Step 6: utilization is more than or equal to numP+1 time domain discrete InSAR and interferes right solution twine phase place and corresponding radar wavelength λ, incidence angle θ, oblique distance r and respectively interfere right vertical parallax length B _⊥respectively in the formula (6) described in substitution step 5, by a plurality of equations simultaneousness system of equations that obtain, calculate the solve for parameter P of dynamic settling model and the elevation residual delta h of high coherent point, described time domain discrete InSAR interferes right solution to twine solution that phase place obtains from step 1 to twine phase diagram and obtain; The mining area dynamic settling model that the solve for parameter P substitution step 4 of dynamic settling model is chosen, calculates the earth's surface sequential sedimentation of any time, realizes and interferes right mining area sequential deformation monitoring based on time domain discrete InSAR.

Known according to simulated data, the InSAR of the 2nd, 7,12 phases interferes interfering.Therefore, all the other time domain discrete InSAR interfere the equation suc as formula (6) that all can set up locating in high coherent point (i, j), and are used matrix form to be expressed as:

A _M×2P′ _2×M＝δφ _M×1(7)

In formula: the matrix of coefficients that A is system of equations, M can with time domain discrete InSAR interfere for number (for this simulated experiment, M=10), P ' is unknown number battle array, and δ φ is that solution twines phase value matrix, and its concrete form is as follows:

$B={\left[\begin{array}{cc}\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}& \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;1}}{r_1\sin \mathrm{\&theta;}(i,j)}\\ \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}& \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;3}}{r_3\sin \mathrm{\&theta;}(i,j)}\\ .& .\\ .& .\\ .& .\\ \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos \mathrm{\&theta;}(i,j)}& \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;n-2}}{r\sin \mathrm{\&theta;}(i,j)}\\ \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}\cos (i,j)}& \frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B_{\&perp;n}}{r\sin \mathrm{\&theta;}(i,j)}\end{array}\right]}_{M\&times;2}$ $P^{\&prime;}={\left[\begin{array}{cc}\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_1}}-\frac{W_0(i,j)}{1+a(i,j)e^0}& \mathrm{\&Delta;h}(i,j)\\ \frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_3}}-\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_2}}& \mathrm{\&Delta;h}(i,j)\\ .& .\\ .& .\\ .& .\\ \frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_{n-2}}}-\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_{n-3}}}& \mathrm{\&Delta;h}(i,j)\\ \frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_n}}-\frac{W_0(i,j)}{1+a(i,j)e^{-b(i,j)\mathrm{\&Delta;}{t}_{n-1}}}& \mathrm{\&Delta;h}(i,j)\end{array}\right]}_{2\&times;M}^T$

${\mathrm{\&delta;\&phi;}}_{M\&times;1}={\left[\begin{array}{c}{\mathrm{\&delta;\&phi;}}_1(i,j)\\ {\mathrm{\&delta;\&phi;}}_3(i,j)\\ .\\ .\\ .\\ {\mathrm{\&delta;\&phi;}}_{n-2}(i,j)\\ {\mathrm{\&delta;\&phi;}}_n(i,j)\end{array}\right]}_{M\&times;1}$

In system of equations (7), its unknown number to be asked is model parameter W ₀, a, b and elevation residual delta h.That is to say, if the number of equation surpasses 4 in (7), can solve all unknown numbers.Yet, from formula (7), can find out, each equation is non-linear, and direct solution is difficulty comparatively.Therefore, the present invention utilizes the globally optimal solution of the Genetic algorithm searching equation unknown number in intelligent algorithm.For other the high coherent point in interferogram, all can utilize aforesaid way to estimate each high relevant globally optimal solution of pointing out unknown number.Bring the globally optimal solution of each high coherent point into Logistic model, can estimate the earth's surface sequential deformation of any time of all SAR images, thereby realize, based on time domain discrete InSAR, interfere right mining area surface sequential deformation monitoring.

In order to verify reliability of the present invention, first, bring the correction time of the globally optimal solution of the high coherent point of above-mentioned estimation and SAR image into Logistic model, then estimate the sequential dynamic settling of 14 phases of each point, as shown in Figure 4.Afterwards, the present invention has added up the histogram of difference between the 14 phases sequential sedimentation estimated and the analogue value, and its result as shown in Figure 5.As can be seen from the figure, both differences mainly concentrate near 0, and its average is 0.7mm, and root-mean-square error is 37.3mm.This result shows, the present invention proposes, and based on the discontinuous InSAR of time domain, to interfere obtaining the method for mining area surface sequential deformation be accurately feasible.

Claims (4)

1. A method for monitoring time-series deformation of the mining area surface based on time-domain discrete InSAR interferometry, characterized in that it comprises the following steps: Step 1: Obtain time-domain discrete InSAR interferometric pairs that do not cover the entire SAR image timing process; All the SAR images covering the mining area to be monitored are time-sequenced, and the differential synthetic aperture radar interferometry D-InSAR is used to perform pairwise differential interference to obtain the coherent image group and the unwrapped image group, and find out the InSAR interference that cannot be interfered from the coherent image group Yes, and remove the corresponding coherence map and unwrapped phase map, so as to obtain the time-domain discrete InSAR interferometric pair and the corresponding coherent map group and unwrapped map group that do not cover the timing process; The coherence map is the basis for evaluating the similarity of two SAR images, and is generated during interference processing; The non-interference means that the coherence of the coherence map is less than the set coherence threshold; Step 2: Obtain the highly coherent points in the time-domain discrete InSAR interferometric pair; According to the set coherence threshold, extract the high-coherence points in each time-domain discrete InSAR interferometric pair from the coherence map group corresponding to the time-domain discrete InSAR interferometric pair that does not cover the time series process; The high coherence point is the coherence threshold set for the coherence of all time-domain discrete coherence maps at this point; Step 3: Use low-pass filtering to weaken the atmospheric delay and noise at the high-coherence point, and ignore the contribution of horizontal movement to the radar line-of-sight deformation. The unwrapped phase δφ of the high-coherence point (i, j) is: $\mathrm{<span\; class="notranslate">\; \&delta;\&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 4</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}\mathrm{<span\; class="notranslate">\; cos</span>}\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; W</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; B</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; W</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 4</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}\frac{{\mathrm{<span\; class="notranslate">\; B</span>}}_{<span\; class="notranslate">\; \&perp;</span>}\mathrm{<span\; class="notranslate">\; \&Delta;h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; sin</span>}\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}$ Among them, λ is the radar wavelength, θ(i,j) is the radar incidence angle of the high coherence point, r is the distance from the radar satellite to the target, B _⊥ is the vertical baseline length of the two SAR images, t _B and t _A are respectively The acquisition time of the two SAR images, the above parameters are directly obtained from the header files of the corresponding SAR images; W is the surface subsidence value of the high coherence point, obtained from the unwrapped phase map of the time-domain discrete InSAR interferometric pair, Δh is the elevation residual at the high coherence point, and is the coefficient to be obtained; Step 4: Select the dynamic settlement model of the mining area as W(Δt)=f(Δt,P); In the formula, Δt is the interval time relative to the initial settlement moment t ₀ of the mining area; f is the dynamic settlement model mapping function; P is the model parameter to be estimated, and the number is numP; Step 5: Establish the relationship equation between the estimated parameter P and the unwrapping phase of the dynamic settlement model; $\mathrm{<span\; class="notranslate">\; \&delta;\&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 4</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}\mathrm{<span\; class="notranslate">\; cos</span>}\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&Delta;t</span>}}_{\mathrm{<span\; class="notranslate">\; B</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; A</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 4</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}}{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}\frac{{\mathrm{<span\; class="notranslate">\; B</span>}}_{<span\; class="notranslate">\; \&perp;</span>}\mathrm{<span\; class="notranslate">\; \&Delta;h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; r</span>}\mathrm{<span\; class="notranslate">\; sin</span>}\mathrm{<span\; class="notranslate">\; \&theta;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; )</span>}$ Among them, Δt _B =t _B -t ₀ and Δt _A =t _A -t ₀ are the corrected SAR image time; Step 6: Use the unwrapped phases of time-domain discrete InSAR interferometric pairs greater than or equal to numP+1 and the corresponding radar wavelength λ, incident angle θ, slant distance r, and the vertical baseline length B _⊥ of each interferometric pair to be substituted into step 5 respectively In the formula, the multiple equations obtained by the simultaneous equations are used to calculate the estimated parameter P of the dynamic settlement model and the elevation residual Δh of the high coherence point, and the unwrapped phase of the time-domain discrete InSAR interferometric pair is obtained from step 1 Obtained from the obtained unwrapped phase map; Substitute the estimated parameter P of the dynamic subsidence model into the dynamic subsidence model of the mining area selected in step 4 to calculate the time-series subsidence of the surface at any time, and realize the time-series deformation monitoring of the mining area based on the time-domain discrete InSAR interferometric pair . 2. the mining area surface time-series deformation monitoring method based on time domain discrete InSAR interference pair according to claim 1, is characterized in that, described step 4 selects mining area dynamic model to comprise Knothe model, Gompertz model, Logistic model, Richards model or Any of the Weibull models. 3. the mining area surface time series deformation monitoring method based on time-domain discrete InSAR interference pair according to claim 2, it is characterized in that, when solving the model to be evaluated parameter P in the described step 6, utilize genetic algorithm to search for the pending value of the dynamic settlement model The global optimal solution of estimated parameters is substituted into the dynamic subsidence model of the mining area selected in step 4. 4. The time-series surface deformation monitoring method in mining areas based on time-domain discrete InSAR interferometric pairs according to any one of claims 1-3, wherein the coherence threshold in step 2 is 0.2-0.4.