CN107389029A - A kind of surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology - Google Patents

Abstract

The invention discloses an integrated land subsidence monitoring method based on the fusion of multi-source monitoring technologies, comprising: laying out leveling points and GPS monitoring points in the key monitoring area of land subsidence, and laying out CR-GPS-leveling integrated points in the stable area of the ground surface; and Under the same time series, obtain GPS data, surface meteorological data, MODIS data, SAR images and leveling data; combine the synchronous observation data of surrounding CGPS stations and IGS stations and surface meteorological data to solve GPS data; then combine GPS data and surface meteorological data Solve the atmospheric delay phase information with MODIS data; extract the stable PS point and PS point deformation information from the initial differential interferometric phase map; then use the deformation information corresponding to the PS point, GPS point, benchmark point, and integrated point to construct the vertical deformation of land subsidence Field, the horizontal deformation field and the vertical deformation field are combined to construct a spatial data field, and the three-dimensional deformation field information of land subsidence with high temporal and spatial resolution can be obtained. The invention can obtain large-scale, high-precision, high-time-space resolution three-dimensional deformation information of the earth surface.

Description

An integrated monitoring method of land subsidence based on fusion of multi-source monitoring technology

technical field

The invention relates to the field of land subsidence monitoring, in particular to an integrated monitoring method of land subsidence based on fusion of multi-source monitoring technologies.

Background technique

At present, the monitoring methods of land subsidence mainly include precision leveling survey, bedrock marker-layer marker survey, GPS survey and synthetic aperture radar differential interferometry (InSAR).

Among them, the precision leveling survey obtains the surface deformation information by laying out graded leveling network, adjustment calculation and spatial interpolation. The land subsidence information obtained by this method has high accuracy and reliability, but due to its long re-measurement cycle, manpower The huge consumption of material resources, the inability to meet the requirements of real-time dynamic monitoring of land subsidence, and the discontinuous monitoring information obtained limit the wide application of this method. However, as far as the existing land subsidence monitoring technology is concerned, precision leveling is still unmatched by other monitoring technologies due to its high precision, and it is often used to verify the accuracy of new land subsidence monitoring technology.

The bedrock mark-layer mark monitoring method can obtain the information of vertical layered ground subsidence and deformation with high precision, and its precision can reach 0.01-0.1mm. However, due to the complicated operation, high construction technology and high cost, the wide application of this method in regional land subsidence monitoring is limited, and it is often used in the research of land subsidence mechanism.

GPS measurement technology has played an important role in land subsidence monitoring with the continuous improvement of instruments and unwrapping algorithms. GPS measurement has the advantages of short period, high positioning accuracy, rapid network deployment, and all-weather weather. However, due to the limitation of the winding algorithm, the accuracy of vertical deformation monitoring is still an unavoidable defect. Moreover, what GPS measurement acquires is point-like distribution of ground monitoring point deformation information. In areas with poor signal or obstacles, it is difficult to obtain the elevation value of monitoring points, which limits the use of this method.

Synthetic aperture radar differential interferometry technology is a new type of space earth observation technology developed in the past two decades. It is characterized by real-time fast, large-scale, high-precision, and its vertical deformation monitoring accuracy can reach mm level. However, in terms of horizontal deformation monitoring, its detection ability is limited and it is not sensitive to horizontal deformation. In addition, the phase unwrapping is seriously affected by atmospheric delay and space-time loss of correlation, so it is necessary to eliminate the influence of these errors when solving.

It can be seen that through the analysis of the characteristics of the above-mentioned land subsidence monitoring methods, it is found that the current land subsidence monitoring technology has its own advantages and disadvantages. How to create a new integrated monitoring method of land subsidence based on the fusion of multi-source monitoring technologies, so that it can organically integrate existing subsidence monitoring technologies, perform data fusion on subsidence information obtained by various monitoring methods, and break through a single monitoring technology It is one of the important research and development topics in the field of land subsidence monitoring technology research at present to make use of the respective monitoring advantages of various monitoring methods to obtain large-scale, high-precision, and high-temporal-resolution three-dimensional deformation information of the surface.

Contents of the invention

The technical problem to be solved by the present invention is to provide an integrated land subsidence monitoring method based on the fusion of multi-source monitoring technologies, so that the existing subsidence monitoring technologies can be organically integrated, and the subsidence information obtained by various monitoring means can be data fused. Obtain large-scale, high-precision, and high-spatial-resolution three-dimensional deformation information of the surface, thereby overcoming the shortcomings of existing land subsidence monitoring methods.

In order to solve the above technical problems, the present invention provides an integrated monitoring method for land subsidence based on the fusion of multi-source monitoring technologies, which includes the following steps:

(1) Deploy benchmarking points and GPS monitoring points for land subsidence monitoring in key monitoring areas of land subsidence, and deploy CR-GPS-leveling integrated points in areas with stable ground surfaces;

(2) Under the same time series, obtain GPS data, surface meteorological data, MODIS data, SAR images and leveling data;

(3) Combining the synchronous observation data of surrounding CGPS stations and IGS stations and the above-mentioned surface meteorological data, using the open source software GAMIT to jointly solve the GPS baseline vector, and then using the network adjustment software to perform network adjustment calculation on the GPS data to obtain high Three-dimensional coordinate information of precision ground subsidence GPS monitoring points and CR-GPS-leveling integration points, the three-dimensional coordinate information includes plane position and elevation value;

(4) Combine the GPS data, surface meteorological data and MODIS data to solve the atmospheric delay phase information;

(5) Using DORS software or GAMMA software to carry out differential interferometric processing on the SAR image to obtain an initial differential interferometric phase map;

(6) Using the amplitude dispersion index and spatial phase correlation feature to extract the stable PS points in the initial differential interferogram, estimate the linear deformation and DEM error on each PS point, from the initial differential interferogram Subtract the linear deformation and DEM error on each PS point to obtain the PS-InSAR residual phase;

The PS-InSAR residual phase includes nonlinear deformation phase, atmospheric delay phase and noise. The PS-InSAR residual phase is solved by using a three-dimensional unwrapping algorithm, and the nonlinear deformation phase is separated by high-pass and low-pass filtering techniques. bit and atmospheric delay phase;

(7) The atmospheric delay phase separated from the PS-InSAR residual phase of the step (6) and the atmospheric delay phase obtained by the joint inversion of GPS/MODIS data in the step (4) are subjected to mean value fusion processing to establish a high Atmospheric delay mean model with high precision and high spatio-temporal resolution;

(8) subtracting the atmospheric delay mean phase part after step (7) fusion from the initial differential interferogram, and then obtaining a high-precision PS-InSAR differential interferogram;

(9) With the CR-GPS-leveling integrated point as a reference, phase unwrapping is carried out to the PS-InSAR differential interferogram, and stable PS point deformation information is extracted;

(10) geocoding is carried out to the PS point deformation information extracted in the step (9), unified in the geodetic coordinate reference frame;

(11) Utilize the deformation information corresponding to the PS points, GPS points, benchmarking points, and CR-GPS-leveling integrated points, and use kriging spatial interpolation technology to perform interpolation calculations in the network to construct high-precision, high-spatial-resolution The vertical deformation field of land subsidence realizes the fusion of GPS, InSAR and leveling data in the vertical deformation field;

(12) Interpolate the horizontal monitoring results of the GPS network in the space domain to form the horizontal deformation field of land subsidence, and use the ensemble Kalman filter algorithm to fuse the horizontal deformation field with the vertical deformation field, estimate and predict each point in the network, and construct Unified spatial data field, and then obtain three-dimensional deformation field of land subsidence with high spatial resolution;

(13) Utilizing the characteristics of high time resolution of GPS, based on the GIS platform, interpolation calculation is performed on the three-dimensional deformation field of land subsidence in the time domain, so as to realize the continuous encryption of the three-dimensional deformation field of land subsidence in the time domain, and then obtain 3D deformation field information of land subsidence with high spatio-temporal resolution.

As an improvement of the present invention, the step (4) comprises the step of utilizing GPS data to correct the MODIS data, specifically:

A. Use the GPS calculation results obtained in step (3) to combine with surrounding CGPS stations, use GAMIT software to calculate the high-precision tropospheric zenith total delay ZTD within the SAR satellite transit time, and then use the ground meteorological data to calculate the zenith static The mechanical delay ZHD, and then the zenith wet delay ZWD is obtained by subtracting the zenith static delay ZHD from the total tropospheric zenith delay ZTD, and the calculation formula is as follows:

In the formula, P _s is the surface atmospheric pressure value, is the latitude of the GPS site, H is the elevation value of the GPS site;

B. Convert the atmospheric precipitable water vapor content PWV retrieved by MODIS into the zenith wet delay ZWD', the relationship between the two is:

Among them, ρ _w is the density of liquid water, T _M is the average atmospheric temperature, R ₀ is the universal gas constant, M _W is the molar mass of liquid water, k2 and k3 are atmospheric refraction constants, and the value range of Π is 6.0-6.5;

C. Regressively fit the zenith wet delay ZWD calculated in step A and the zenith wet delay ZWD' obtained in step B to realize the correction of GPS data to MODIS data.

As a further improvement, the step (4) also includes the step of eliminating cloud pollution in the MODIS data, specifically: using the cloud product of MODIS as a mask, removing the pixels that have clouds in the MODIS inversion water vapor, and using spatial interpolation at the same time The technology interpolates the water vapor content value of the cloud-polluted area from the surrounding pixels.

As a further improvement, the space interpolation technique is reciprocal distance weighting method, spline function interpolation method or Kriging interpolation method.

As a further improvement, the step (4) includes a delay phase The calculation formula is: where θ _inc is the radar incidence angle, and ZWD is the corrected zenith wet delay;

Then calculate the atmospheric delay phase of the SAR main and auxiliary image acquisition time The calculation formula is:

Further improvement, the delayed phase Perform low-pass filtering.

As a further improvement, after the step (4), it also includes the step of correcting the satellite orbit error using CR-GPS-leveling integration point GPS measurement data, wherein the CR-GPS-leveling integration point GPS measurement data is used to accurately obtain the CR-GPS-leveling integration point GPS measurement data. The exact position of the GPS-leveling point in the SAR image, and the reverse operation of the geocoding to obtain the precise orbit information of the SAR satellite;

Or the step of correcting the satellite orbit error by using the SAR satellite precise orbit ephemeris file.

After adopting such design, the present invention has the following advantages at least:

1. The present invention can unify the GPS measurement data and InSAR data into the same reference coordinate system by setting the CR-GPS-leveling integrated point, correct the SAR satellite orbit error with the GPS measurement data, and use the GPS data on the integrated point to invert the atmosphere At the same time, the integration of the integrated point with GPS points and benchmarking points can achieve the purpose of encrypted monitoring sites.

2. The present invention proposes to comprehensively utilize the synchronous observation data of surrounding CGPS stations and IGS stations and ground meteorological data for joint calculation, and obtain high-precision GPS calculation results within a unified reference frame. This step is obviously different from the existing GPS calculation. The height calculation accuracy of the invention can reach the second-class leveling accuracy, and fully meet the requirements of land subsidence monitoring in terms of vertical monitoring accuracy.

3. The present invention calculates atmospheric delay phase information by using external data (high-precision GPS data, surface meteorological data and MODIS data), and at the same time combines the advantages of PS-InSAR technology to eliminate part of atmospheric delay, and jointly inverts GPS/MODIS data Atmospheric delay phase separated from the PS-InSAR residual phase and the atmospheric delay phase separated from the PS-InSAR residual phase are subjected to mean fusion processing, and an atmospheric delay mean model with high precision and high temporal and spatial resolution is established. After the coordinate system is unified, the fused atmospheric delay average phase is subtracted from the initial differential interferometric phase to eliminate the influence of atmospheric delay, thereby obtaining high-precision differential interferometric phase information.

4. The present invention fully considers that the ground deformation results acquired by different monitoring techniques have different spatial dimensions and directions. Therefore, the line-of-sight deformation value on the PS point is projected to the vertical direction to achieve unity with the vertical deformation results of GPS and leveling. Then construct a Delaunay triangle network with high-density PS points, GPS points, benchmarking points and integrated points, and evaluate the accuracy and stability of each point in the network. At the same time, kriging spatial interpolation is carried out in the network to build a high-precision, high-spatial-resolution vertical deformation monitoring network for land subsidence. The horizontal deformation field of land subsidence is obtained by spatial interpolation using the horizontal deformation values obtained from GPS points and integrated points. Then, the ensemble Kalman filter algorithm is used to fuse the horizontal deformation field and the vertical deformation field to construct a unified spatial data field and obtain a three-dimensional deformation field of land subsidence with high spatial resolution. Finally, the GPS real-time monitoring results are interpolated in the time domain in the three-dimensional deformation field, and then the three-dimensional deformation field of land subsidence with high temporal and spatial resolution is obtained, which realizes the effective integration of multi-source monitoring technology and overcomes the limitation of a single monitoring method.

Description of drawings

The above is only an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific embodiments.

Fig. 1 is the overall flowchart of the land subsidence integrated monitoring method of the present invention;

Fig. 2 is the step flow chart of the joint correction of atmospheric delay of GPS/ground meteorological data/MODIS data of land subsidence integrated monitoring method of the present invention;

Fig. 3 is the flow chart of the steps of SAR interferometry correction by the atmospheric delay mean model of the land subsidence integrated monitoring method of the present invention;

Fig. 4 is a schematic diagram of InSAR interferometry imaging in the present invention;

Fig. 5 is a coordinate conversion diagram between GPS and SAR in the present invention.

detailed description

The specific steps of the invention of the land subsidence integrated monitoring method are described in detail in conjunction with the accompanying drawings.

With reference to shown in accompanying drawing 1, the land subsidence integrated monitoring method of the present invention comprises the steps:

1) Lay out benchmarking points, GPS points and artificial corner reflector (CR)-GPS-leveling integrated points

Firstly, precision leveling points and GPS monitoring points are arranged in the key monitoring areas of land subsidence for ground subsidence monitoring; CR-GPS-leveling integrated points are arranged in areas with relatively stable ground surface, which can be used as reference points for subsidence monitoring and air-ground monitoring points. Integral connection reference point. And the CR-GPS-leveling integrated point can unify the subsequent GPS measurement data and InSAR data into the same reference coordinate system, and the GPS data measured at this point can be used to correct SAR satellite orbit errors and invert atmospheric delays, eliminating InSAR phase Atmospheric effects in unwrapping.

2) Obtain GPS data, surface meteorological data, MODIS data, SAR images, leveling data

Specifically, the collected GPS data, surface meteorological data, MODIS data, and SAR images should have the same time series, while the time series of leveling can be appropriately relaxed due to the slowness of its measurement.

3) Solving GPS data

Combining the synchronous observation data of surrounding CGPS stations and IGS stations and surface meteorological data, the open source software GAMIT is used to jointly solve the GPS baseline vector, and the network adjustment software is used to perform network adjustment calculation on the GPS measurement data to obtain high-precision ground subsidence GPS monitoring points and CR-GPS- leveling point three-dimensional coordinate information, including plane position and elevation value.

4) Calculate the atmospheric delay phase value by combining GPS data, surface meteorological data and MODIS data

Referring to Figure 2, ① using GPS monitoring points and CR-GPS-leveling integration points to obtain GPS solution results combined with surrounding CGPS stations, using GAMIT software to calculate the high-precision tropospheric zenith total delay during the transit time of SAR satellites ( ZTD). The zenith static delay (ZHD) is calculated using the surface meteorological data, and then the zenith wet delay (ZWD) is obtained through ZTD-ZHD. Calculated as follows:

In the formula, P _s is the surface atmospheric pressure value, is the latitude of the GPS site, and H is the elevation value of the GPS site.

② Since the MODIS inversion obtains the atmospheric precipitable water vapor content (PWV), when using GPS data to correct the MODIS water vapor value, the PWV obtained by MODIS inversion should be converted into ZWD. The relationship between the two is:

Among them, ρ _w is the density of liquid water, T _M is the average atmospheric temperature, R ₀ is the universal gas constant, M _W is the molar mass of liquid water, k2 and k3 are atmospheric refraction constants, and the value range of Π is 6.0-6.5.

The zenith wet delay GPS (ZWD) calculated by GPS combined with surface meteorological data and the zenith wet delay MODIS (ZWD) obtained by inversion of MODIS data are used for regression fitting, and then the correction of MODIS (ZWD) by GPS (ZWD) is achieved. Purpose.

③When there are clouds in the atmosphere, the MODIS data will not correctly reflect the water vapor content in the atmosphere. Therefore, it is necessary to use the MODIS cloud products as a mask to remove the pixels with clouds in the MODIS inversion water vapor. At the same time, the spatial interpolation technology (IDW, Spline interpolation or Kriging interpolation) is used to interpolate the water vapor content of the cloud-polluted area from the surrounding pixels.

④Since the SAR differential interferogram contains phase information, in order to remove the atmospheric delay calculated jointly by GPS and MODIS from the interferogram, it is necessary to convert the path delay into a delay phase delay phase The calculation formula is:

In the formula, θ _inc is the radar incident angle, and λ is the wavelength. In order to weaken the influence of noise and operating errors, it is necessary to Perform low-pass filtering.

⑤ Calculate the differential atmospheric delay phase at the acquisition time of the SAR main and auxiliary images. The formula is:

5) Use CR-GPS-leveling point GPS measurement data or SAR satellite precise orbit ephemeris files to correct satellite orbit errors

Referring to Figure 3, due to the need for rough registration and precise registration of the main and auxiliary images in the process of repeated orbit SAR interferometry, and satellites will have orbit deviations during their revisit periods, eliminating satellite orbit errors is essential for SAR Successful differential interferometry of images is of great significance. Therefore, for SAR data that cannot obtain precise orbit files, the present invention uses the established CR-GPS-leveling integrated point GPS measurement data to accurately obtain the exact position of the integrated point in the SAR image, and uses the inverse operation of geocoding to reverse the precise orbit of the SAR satellite information. For SAR data that can obtain precise orbit ephemeris files, the precise orbit files can be directly used to correct satellite orbit errors.

6) SAR image differential interference processing

Use DORS software or GAMMA software to perform differential interferometric processing on time series SAR images to obtain differential interferometric phase maps. After differential interference processing, each pixel contains the following components:

In the formula: is the interferometric phase of the point target; Deformation phase for the radar line of sight direction; is the terrain phase; is the atmospheric delay phase; is the orbit error phase; is the noise phase;

Among them, the principle of differential interference:

Synthetic Aperture Radar Interferometry (InSAR) obtains terrain or deformation phase by performing differential interferometry on twice acquired SAR data. Because it is difficult to set up dual antennas on a single satellite, spaceborne SAR generally uses repeated orbits for interferometry. Synthetic aperture radar interferometry is by measuring the echo phase information of the target object, using the different spatial position relationship of the two radar imaging, and according to the triangular similarity principle, to obtain the deformation information of the ground object - the elevation or motion state of the target object ( speed, attitude, etc.). The interferometric radar measurement system transmits radar signals to the ground through a single antenna, and then uses dual antennas to simultaneously receive the reflected echoes of ground objects. Since the dual antennas have a time difference when receiving echo signals, interferometry results in different time periods can be obtained. Therefore, the elevation information of surface objects can be obtained through synthetic aperture radar differential interferometry. Repeated orbit radar interferometry has two main purposes, one is to measure surface elevation information, and the other is to monitor surface deformation information. Since the satellite orbits at different times do not completely coincide during repeated orbit interferometry, the phase signal obtained by interferometry contains both terrain phase information and line-of-sight direction displacement information.

The principle of InSAR differential interferometry is shown in Figure 4. In the figure, A ₁ and A ₂ represent the positions of the dual antennas respectively, R ₁ and R ₂ are the paths from the two ends of the antenna to a certain target on the ground surface, θ ₁ and θ ₂ is the incident angle, the baseline B is the spatial distance of the antenna between the two acquisitions of ground SAR images, B _|| is the parallel component of the baseline, and B _⊥ is the vertical component of the baseline. The angle between the baseline B and the horizontal direction is α, H is the height of the sensor, and Z is the elevation value of the surface terrain. Wherein, the SAR signals received by the antennas A1 and A2 are represented by the following equations ( ₁ ) and ( ₂ ):

Since the spatial position of the radar satellite is not the same during the two transits, the SAR images of the same area acquired twice do not completely overlap, and precise orbit files and main images need to be used for registration. Complex conjugate multiplication is performed on the two registered SAR images to generate an interferogram. The result of the interference is as follows (3):

The phase caused by the deformation Δ _r in the line of sight direction (Δ _r = R ₁ -R ₂ , which is the path length difference) can be calculated as follows (4) and (5):

or

Among them, φ _d is the deformation phase; λ is the wavelength; Δ _r is the difference between the two path lengths; R ₁ is the path length of the first transit of the satellite; R ₂ is the path length of the second transit of the satellite. The deformation amount Δ _r of the target point can be calculated by formula (4).

7) Extract PS point phase information and remove error components

Using amplitude dispersion index and spatial phase correlation features to extract stable PS points, after estimating the linear deformation and DEM error on each PS point, subtract them from the initial differential interferometric phase to obtain the residual phase, which is mainly Includes nonlinear deformation phase, atmospheric phase, and noise. A three-dimensional unwrapping algorithm is used to solve the residual phase, and the nonlinear deformation and atmospheric phase are separated by high-pass and low-pass filtering techniques.

8) Atmospheric delay phase separated from PS-InSAR residual phase Atmospheric delay phase obtained by joint inversion with GPS/MODIS data Perform mean fusion processing to establish an atmospheric delay mean model with high precision and high spatio-temporal resolution.

Since the atmospheric delay phase separated from the PS-InSAR residual phase is inconsistent with the coordinate system of the atmospheric delay phase obtained by joint inversion of GPS/MODIS data, it is necessary to use the coordinate conversion formula to unify the coordinate system and unify the two to the radar coordinate system. Then use the raster calculation tool to calculate the mean value of the two at the pixel scale to obtain The formula is:

The coordinate conversion formula is as follows:

Surface deformation rate obtained by GPS monitoring is the three-dimensional deformation information, based on the three-dimensional unit vector It can be decomposed into deformation values of three vector directions (due east, true north, vertical direction). As shown in Figure 5.

When the SAR satellite is deorbited and passed through the scan θ and α are the incident angle and azimuth angle of the SAR satellite, respectively.

In SAR images, the ground deformation rate can be decomposed into two-dimensional unit vectors Where i∈{descending, ascending}; that is, i∈{descending orbit, ascending orbit}:

in: with Represent the line-of-sight and azimuth deformations, respectively.

The deformation values in three directions acquired by GPS can be projected onto the SAR satellite geometric space through projection transformation:

(1) SAR descending orbit scanning:

Convert GPS satellite B _3d coordinate system to The corresponding coordinate system is as follows:

with

According to the above formula, the conversion relationship between GPS and SAR images is obtained:

(2) SAR ascending orbit scanning:

Convert GPS satellite B _3d coordinate system to with The corresponding coordinate system is as follows:

with

According to the above formula, the conversion relationship between GPS and SAR images is obtained:

9) Subtract the fused atmospheric delay mean phase part from the initial differential interferophase.

Since the removal of the atmospheric delay phase is carried out for the pixel, it is also necessary to The phase map is registered with the initial differential interferogram to achieve the unity of the pixel scale, and then the grid calculation can be performed to obtain a high-precision differential interferogram that eliminates the influence of atmospheric delay.

In the formula, is the differential interferogram after atmospheric delay correction, is the initial differential interferogram.

10) Taking the CR-GPS-leveling integrated point as a reference, the phase unwrapping is performed on the atmospherically corrected PS-InSAR differential interferogram, and the stable PS point deformation information is extracted.

11) Geocode the PS point deformation information extracted by PS-InSAR, and unify it into the geodetic coordinate reference frame.

12) Use high-density PS points, GPS points, benchmarking points, and CR-GPS-leveling integrated points to fuse vertical deformation monitoring results.

Since the deformation rate of the PS point obtained by single-track SAR differential interferometry is the deformation value of the radar line of sight, the deformation rate obtained by GPS is the three-dimensional deformation information, which can be decomposed into due east, due north and vertical directions. Therefore, when merging the line-of-sight deformation value of PS point with the vertical deformation value of GPS, it is necessary to project the line-of-sight deformation value of PS point to the vertical direction. The calculation formula is: U _u =d _los /cosθ where d _los is the deformation value in the direction of the radar line of sight, θ is the incident angle of the radar wave, and U _u is the deformation value in the vertical direction of the radar. Then, the high-density PS points, GPS points, benchmarking points and integrated points are used to construct a Delaunay triangle network using the GIS spatial analysis method, and the accuracy and stability of each point in the network are evaluated. Utilizing the deformation information corresponding to the above monitoring points, kriging spatial interpolation technology is used to perform interpolation calculations in the network to construct a high-precision, high-spatial-resolution vertical deformation monitoring network for land subsidence, and to realize GPS, InSAR and leveling data in the network. Fusion of vertical deformation fields.

13) The horizontal monitoring results of the GPS network are interpolated in the spatial domain to form the horizontal deformation field of land subsidence, and at the same time, the ensemble Kalman filter algorithm is used to integrate the horizontal deformation field with the above-mentioned vertical deformation field, estimate and predict each point in the network, and build a unified The spatial data field is used to obtain the three-dimensional deformation field of land subsidence with high spatial resolution.

14) Utilizing the characteristics of high time resolution of GPS, based on the GIS platform, the above-mentioned three-dimensional deformation field of land subsidence is interpolated in the time domain, so as to realize the continuous encryption of the three-dimensional deformation field of land subsidence in the time domain, and then obtain a high-temporal-spatial High-resolution 3D deformation field information of land subsidence.

The present invention utilizes land subsidence data obtained by various monitoring methods and the characteristics of target ground objects. Based on the principle of information optimization, more continuous, comprehensive and comprehensive monitoring information can be obtained through organic combination, which can not only improve the accuracy of land subsidence information extraction It also reduces the ambiguity and uncertainty of the overall understanding of land subsidence phenomena.

The present invention comprehensively utilizes multi-source ground subsidence monitoring means, integrates multi-source ground subsidence monitoring technologies, performs data fusion on subsidence information acquired by various monitoring means, overcomes the defects of existing monitoring technologies, and breaks through the limitations of a single monitoring technology. The land subsidence comprehensive monitoring method based on the integration of multi-source monitoring technology can obtain large-scale, high-precision, and high-spatial-resolution three-dimensional deformation information of the ground surface.

The above is only a preferred embodiment of the present invention, and does not limit the present invention in any form. Those skilled in the art make some simple modifications, equivalent changes or modifications by using the technical content disclosed above, all of which fall within the scope of the present invention. within the scope of protection of the invention.

Claims (7)

1. A kind of 1. surface subsidence integrated monitor method based on the fusion of multi-source monitoring technology, it is characterised in that comprise the following steps：

   (1) bench mark and the GPS monitoring points for Ground Subsidence Monitoring are laid in surface subsidence emphasis monitored area, it is steady in earth's surface Lay CR-GPS- levels one point in fixed region；

   (2) under identical time series, gps data, ground meteorological data, MODIS data, SAR images and level is obtained and is surveyed Measure data；

   (3) combine periphery CGPS stations and IGS station simultaneous observation data and the ground meteorological data, utilize open source software GAMIT Combined Calculation GPS basic lineal vectors, then balancing calculation of GPS net is carried out to the gps data using net adjusted data software, obtain high-precision The three-dimensional coordinate information that degree surface subsidence GPS monitoring points and CR-GPS- levels are integrally put, the three-dimensional coordinate information include plane Position and height value；

   (4) gps data, ground meteorological data and the MODIS data calculation atmosphere delay phase informations are combined；

   (5) differential interferometry processing is carried out to the SAR images using DORS softwares or GAMMA softwares, obtains initial differential interference Phase diagram；

   (6) extracted using amplitude dispersion index and space phase correlative character stable in the initial differential interferometric phase image PS points, estimate the linear deformation on each PS points and DEM errors, will be described each from the initial differential interferometric phase image Linear deformation and DEM errors on PS points subtract, and produce PS-InSAR residual phases；

   The PS-InSAR residual phases include non-linear deformation phase, atmosphere delay phase and noise, to the PS-InSAR Residual phase twines algorithm using Three-Dimensional Solution and resolved, and non-linear deformation phase and big is isolated using high and low pass filtering technique Gas postpones phase；

   (7) by the atmosphere delay phase separated in the PS-InSAR residual phases of the step (6) and the step (4) The atmosphere delay phase that GPS/MODIS data aggregate invertings obtain does average fusion treatment, establishes high accuracy, high-spatial and temporal resolution Atmosphere delay mean value model；

   (8) the atmosphere delay average phase bit position after step (7) fusion is subtracted from the initial differential interferometric phase image, and then Obtain high-precision PS-InSAR differential interferometries phase diagram；

   (9) with the CR-GPS- levels, integrally for reference data, phase is carried out to the PS-InSAR differential interferometries phase diagram for point Solution twines, and extracts stable PS point deformation datas；

   (10) geocoding is carried out to the PS points deformation data of extraction in the step (9), it is unified to arrive geodetic coordinates reference frame It is interior；

   (11) using deformation data corresponding to the PS points, GPS point, bench mark, CR-GPS- levels one point, using Ke Lijin Spatial interpolation technology carries out interpolation calculation in net, and structure high accuracy, the surface subsidence vertical deformation field of high spatial resolution are real Show the fusion of GPS, InSAR and measurement of the level data in vertical deformation field；

   (12) GPS network level monitoring result is subjected to spatial domain interpolation and forms surface subsidence horizontal deformation field, while utilize set Kalman filtering algorithm is merged horizontal deformation field with the vertical deformation field, and estimation prediction, structure are carried out to each point in net Unified spatial data field is built, and then obtains high spatial resolution surface subsidence three-dimensional shaped variable field；

   (13) there is the characteristic of high time resolution using GPS, it is three-dimensional to the surface subsidence in time-domain based on GIS platform Deformation Field carries out interpolation calculating, so as to realize continuous encryption of the surface subsidence three-dimensional shaped variable field in time-domain, and then is had There is the surface subsidence three-dimensional shaped variable field information of high-spatial and temporal resolution.
2. 2. the surface subsidence integrated monitor method according to claim 1 based on the fusion of multi-source monitoring technology, its feature exist In, in the step (4) using gps data to MODIS Data corrections the step of, be specially：

   A, the GPS calculation results joint periphery CGPS stations obtained using step (3), SAR satellite mistakes are calculated using GAMIT softwares High-precision tropospheric zenith total delay ZTD in the time of border, recycles the ground meteorological data to calculate zenith hydrostatic and prolongs Slow ZHD, and then Zenith hydrostatic delay ZHD is subtracted by tropospheric zenith total delay ZTD and draws Zenith wet delay ZWD, calculate Formula is as follows：

   In formula, P_sFor surface air pressure value,For GPS website latitudes, H is GPS website height values；

   B, the precipitable water vapour content PWV that MODIS invertings obtain is changed into Zenith wet delay ZWD ', both sides relation is：

   <mrow> <mo>&amp;Pi;</mo> <mo>=</mo> <mfrac> <mrow> <msup> <mi>ZWD</mi> <mo>&amp;prime;</mo> </msup> </mrow> <mrow> <mi>P</mi> <mi>W</mi> <mi>V</mi> </mrow> </mfrac> <mo>=</mo> <msup> <mn>10</mn> <mrow> <mo>-</mo> <mn>6</mn> </mrow> </msup> <msub> <mi>&amp;rho;</mi> <mi>w</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mn>2</mn> </msub> <mo>+</mo> <mfrac> <msub> <mi>k</mi> <mn>3</mn> </msub> <msub> <mi>T</mi> <mi>M</mi> </msub> </mfrac> <mo>)</mo> </mrow> <mfrac> <msub> <mi>R</mi> <mn>0</mn> </msub> <msub> <mi>M</mi> <mi>W</mi> </msub> </mfrac> </mrow>

   Wherein, ρ_wFor liquid water density, T_MFor Zenith Distance temperature, R₀For universal gas constant, M_WFor aqueous water molal weight, K2, k3 are atmospheric refraction constant, and П spans are 6.0~6.5；

   C, the Zenith wet delay ZWD ' regression fits that will be obtained in the Zenith wet delay ZWD calculated in step A and step B, it is real Existing correction of the gps data to MODIS data.
3. 3. the surface subsidence integrated monitor method according to claim 2 based on the fusion of multi-source monitoring technology, its feature exist In the step of also including eliminating the pollution of MODIS data medium cloud in the step (4), specially：Using MODIS cloud product as Mask, there will be pixel existing for cloud to remove in MODIS inverting steam, while using spatial interpolation technology by by cloud Polluted area Moisture content value is generated by the picture element interpolation of surrounding.
4. 4. the surface subsidence integrated monitor method according to claim 3 based on the fusion of multi-source monitoring technology, its feature exist In the spatial interpolation technology is distance-reverse weighting function, spline interpolation method or Kriging interpolation methods.
5. 5. the surface subsidence integrated monitor method according to claim 2 based on the fusion of multi-source monitoring technology, its feature exist In the step (4) includes postponing phaseCalculating, calculation formula is：In formula, θ_inc For radar incidence angle, ZWD is the Zenith wet delay after correction；

   And then calculate the atmosphere delay phase at SAR main and auxiliary image capturing momentCalculation formula is:
6. 6. the surface subsidence integrated monitor method according to claim 5 based on the fusion of multi-source monitoring technology, its feature exist In the delay phaseCarry out low-pass filtering treatment.
7. 7. the surface subsidence integrated monitor method according to claim 1 based on the fusion of multi-source monitoring technology, its feature exist In, the step of also including integrally putting gps measurement data correction satellite orbital error using CR-GPS- levels after the step (4), Wherein gps measurement data is integrally put using the CR-GPS- levels of foundation accurately obtain the CR-GPS- levels one point in SAR Accurate location in image, and utilize geocoding inverse operation reverse SAR satellite precise orbit information；

   Or the step of using SAR satellite precise orbit ephemeris file correction satellite orbital error.