WO2022214114A2 - 融合GNSS数据与InSAR技术的桥梁变形监测方法 - Google Patents

融合GNSS数据与InSAR技术的桥梁变形监测方法 Download PDF

Info

Publication number
WO2022214114A2
WO2022214114A2 PCT/CN2022/106155 CN2022106155W WO2022214114A2 WO 2022214114 A2 WO2022214114 A2 WO 2022214114A2 CN 2022106155 W CN2022106155 W CN 2022106155W WO 2022214114 A2 WO2022214114 A2 WO 2022214114A2
Authority
WO
WIPO (PCT)
Prior art keywords
image
phase
sar
deformation
unwrapped
Prior art date
Application number
PCT/CN2022/106155
Other languages
English (en)
French (fr)
Other versions
WO2022214114A3 (zh
Inventor
崔玉萍
Original Assignee
中咨数据有限公司
中国公路工程咨询集团有限公司
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 中咨数据有限公司, 中国公路工程咨询集团有限公司 filed Critical 中咨数据有限公司
Publication of WO2022214114A2 publication Critical patent/WO2022214114A2/zh
Publication of WO2022214114A3 publication Critical patent/WO2022214114A3/zh
Priority to ZA2023/00146A priority Critical patent/ZA202300146B/en

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B15/00—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons
    • G01B15/06—Measuring arrangements characterised by the use of electromagnetic waves or particle radiation, e.g. by the use of microwaves, X-rays, gamma rays or electrons for measuring the deformation in a solid
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques

Definitions

  • the invention belongs to the technical field of bridge deformation monitoring, in particular to a bridge deformation monitoring method integrating GNSS data and InSAR technology.
  • bridges As an important part of transportation infrastructure, bridges have become an indispensable part of urban economic development. However, the bridge itself is heavy and the load is complex and changeable. Factors such as changes in geological conditions, overloading, and bad weather will have many adverse effects on the bridge, bringing hidden dangers to the safety of the bridge, affecting the safety of traffic, and even threatening the safety of people's lives and property. Therefore, strengthening the deformation monitoring of bridges is of great significance to the safe operation of bridges.
  • Synthetic aperture radar interferometry technology uses the phase information of SAR images to achieve high-precision monitoring of surface deformation, and its deformation monitoring accuracy can reach centimeters or even millimeters.
  • the technical characteristics of range, high density, long time series and traceability have unique advantages compared with traditional GNSS, leveling and other methods.
  • the present invention provides a bridge deformation monitoring method integrating GNSS data and InSAR technology, which can effectively solve the above problems.
  • the present invention provides a bridge deformation monitoring method integrating GNSS data and InSAR technology, comprising the following steps:
  • Step 1 Determine the bridge monitoring object; obtain g high-resolution SAR images covering the bridge monitoring object;
  • g high-resolution SAR images are SAR images obtained at different monitoring times; each high-resolution SAR image has imaging time attributes, image backscatter information and imaging orbit information;
  • Step 2 in each high-resolution SAR image, determine one high-resolution SAR image as the SAR master image, and the other g-1 high-resolution SAR images as the SAR slave image;
  • Step 3 According to the imaging orbit information and image backscattering information of each high-resolution SAR image, rotate, zoom, and sample each SAR slave image to realize the image registration of each SAR slave image and the SAR master image. This obtains g-1 registered SAR slave images;
  • Step 4 For the SAR master image and the g-1 registered SAR slave images, any two SAR images form an interferometric image pair, thereby generating all interferometric image pairs;
  • the two SAR images contained in it are conjugately multiplied to generate an interferometric phase map
  • Step 5 for each interference phase map in the set of interference phase maps, perform preprocessing, mask processing and phase unwrapping processing to obtain an unwrapped differential interference phase map;
  • the set of interferometric phase maps is converted into a set of unwrapped differential interference phase maps
  • Step 6 for the set of unwrapped differential interference phase images obtained in step 5, it is assumed that there are M sets of unwrapped differential interference phase images, respectively expressed as: unwrapped differential interference phase image tu 1 , unwrapped differential interference phase image tu 2 ,. .., the unwrapped differential interferometric phase map tu M ;
  • M interferometric image pairs there are M interferometric image pairs in total, that is, when the interferometric image pairs are not considered to contain common SAR images, the M interferometric image pairs contain 2M SAR images; and then the 2M SAR images are deduplicated processing, assuming that N+1 SAR images are finally obtained;
  • the N+1 SAR images are sorted by imaging time and expressed as: SAR image I 0 , SAR image I 1 , ..., SAR image I N , and the corresponding imaging times are: t 0 , t 1 , . .., t N ; that is: the imaging time of the SAR image I 0 is t 0 , the imaging time of the SAR image I 1 is t 1 , . . , and the imaging time of the SAR image I N is t N ;
  • the unwrapped differential interference phase map tu j (I c , I d ) the The unwrapped phase difference can be directly obtained by reading the unwrapped differential interference phase image;
  • Step 7 the M unwrapped differential interferometric phase images and the N+1 SAR images have the same image size and the same number of pixels;
  • Step 7.1 establish the following equation:
  • Step 7.2 solve the equation established in step 7.1, and obtain the relative phase change sequence of the InSAR radar line-of-sight deformation of the bridge at the pixel (l, p) position in the image coordinate system That is: the relative phase change value of the InSAR radar line-of-sight deformation of the InSAR radar in the image coordinate system at the time t 1 , ..., t N relative to the time t 0 at the position of the pixel (l, p);
  • Step 7.3 according to the SAR imaging geometric parameters and wavelength information, the relative phase change sequence of the InSAR radar line-of-sight deformation of the bridge at the pixel (l, p) position in the image coordinate system Converted to the relative displacement change sequence of the InSAR radar line - of - sight deformation in the image coordinate system
  • the relative displacement change sequence is geocoded to obtain the relative displacement change sequence of the InSAR radar line-of-sight deformation variable in the geographic coordinate system; it is expressed as:
  • Step 8 Arrange num GNSS measuring points on the bridge monitoring object
  • each GNSS measuring point can continuously measure the three-dimensional deformation and displacement value of the location of the measuring point according to the sampling interval;
  • Deformation map gr i the corrected bridge deformation map gr' i is obtained, which represents the displacement change value of the InSAR radar line-of-sight deformation of each pixel in the geographic coordinate system at time t i relative to the initialization time.
  • the purpose of the invention is to overcome the low spatial resolution of deformation monitoring caused by the sparse layout of traditional GNSS measurement and other bridge deformation monitoring means, the InSAR technical means mostly obtain the deformation variables relative to the reference point of unknown deformation situation, and the linear deformation model commonly used by InSAR in bridges. Due to the large error of deformation monitoring, a bridge deformation monitoring method that integrates GNSS data and InSAR technology is proposed. The existing measurement results of GNSS and InSAR high spatial resolution deformation monitoring capability are combined to greatly improve the deformation of bridge infrastructure. Spatial resolution of monitoring.
  • Fig. 1 is the schematic flow chart of the bridge deformation monitoring method of fusion GNSS data and InSAR technology provided by the present invention
  • Fig. 2 is the interferometric image pair combination diagram of the SAR image
  • FIG. 3 is a schematic diagram of a processed differential interference phase diagram
  • Figure 4 is a schematic diagram of the unwrapped differential interference phase diagram with the correct unwrapping result
  • Figure 5 is the InSAR deformation distribution map of the bridge
  • FIG. 6 is a schematic diagram of a method of determining the coefficient matrix.
  • the purpose of the invention is to overcome the low spatial resolution of deformation monitoring caused by the sparse layout of traditional GNSS measurement and other bridge deformation monitoring means, the InSAR technical means mostly obtain the deformation variables relative to the reference point of unknown deformation situation, and the linear deformation model commonly used by InSAR in bridges. Due to the large error of deformation monitoring, a bridge deformation monitoring method that integrates GNSS data and InSAR technology is proposed. The existing measurement results of GNSS and InSAR high spatial resolution deformation monitoring capability are combined to greatly improve the deformation of bridge infrastructure. Spatial resolution of monitoring.
  • the bridge deformation monitoring method of fusion GNSS data and InSAR technology comprises the following steps:
  • Step 1 Determine the bridge monitoring object; obtain g high-resolution SAR (Synthetic Aperture Radar, Synthetic Aperture Radar) images covering the bridge monitoring object.
  • g high-resolution SAR images are SAR images obtained at different monitoring times; each high-resolution SAR image has imaging time attributes, image backscatter information and imaging orbit information.
  • high-resolution SAR images require ground pixels of SAR data to be better than 3 meters, which can be high-resolution radar satellites such as COSMO-SkyMed radar satellite constellation, TerraSAR-X/TanDEM-X radar satellite, and RadarSAT-2 radar satellite.
  • the beamforming mode of the system and the strip mode of high-resolution imaging provide data sources for obtaining high-resolution spaceborne SAR data sets of repeated bridge orbits; at the same time, the vertical baseline between the selected SAR images needs to be as short as possible.
  • -X as an example, the vertical baseline distribution of the SAR data set is required to be within 300 meters to minimize the impact of DSM elevation errors.
  • Step 2 in each high-resolution SAR image, determine one high-resolution SAR image as the SAR master image, and the other g-1 high-resolution SAR images as the SAR slave image.
  • Step 3 According to the imaging orbit information and image backscattering information of each high-resolution SAR image, rotate, zoom, and sample each SAR slave image to realize the image registration of each SAR slave image and the SAR master image. This results in g-1 registered SAR slave images.
  • the registration offset of the master-slave SAR image can be calculated based on the SAR image imaging orbit information, and the master-slave SAR image can be roughly registered; then the intensity-normalized cross-correlation algorithm can be used to accurately register the master-slave SAR image. , the registration accuracy needs to be better than 1/8 pixel.
  • Step 4 for the SAR master image and the g-1 registered SAR slave images, any two SAR images form an interferometric image pair, thereby generating all interferometric image pairs.
  • the two SAR images it contains are conjugately multiplied to generate an interferometric phase map. Therefore, for all interferometric image pairs, multiple interferometric phase maps are generated, thereby forming a set of interferometric phase maps.
  • Step 5 for each interference phase map in the set of interference phase maps, perform preprocessing, mask processing and phase unwrapping processing to obtain an unwrapped differential interference phase map.
  • the set of interferometric phase maps is converted into a set of unwrapped differential interference phase maps.
  • Step 5 is specifically:
  • Step 5.1 preprocessing:
  • Each interferometric phase image is preprocessed, including ground phase removal, terrain phase removal, and filtering, to obtain a processed differential interferometric phase image.
  • high-precision DSM data covering bridge monitoring objects; perform terrain phase removal based on high-precision DSM data.
  • high-precision DSM data requires strong current situation, and the elevation information of bridge monitoring objects provided by it is required after the bridge is built; the accuracy and spatial resolution of high-precision DSM data are higher than AW3D30 DSM data.
  • the Gaofen-7 stereo pair, the WorldView-3 stereo pair, the TerraSAR-X/TanDEM-X double star system, and the UAV oblique photogrammetry can provide high-resolution and high-precision DSM data for the study area.
  • the processed differential interference phase map is schematic diagram.
  • Step 5.2 mask processing:
  • Step 5.2.1 making a bridge monitoring object mask; wherein, the bridge monitoring object mask means: in the mask image, the pixel point of the bridge monitoring object is 1, and the pixel points of all other positions are 0.
  • Step 5.2.2 using the bridge monitoring object mask to perform mask processing on each processed differential interference phase image to obtain a differential interference phase image after mask processing; wherein, the differential interference phase image after mask processing , the phase of the pixel point of the bridge monitoring object remains unchanged, and the phase of the pixel point in other positions becomes 0 uniformly.
  • phase unwrapping algorithms are used to perform phase unwrapping operations, thereby obtaining multiple unwrapped differential interferometric phase images, and a screening algorithm is used to eliminate the unwrapping error.
  • Differential interferometric phase map preserving the unwrapped differential interferometric phase map with the correct phase unwrapping result.
  • the unwrapped differential interferometric phase maps with correct unwrapping results obtained by multiple phase unwrapping algorithms are combined, and finally multiple unwrapped differential interferometric phase maps are obtained.
  • phase unwrapping for each differential interferometric phase image processed by the mask, a reference point is selected in the phase stable area of the bridge, and various phase unwrapping algorithms such as the minimum cost flow method and the branch-cut method are used respectively, and all differential interferometric phase images are processed. Phase unwrapping.
  • Step 6 for the set of unwrapped differential interference phase images obtained in step 5, it is assumed that there are M sets of unwrapped differential interference phase images, respectively expressed as: unwrapped differential interference phase image tu 1 , unwrapped differential interference phase image tu 2 ,. .., the unwrapped differential interferometric phase map tu M .
  • M unwrapped differential interferometric phase maps there are M interferometric image pairs in total, that is, when the interferometric image pairs are not considered to contain common SAR images, the M interferometric image pairs contain 2M SAR images; and then the 2M SAR images are deduplicated Processing, it is assumed that N+1 SAR images are finally obtained.
  • the N+1 SAR images are sorted by imaging time and expressed as: SAR image I 0 , SAR image I 1 , ..., SAR image I N , and the corresponding imaging times are: t 0 , t 1 , . .., t N ; that is, the imaging time of the SAR image I 0 is t 0 , the imaging time of the SAR image I 1 is t 1 , . . . , and the imaging time of the SAR image IN is t N .
  • the unwrapped differential interference phase map tu j (I c , I d ) the The unwrapped phase difference can be obtained directly by reading the unwrapped differential interferometric phase map.
  • Step 7 the M unwrapped differential interferometric phase images and the N+1 SAR images have exactly the same size and have the same number of pixels.
  • Step 7.1 establish the following equation:
  • Step 7.2 solve the equation established in step 7.1, and obtain the relative phase change sequence of the InSAR radar line-of-sight deformation of the bridge at the pixel (l, p) position in the image coordinate system That is: the relative phase change value of the InSAR radar line-of - sight deformation variable in the image coordinate system at the time t 1 , .
  • step 6 For the convenience of understanding step 6 and step 7, an embodiment is listed below:
  • step 3 there are a total of 4 SAR images, namely SAR image I0, SAR image I1, SAR image I2, and SAR image I3.
  • step 4 any two SAR images form an interferometric image pair, thereby generating 6 interferometric image pairs.
  • an interferometric phase diagram is generated, so a total of 6 interferometric phase diagrams are obtained.
  • step 5 by performing preprocessing, mask processing and phase unwrapping processing on each interferometric phase image, the differential interferometric phase images with unwrapping errors are eliminated, and the unwrapped differential interferometric phase images with correct phase unwrapping results are retained. Therefore, in this step, three unwrapped differential interference phase images are finally obtained, which are expressed as:
  • Unwrapped differential interferometric phase map tu 3 the corresponding interferometric image pairs are: SAR image I0 and SAR image I3; the unwrapped phase difference of any pixel is expressed as is a known value.
  • Unwrapped differential interferometric phase map tu 4 the corresponding interferometric image pairs are: SAR image I0 and SAR image I2; the unwrapped phase difference of any pixel is expressed as is a known value.
  • Unwrapped differential interferometric phase map tu 5 the corresponding interferometric image pairs are: SAR image I1 and SAR image I3; the unwrapped phase difference of any pixel is expressed as is a known value.
  • SAR image I0 SAR image I0
  • SAR image I1 SAR image I2
  • the coefficient matrix is thus determined.
  • N+1 SAR images are involved.
  • the position of the pixel on the unwrapped differential interferometric phase image can be described by the azimuth coordinate l and the range coordinate p.
  • the unwrapped phase of the jth unwrapped differential interferometric phase image at the pixel (l, p) can be expressed as:
  • the jth unwrapped differential interferometric phase image is generated from the SAR image with imaging time t B and the SAR image with imaging time t A after conjugation, preprocessing, mask processing and phase unwrapping.
  • d(t B , l, p) and d(t A , l, p) are pixel (l, p) at time t B and t A relative to the start time of the SAR dataset The amount of deformation of the radar line of sight.
  • ⁇ is the radar wavelength
  • r is the distance between the pixel (l, p) and the satellite
  • ⁇ is the local incident angle of the pixel
  • B j ⁇ is the time baseline of the interference pair
  • ⁇ z is the elevation residual
  • ⁇ n j is the noise phase.
  • the phase of the elevation residual caused by the elevation residual can be ignored; at the same time, since the coverage of bridges is often small, it can be considered that within this range
  • the atmospheric distribution is relatively uniform, so the atmospheric delay phase in the unwrapped differential interferometric phase diagram can also be ignored; the noise phase can be weakened by filtering, and the unwrapped phase of the jth differential interferometric phase of the bridge can be optimized as:
  • phase of a certain pixel of each SAR image at the time point t 1 ... t N is calculated as express, Expressed as the phase value of the pixel after each differential interference phase unwrapping, there are:
  • This step is specifically: using the least squares or singular value decomposition method to solve the equation established in step 7.1.
  • the solution is solved by the least squares method.
  • the superscript T represents the rank of the matrix.
  • the singular value decomposition method is used to solve it.
  • Step 7.3 according to the SAR imaging geometric parameters and wavelength information, the relative phase change sequence of the InSAR radar line-of-sight deformation of the bridge at the pixel (l, p) position in the image coordinate system Converted to the relative displacement change sequence of InSAR radar line - of - sight deformation variables in the image coordinate system
  • the relative displacement change sequence is geocoded to obtain the relative displacement change sequence of the InSAR radar line-of-sight deformation variable in the geographic coordinate system; it is expressed as:
  • a bridge deformation graph gr 1 is formed, as shown in Figure 5, which is a bridge deformation graph.
  • the relative displacement change value of the InSAR radar line-of-sight deformation variable of each pixel in the geographic coordinate system is drawn to form the bridge deformation map gr N .
  • Step 8 Arrange num GNSS measuring points on the bridge monitoring object.
  • each GNSS measuring point can continuously measure the three-dimensional deformation and displacement value of the location of the measuring point according to the sampling interval.
  • Deformation map gr i the corrected bridge deformation map gr' i is obtained, which represents the displacement change value of the InSAR radar line-of-sight deformation of each pixel in the geographic coordinate system at time t i relative to the initialization time.
  • the displacement values include: d Xi , d Yi , and d Ui , which represent the displacement values in the east-west, north-south, and vertical directions, respectively.
  • Step 8.2 convert the three-dimensional deformation and displacement value measured by the zth GNSS measuring point into the measured value of the radar line-of-sight deformation using the following formula:
  • ⁇ is the local incident angle of the radar wave corresponding to the SAR main image in step 2.
  • ⁇ is the satellite orbit azimuth of the SAR main image in step 2.
  • Step 8.3 for the zth GNSS measuring point, look up the bridge deformation graph gr i , and obtain the relative displacement change value of the InSAR radar line-of-sight deformation of the pixel where the zth GNSS measuring point is located.
  • Step 8.4 use the following formula to get the correction offset
  • step 8.5 the bridge deformation graph gr i is deformed as a whole to the corrected bridge deformation graph gr' i .
  • step 8 The idea of step 8 is:
  • Extract the position of the GNSS equipment that is, the radar line-of-sight deformation measurement value of the bridge InSAR time series with the same name at the location of the GNSS measuring point, and solve the correction value of the radar line-of-sight deformation measurement value to the InSAR bridge deformation map of the bridge at each time.
  • the InSAR bridge deformation map of the bridge at each time is corrected to obtain the bridge deformation distribution and deformation history after InSAR correction.
  • the bridge deformation distribution and deformation history integrated with GNSS monitoring data can be obtained.
  • the invention is a bridge deformation monitoring method that integrates GNSS data and InSAR technology, utilizes the characteristics that bridges still retain good coherence in a long period of time, adopts high-precision DSM and high-resolution SAR images to obtain bridge time-series deformation distribution, and uses
  • the GNSS simultaneous deformation monitoring results correct the bridge InSAR time series results, which overcomes the low spatial resolution of deformation monitoring caused by the sparse layout of traditional GNSS measurement, leveling and other bridge deformation monitoring methods, and the InSAR technical means obtains relatively unknown deformation reference points
  • the deformation amount and the linear deformation model of the model have poor applicability in bridge deformation monitoring, which greatly improves the spatial resolution of bridge infrastructure deformation monitoring, and provides more comprehensive data support for bridge safety operation. In the field of bridge deformation monitoring Broad application prospects.

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Image Processing (AREA)

Abstract

本发明提供一种融合GNSS数据与InSAR技术的桥梁变形监测方法，包括：获取覆盖桥梁监测对象的高分辨率SAR影像；经配准后，生成干涉相位图集合；对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图；计算得到在地理坐标系下的InSAR雷达视线向形变量相对位移变化序列；采用各个GNSS测点测得的三维形变位移值，校正桥梁形变图，得到校正后的桥梁形变图。本发明将GNSS的已有测量结果和InSAR高空间分辨率形变监测能力结合，极大改善了桥梁基础设施形变监测的空间分辨率。

Description

融合GNSS数据与InSAR技术的桥梁变形监测方法
本申请要求于2021年8月10日提交中国专利局、申请号为202110911658.8、发明名称为“融合GNSS数据与InSAR技术的桥梁变形监测方法”的中国专利申请的优先权，其全部内容通过引用结合在本申请中。
技术领域
本发明属于桥梁变形监测技术领域，具体涉及一种融合GNSS数据与InSAR技术的桥梁变形监测方法。
背景技术
桥梁作为交通基础设施的重要组成部分，已经成为城市经济发展不可缺少的一部分。但桥梁本身重量大、荷载复杂多变，地质条件变化、超载、恶劣天气等因素会对桥梁产生诸多不良影响，给桥梁安全带来隐患，影响行车安全，甚至威胁人们生命财产安全。因此，加强对桥梁的变形监测，对桥梁安全运营具有重要的意义。
传统桥梁形变测量方式，如GNSS、水准测量等方式，虽然精度较高，但需要现场布设设备，耗费人力物力、成本较高，且监测结果为稀疏点形式，不利于对桥梁开展整体分析。合成孔径雷达干涉测量技术，利用SAR影像的相位信息实现对地表变形高精度的监测，其形变监测精度可以达到厘米甚至毫米量级；同时，SAR卫星全天时全天候的影像获取能力，以及InSAR大范围、高密度、长时序、可回溯的技术特点，相比传统GNSS、水准测量等方式具有独特的优势。
但是，在InSAR技术应用过程中，获取的往往是相对于未知形变情况的参考点的相对形变量；同时，由于桥梁的变形影响因素较多，其在时间上往往呈现出非线性的特征，而InSAR技术中通常应用线性形变模型求解形变，导致结果可能会引入较大的误差。
发明内容
针对现有技术存在的缺陷，本发明提供一种融合GNSS数据与InSAR技术的桥梁变形监测方法，可有效解决上述问题。
本发明采用的技术方案如下：
本发明提供一种融合GNSS数据与InSAR技术的桥梁变形监测方法，包括以下步骤：
步骤1，确定桥梁监测对象；获取覆盖桥梁监测对象的g幅高分辨率SAR影像；
其中，g幅高分辨率SAR影像为不同监测时刻得到的SAR影像；每幅高分辨率SAR影像具有成像时间属性、影像后向散射信息和成像轨道信息；
步骤2，在各幅高分辨率SAR影像中，确定一幅高分辨率SAR影像作为SAR主影像，其他g-1幅高分辨率SAR影像作为SAR从影像；
步骤3，根据各幅高分辨率SAR影像的成像轨道信息和影像后向散射信息，对各幅SAR从影像进行旋转缩放采样操作，实现各幅SAR从影像与SAR主影像的影像配准，由此得到g-1幅配准后的SAR从影像；
步骤4，对于SAR主影像和g-1幅配准后的SAR从影像，任意两幅SAR影像形成一个干涉像对，由此生成所有干涉像对；
对于每个干涉像对，其包含的两幅SAR影像共轭相乘，生成一幅干涉相位图；
因此，对于所有干涉像对，生成多幅干涉相位图，由此组成干涉相位图集合；
步骤5，对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图；
由此将干涉相位图集合，转化为解缠差分干涉相位图集合；
步骤6，对于步骤5得到的解缠差分干涉相位图集合，假设共有M幅解缠差分干涉相位图，分别表示为：解缠差分干涉相位图tu 1，解缠差分干涉相位图tu 2，...，解缠差分干涉相位图tu M；
其中：对于任意解缠差分干涉相位图tu j，其中，j＝1，2，...，M，对应的干涉像对表示为：SAR影像I j0和SAR影像I j1，即：SAR影像I j0和SAR影像I j1共轭相乘后，经过预处理、掩膜处理和相位 解缠处理后，得到解缠差分干涉相位图tu j；
对于M幅解缠差分干涉相位图，共对应M个干涉像对，即：不考虑干涉像对包含共同SAR影像时，M个干涉像对包含2M幅SAR影像；再对2M幅SAR影像去重处理，假设最终得到N+1幅SAR影像；
将N+1幅SAR影像，按成像时间排序，依次表示为：SAR影像I 0，SAR影像I 1，...，SAR影像I N，对应的成像时间依次为：t 0，t 1，...，t N；即：SAR影像I 0的成像时间为t 0，SAR影像I 1的成像时间为t 1，...，SAR影像I N的成像时间为t N；
因此，将M幅解缠差分干涉相位图中的任意解缠差分干涉相位图，表示为解缠差分干涉相位图tu j(I c，I d)，其中，c＞d，c≠d，c＝1，2，...，N；d＝0，1，2，...，N-1，含义为：解缠差分干涉相位图tu j(I c，I d)，是由SAR影像I c和SAR影像I d共轭相乘后，经过预处理、掩膜处理和相位解缠处理后生成，解缠差分干涉相位图tu j(I c，I d)中，任意一个像元的解缠相位差，可直接通过解缠差分干涉相位图读取获得；
步骤7，M幅解缠差分干涉相位图和N+1幅SAR影像的图像大小完全相同，具有相同数量的像元；
对于任意解缠差分干涉相位图中的像元(l，p)，其中，l，p分别为像元的行号和列号，执行步骤7.1-步骤7.3，计算得到像元(l，p)对应的在地理坐标系下的InSAR雷达视线向形变量相对位移变化序列：
步骤7.1，建立以下方程：
Figure PCTCN2022106155-appb-000001
其中：
Figure PCTCN2022106155-appb-000002
分别表示SAR影像I 1相对于SAR影像I 0在像元(l，p)位置的解缠相位差，SAR影像I 2相对于SAR影像I 0在像元(l，p)位置的解缠相位差，...，SAR影像I N相 对于SAR影像I 0在像元(l，p)位置的解缠相位差；
Figure PCTCN2022106155-appb-000003
Figure PCTCN2022106155-appb-000004
为未知待求值；
Figure PCTCN2022106155-appb-000005
分别表示解缠差分干涉相位图tu 1在像元(l，p)位置的解缠相位差，解缠差分干涉相位图tu 2在像元(l，p)位置的解缠相位差，...，解缠差分干涉相位图tu M在像元(l，p)位置的解缠相位差；
Figure PCTCN2022106155-appb-000006
可通过对应的解缠差分干涉相位图直接读取到；
Figure PCTCN2022106155-appb-000007
为M*N的系数矩阵A，为已知值，通过以下方式确定：
假设对于系数矩阵A中的任意第j行，其中，j＝1，2，...，M，则有以下关系式：
Figure PCTCN2022106155-appb-000008
Figure PCTCN2022106155-appb-000009
为解缠差分干涉相位图tu j(I c，I d)在像元(l，p)位置的解缠相位差；解缠差分干涉相位图tu j(I c，I d)对应SAR影像I c和SAR影像I d；
如果d＝0，则：在
Figure PCTCN2022106155-appb-000010
中，令
Figure PCTCN2022106155-appb-000011
前面的系数a jc为1，其他N-1个系数均为0，由此确定了系数a j1，...，a jN的值；
如果d≠0，则：在
Figure PCTCN2022106155-appb-000012
中，令
Figure PCTCN2022106155-appb-000013
前面的系数a jc为1，令
Figure PCTCN2022106155-appb-000014
前面的系数a jd为-1，其他N-2个系数均为0，由此确定了系数a j1，...，a jN的值；
步骤7.2，求解步骤7.1建立的方程，得到桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序列
Figure PCTCN2022106155-appb-000015
即：桥梁在像元(l，p)位置，相对于t 0时刻，分别在t 1，...，t N时刻在图像坐标系下的InSAR雷达视线向形变量相对相位变化值；
步骤7.3，依据SAR成像几何参数和波长信息，将桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序 列
Figure PCTCN2022106155-appb-000016
转换为图像坐标系下InSAR雷达视线向形变量相对位移变化序列Δd(t 1)，Δd(t 2)，...，Δd(t N)，并将图像坐标系下InSAR雷达视线向形变量相对位移变化序列地理编码，得到地理坐标系下的InSAR雷达视线向形变量相对位移变化序列；表示为：
Figure PCTCN2022106155-appb-000017
通过以上方式，假设任意解缠差分干涉相位图中，桥梁监测对象覆盖N 0个像元，则t 1时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 1；
t 2时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 2；
依此类推
t N时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr N；
步骤8，在桥梁监测对象布置num个GNSS测点；
从初始化时刻开始，各个GNSS测点按采样间隔，可持续测得测点所在位置的三维形变位移值；
因此，获得t i时刻各个GNSS测点测得的三维形变位移值，其中，i＝1，2，...，N，采用t i时刻各个GNSS测点测得的三维形变位移值，校正桥梁形变图gr i，得到校正后的桥梁形变图gr’ i，表示t i时刻相对于初始化时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量位移变化值。
本发明提供的融合GNSS数据与InSAR技术的桥梁变形监测方法具有以下优点：
本发明的目的在于克服传统GNSS测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、InSAR技术手段获取的多是相对未知形变情况参考点的形变量以及InSAR常用的线性形变模型在桥梁形变监测的误差较大的不足，提出一种融合GNSS数据与InSAR技术的桥梁变形监测方法，将GNSS的已有测量结果和InSAR高空间分辨率形变监测能力结合，极大改善了桥梁基础设施形变监测的空间分辨率。
说明书附图
图1为本发明提供的融合GNSS数据与InSAR技术的桥梁变形监测方法的流程示意图；
图2为SAR影像的干涉像对组合图；
图3为处理后的差分干涉相位图的示意图；
图4为解缠结果正确的解缠差分干涉相位图示意图；
图5是桥梁InSAR形变分布图；
图6是系数矩阵确定方式的原理图。
具体实施方式
本发明的目的在于克服传统GNSS测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、InSAR技术手段获取的多是相对未知形变情况参考点的形变量以及InSAR常用的线性形变模型在桥梁形变监测的误差较大的不足，提出一种融合GNSS数据与InSAR技术的桥梁变形监测方法，将GNSS的已有测量结果和InSAR高空间分辨率形变监测能力结合，极大改善了桥梁基础设施形变监测的空间分辨率。
参考图1，本发明提供的融合GNSS数据与InSAR技术的桥梁变形监测方法，包括以下步骤：
步骤1，确定桥梁监测对象；获取覆盖桥梁监测对象的g幅高分辨率SAR(SyntheticAperture Radar，即合成孔径雷达)影像。其中，g幅高分辨率SAR影像为不同监测时刻得到的SAR影像；每幅高分辨率SAR影像具有成像时间属性、影像后向散射信息和成像轨道信息。
实际应用中，高分辨率SAR影像要求SAR数据地面像元优于3米，可以为COSMO-SkyMed雷达卫星星座、TerraSAR-X/TanDEM-X雷达卫星、RadarSAT-2雷达卫星等高分辨率雷达卫星系统的聚束模式，高分辨率成像的条带模式，为获取桥梁重复轨道的高分辨率星载SAR数据集提供数据源；同时，选取的SAR影像之间的垂直基线需要尽量短，以TerraSAR-X为例，SAR数据集垂直基线分布要求在300米以内，以尽量降低DSM高程误差带来的影响。
步骤2，在各幅高分辨率SAR影像中，确定一幅高分辨率SAR影像 作为SAR主影像，其他g-1幅高分辨率SAR影像作为SAR从影像。
步骤3，根据各幅高分辨率SAR影像的成像轨道信息和影像后向散射信息，对各幅SAR从影像进行旋转缩放采样操作，实现各幅SAR从影像与SAR主影像的影像配准，由此得到g-1幅配准后的SAR从影像。
具体的，可以基于SAR影像成像轨道信息求解主从SAR影像的配准偏移量，对主从SAR影像进行粗配准；再使用强度归一化互相关算法对主从SAR影像进行精确配准，配准精度需要优于1/8个像素。
步骤4，对于SAR主影像和g-1幅配准后的SAR从影像，任意两幅SAR影像形成一个干涉像对，由此生成所有干涉像对。对于每个干涉像对，其包含的两幅SAR影像共轭相乘，生成一幅干涉相位图。因此，对于所有干涉像对，生成多幅干涉相位图，由此组成干涉相位图集合。
具体的，忽略时间基线和空间基线，对所有高分辨率SAR影像进行两两干涉组合，生成尽可能多的干涉像对，对于每个干涉像对，其包含的两幅SAR影像共轭相乘，由此组成干涉相位图集合；参照附图2，为SAR影像的干涉像对组合图。
步骤5，对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图。由此将干涉相位图集合，转化为解缠差分干涉相位图集合。
步骤5具体为：
步骤5.1，预处理：
对每幅干涉相位图进行预处理，包括平地相位去除、地形相位去除和滤波处理，得到处理后的差分干涉相位图。
本步骤中，采用以下方法，进行地形相位去除：
获取覆盖桥梁监测对象的高精度DSM数据；基于高精度DSM数据，进行地形相位去除。其中：高精度DSM数据要求现势性强，其提供的桥梁监测对象的高程信息在桥梁建成以后；高精度DSM数据的精度和空间分辨率要求高于AW3D30 DSM数据。实际应用中，高分七号立体像对、WorldView-3立体像对、TerraSAR-X/TanDEM-X双星系统、无人机倾斜摄影测量等可为研究区提供高分辨率高精度的DSM数据。因此，利用SAR影像轨道信息和高精度DSM数据去除平地相位和地形相位，并采用 Goldstein等滤波方法对干涉相位图进行滤波，削弱噪声的影响，参照附图3，为处理后的差分干涉相位图的示意图。
步骤5.2，掩膜处理：
步骤5.2.1，制作桥梁监测对象掩膜；其中，所述桥梁监测对象掩膜是指：在掩膜图像中，桥梁监测对象的像素点为1，其他所有位置的像素点为0。
步骤5.2.2，采用桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，得到掩膜处理后的差分干涉相位图；其中，在掩膜处理后的差分干涉相位图中，桥梁监测对象的像素点相位保持不变，其他位置的像素点相位统一变为0。
具体的，通过制作桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，使差分干涉相位图仅保留桥梁监测对象的影像，降低桥梁外其他背景影像的影响。
步骤5.3，相位解缠处理：
对于每幅掩膜处理后的差分干涉相位图，均采用多种相位解缠算法进行相位解缠运算，由此得到多幅解缠差分干涉相位图，并采用筛选算法，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图。
将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，最终得到多个解缠差分干涉相位图。
具体的，对于每幅掩膜处理后的差分干涉相位图，在桥梁相位稳定区域选取参考点，分别采用最小费用流法、枝切法等多种相位解缠算法，对所有差分干涉相位图进行相位解缠。
然后，通过人工目视解译或剖面线分析等方法检查解缠差分干涉相位图是否存在解缠误差，比如明显的相位跳跃，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图。
然后，将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，相同差分干涉组合保留一幅，最大化相位解缠结果正确的解缠差分干涉相位图，参照附图4，为解缠结果正确的解缠差分干涉相位图示意图。
步骤6，对于步骤5得到的解缠差分干涉相位图集合，假设共有M幅解缠差分干涉相位图，分别表示为：解缠差分干涉相位图tu 1，解缠差分干涉相位图tu 2，...，解缠差分干涉相位图tu M。
其中：对于任意解缠差分干涉相位图tu j，其中，j＝1，2，...，M，对应的干涉像对表示为：SAR影像I j0和SAR影像I j1，即：SAR影像I j0和SAR影像I j1共轭相乘后，经过预处理、掩膜处理和相位解缠处理后，得到解缠差分干涉相位图tu j。
对于M幅解缠差分干涉相位图，共对应M个干涉像对，即：不考虑干涉像对包含共同SAR影像时，M个干涉像对包含2M幅SAR影像；再对2M幅SAR影像去重处理，假设最终得到N+1幅SAR影像。
将N+1幅SAR影像，按成像时间排序，依次表示为：SAR影像I 0，SAR影像I 1，...，SAR影像I N，对应的成像时间依次为：t 0，t 1，...，t N；即：SAR影像I 0的成像时间为t 0，SAR影像I 1的成像时间为t 1，...，SAR影像I N的成像时间为t N。
因此，将M幅解缠差分干涉相位图中的任意解缠差分干涉相位图，表示为解缠差分干涉相位图tu j(I c，I d)，其中，c＞d，c≠d，c＝1，2，...，N；d＝0，1，2，...，N-1，含义为：解缠差分干涉相位图tu j(I c，I d)，是由SAR影像I c和SAR影像I d共轭相乘后，经过预处理、掩膜处理和相位解缠处理后生成，解缠差分干涉相位图tu j(I c，I d)中，任意一个像元的解缠相位差，可直接通过解缠差分干涉相位图读取获得。
步骤7，M幅解缠差分干涉相位图和N+1幅SAR影像的图像大小完全相同，具有相同数量的像元。
对于任意解缠差分干涉相位图中的像元(l，p)，其中，l，p分别为像元的行号和列号，执行步骤7.1-步骤7.3，计算得到像元(l，p)对应的在地理坐标系下的InSAR雷达视线向形变量相对位移变化序列：
步骤7.1，建立以下方程：
Figure PCTCN2022106155-appb-000018
其中：
Figure PCTCN2022106155-appb-000019
分别表示SAR影像I 1相对于SAR影像I 0在像元(l，p)位置的解缠相位差，SAR影像I 2相对于SAR影像I 0在像元(l，p)位置的解缠相位差，...，SAR影像IN相对于SAR影像I 0在像元(l，p)位置的解缠相位差；
Figure PCTCN2022106155-appb-000020
为未知待求值。
Figure PCTCN2022106155-appb-000021
分别表示解缠差分干涉相位图tu 1在像元(l，p)位置的解缠相位差，解缠差分干涉相位图tu 2在像元(l，p)位置的解缠相位差，...，解缠差分干涉相位图tu M在像元(l，p)位置的解缠相位差；
Figure PCTCN2022106155-appb-000022
Figure PCTCN2022106155-appb-000023
可通过对应的解缠差分干涉相位图直接读取到。
Figure PCTCN2022106155-appb-000024
为M*N的系数矩阵A，为已知值，通过以下方式确定：
假设对于系数矩阵A中的任意第j行，其中，j＝1，2，...，M，则有以下关系式：
Figure PCTCN2022106155-appb-000025
Figure PCTCN2022106155-appb-000026
为解缠差分干涉相位图tu j(I c，I d)在像元(l，p)位置的解缠相位差；解缠差分干涉相位图tu j(I c，I d)对应SAR影像I c和SAR影像I d。
如果d＝0，则：在
Figure PCTCN2022106155-appb-000027
中，令
Figure PCTCN2022106155-appb-000028
前面的系数a jc为1，其他N-1个系数均为0，由此确定了系数a j1，...，a jN的值。
如果d≠0，则：在
Figure PCTCN2022106155-appb-000029
中，令
Figure PCTCN2022106155-appb-000030
前面的系数a jc为1，令
Figure PCTCN2022106155-appb-000031
前面的系数a jd为-1，其他N-2个系数均为0，由此确定了系数a j1，...，a jN的值。
步骤7.2，求解步骤7.1建立的方程，得到桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序列
Figure PCTCN2022106155-appb-000032
Figure PCTCN2022106155-appb-000033
即：桥梁在像元(l，p)位置，相对于t 0时刻，分别在t 1，...，t N时刻在图像坐标系下的InSAR雷达视线向形变量相对相位变化值。
为方便对步骤6和步骤7进行理解，下面列举一个实施例：
参考图6：
假设经步骤3处理后，共有4幅SAR影像，分别为SAR影像I0，SAR影像I1，SAR影像I2，SAR影像I3。
在步骤4中，任意两幅SAR影像形成一个干涉像对，由此生成6个干涉像对。对于每个干涉像对，生成一幅干涉相位图，因此，共得到6幅干涉相位图。
在步骤5中，通过对每幅干涉相位图进行预处理、掩膜处理和相位解缠处理，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分干涉相位图，因此，本步骤中，最后得到3幅解缠差分干涉相位图，分别表示为：
解缠差分干涉相位图tu 3，对应的干涉像对为：SAR影像I0和SAR影像I3；其任意像元的解缠相位差表示为
Figure PCTCN2022106155-appb-000034
为已知值。
解缠差分干涉相位图tu 4，对应的干涉像对为：SAR影像I0和SAR影像I2；其任意像元的解缠相位差表示为
Figure PCTCN2022106155-appb-000035
为已知值。
解缠差分干涉相位图tu 5，对应的干涉像对为：SAR影像I1和SAR影像I3；其任意像元的解缠相位差表示为
Figure PCTCN2022106155-appb-000036
为已知值。
因此，共得到M＝3幅解缠差分干涉相位图，对应的SAR影像分别为I0、I1、I2和I3，为4幅，即：N+1＝4，N＝3。
假设4幅SAR影像，按成像时间排序，即为：SAR影像I0、SAR影像I1、SAR影像I2和SAR影像I3。
由此建立以下方程：
Figure PCTCN2022106155-appb-000037
其中：
Figure PCTCN2022106155-appb-000038
为未知待求值；
上述方程可转化为以下方程组：
Figure PCTCN2022106155-appb-000039
Figure PCTCN2022106155-appb-000040
Figure PCTCN2022106155-appb-000041
其中：
解缠差分干涉相位图tu 3(I c，I d)＝tu 3(I 3，I 0)，d＝0，因此，a 13＝1，a 11和a 12为0。
解缠差分干涉相位图tu 4(I c，I d)＝tu 3(I 2，I 1)，d≠0，因此，a 22＝1，a 21＝-1，a 23为0。
解缠差分干涉相位图tu 5(I c，I d)＝tu 3(I 3，I 1)，d≠0，因此，a 33＝1，a 31＝-1，a 32为0。
由此确定系数矩阵。
然后，根据方程组中未知待求值与已知值之间的数量关系，选择最小二乘或奇异值分解方法求解。
上述系数确定的原理可以理解为：
假设有M幅解缠差分干涉相位图解缠正确，涉及到的SAR影像N+1幅。解缠差分干涉相位图上像元位置可用方位向坐标l和距离向坐标p来描述，第j幅解缠差分干涉相位图在像元(l，p)的解缠相位
Figure PCTCN2022106155-appb-000042
可表示为：
Figure PCTCN2022106155-appb-000043
其中：
第j幅解缠差分干涉相位图，是由成像时间为t B的SAR影像和成像时间为t A的SAR影像，经过共轭、预处理、掩膜处理和相位解缠处理后生成。
j∈[1，M]，d(t B，l，p)和d(t A，l，p)为像元(l，p)在t B和t A时间相对于SAR数据集起始时间雷达视线向的形变量。
λ为雷达波长，r为像元(l，p)与卫星之间的距离，θ为像元局部入射角，B j⊥为干涉对的时间基线，Δz为高程残差，
Figure PCTCN2022106155-appb-000044
和
Figure PCTCN2022106155-appb-000045
为像元(l，p)在t B和t A时间的大气延迟相位，Δn j为噪声相位。
当使用的DSM高程精度较高时，且SAR影像干涉像对垂直基线较 短时，由高程残差引起的高程残差相位可以忽略；同时，由于桥梁覆盖范围往往较小，可认为该范围内大气分布相对较均匀，因而解缠差分干涉相位图中的大气延迟相位也可以忽略；噪声相位可通过滤波进行削弱，则针对桥梁第j幅差分干涉相位图解缠后的相位可优化为：
Figure PCTCN2022106155-appb-000046
因此，t 1...t N时间点上各SAR影像某像元的相位用
Figure PCTCN2022106155-appb-000047
表示，
Figure PCTCN2022106155-appb-000048
表示为各差分干涉相位图解缠后该像元的相位值，则有：
Figure PCTCN2022106155-appb-000049
Figure PCTCN2022106155-appb-000050
Figure PCTCN2022106155-appb-000051
知
Figure PCTCN2022106155-appb-000052
可组合为：
Figure PCTCN2022106155-appb-000053
其中
Figure PCTCN2022106155-appb-000054
为未知量，
Figure PCTCN2022106155-appb-000055
为已知量，A为一个M×N的矩阵，每一行与
Figure PCTCN2022106155-appb-000056
相乘对应着一幅解缠差分干涉相位图。
Figure PCTCN2022106155-appb-000057
为M个差分干涉对，包含M个方程，
Figure PCTCN2022106155-appb-000058
包含N个未知的SAR影像对应的原始相位值。
本步骤具体为：采用最小二乘或奇异值分解方法，求解步骤7.1建立的方程。
具体的，当M≥N，且矩阵A的秩为N，通过最小二乘法求解。
Figure PCTCN2022106155-appb-000059
其中：
Figure PCTCN2022106155-appb-000060
Figure PCTCN2022106155-appb-000061
上标T代表矩阵的转秩。
当解缠差分干涉相位图集分为L组，并且，矩阵A的秩小于N，为N-L+1，使用奇异值分解方法求解。
步骤7.3，依据SAR成像几何参数和波长信息，将桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序列
Figure PCTCN2022106155-appb-000062
转换为图像坐标系下InSAR雷达视线向形变量相对位移变化序列Δd(t 1)，Δd(t 2)，...，Δd(t N)，并将图像坐标系下InSAR雷达视线向形变量相对位移变化序列地理编码，得到地理坐标系下的InSAR雷达视线向形变量相对位移变化序列；表示为：
Figure PCTCN2022106155-appb-000063
通过以上方式，假设任意解缠差分干涉相位图中，桥梁监测对象覆盖N 0个像元，则t 1时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 1，如图5所示，为桥梁形变图。
t 2时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 2。
依此类推
t N时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr N。
步骤8，在桥梁监测对象布置num个GNSS测点。
从初始化时刻开始，各个GNSS测点按采样间隔，可持续测得测点所在位置的三维形变位移值。
因此，获得t i时刻各个GNSS测点测得的三维形变位移值，其中，i＝1，2，...，N，采用t i时刻各个GNSS测点测得的三维形变位移值，校正桥梁形变图gr i，得到校正后的桥梁形变图gr’ i，表示t i时刻相对于初始化时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量位移变化值。
具体采用以下方法，对桥梁形变图gr i进行校正：
步骤8.1，获得t i时刻各个GNSS测点测得的三维形变位移值；其中，对于任意第z个GNSS测点，z＝1，2，...，num，其在t i时刻的三维形变位移值包括：d Xi、d Yi、d Ui，分别代表东西向、南北向和垂向上的位移值。
步骤8.2，将第z个GNSS测点测得的三维形变位移值，采用以下公式，转化为雷达视线向形变测量值：
Figure PCTCN2022106155-appb-000064
其中：
θ为步骤2中SAR主影像对应的雷达波局部入射角。
α为步骤2中SAR主影像的卫星轨道方位角。
Figure PCTCN2022106155-appb-000065
在t i时刻，第z个GNSS测点测得的雷达视线向形变测量值。
步骤8.3，对于第z个GNSS测点，查找桥梁形变图gr i，获得第z个GNSS测点所在位置的像元的InSAR雷达视线向形变量相对位移变化值
Figure PCTCN2022106155-appb-000066
步骤8.4，采用下式，得到校正偏移量
Figure PCTCN2022106155-appb-000067
Figure PCTCN2022106155-appb-000068
步骤8.5，使桥梁形变图gr i整体发生形变得到校正后的桥梁形变图gr’ i。
步骤8的构思为：
提取GNSS设备所在位置，即GNSS测点所在位置的桥梁InSAR时序同名点的雷达视线向形变测量值，求解得到雷达视线向形变测量值对桥梁各个时间的InSAR桥梁形变图的校正量值，并对桥梁各时间的InSAR桥梁形变图进行校正，获取InSAR校正后的桥梁形变分布和形变历史。
基于以上步骤，可以获取融合了GNSS监测数据的桥梁形变分布和形变历史。
本发明提供的融合GNSS数据与InSAR技术的桥梁变形监测方法具有以下优点：
本发明为一种融合GNSS数据与InSAR技术的桥梁变形监测方法，利用桥梁在长时间段相干性仍保留较好的特点，采用高精度DSM和高分辨率SAR影像获取桥梁时间序列形变分布，利用GNSS同期形变监测结果对桥梁InSAR时序结果进行校正，克服了传统GNSS测量、水准测量等桥梁形变监测手段稀疏布设导致的形变监测空间分辨率低、InSAR技术手段获取的是相对未知形变情况的参考点的变形量和线性形变模型在桥梁形变监测适用性较差的不足，极大改善了桥梁基础设施形变监测的空间分辨率，为桥梁安全运营提供了更加全面的数据支撑，在桥梁形变监测领域有着广阔的应用前景。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (6)

  1. 一种融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，包括以下步骤：
    步骤1，确定桥梁监测对象；获取覆盖桥梁监测对象的g幅高分辨率SAR影像；
    其中，g幅高分辨率SAR影像为不同监测时刻得到的SAR影像；每幅高分辨率SAR影像具有成像时间属性、影像后向散射信息和成像轨道信息；
    步骤2，在各幅高分辨率SAR影像中，确定一幅高分辨率SAR影像作为SAR主影像，其他g-1幅高分辨率SAR影像作为SAR从影像；
    步骤3，根据各幅高分辨率SAR影像的成像轨道信息和影像后向散射信息，对各幅SAR从影像进行旋转缩放采样操作，实现各幅SAR从影像与SAR主影像的影像配准，由此得到g-1幅配准后的SAR从影像；
    步骤4，对于SAR主影像和g-1幅配准后的SAR从影像，任意两幅SAR影像形成一个干涉像对，由此生成所有干涉像对；
    对于每个干涉像对，其包含的两幅SAR影像共轭相乘，生成一幅干涉相位图；
    因此，对于所有干涉像对，生成多幅干涉相位图，由此组成干涉相位图集合；
    步骤5，对于干涉相位图集合中的每幅干涉相位图，进行预处理、掩膜处理和相位解缠处理，得到解缠差分干涉相位图；
    由此将干涉相位图集合，转化为解缠差分干涉相位图集合；
    步骤6，对于步骤5得到的解缠差分干涉相位图集合，假设共有M幅解缠差分干涉相位图，分别表示为：解缠差分干涉相位图tu 1，解缠差分干涉相位图tu 2，...，解缠差分干涉相位图tu M；
    其中：对于任意解缠差分干涉相位图tu j，其中，j＝1，2，...，M，对应的干涉像对表示为：SAR影像I j0和SAR影像I j1，即：SAR影像I j0和SAR影像I j1共轭相乘后，经过预处理、掩膜处理和相位解缠处理后，得到解缠差分干涉相位图tu j；
    对于M幅解缠差分干涉相位图，共对应M个干涉像对，即：不考虑干涉像对包含共同SAR影像时，M个干涉像对包含2M幅SAR影像；再 对2M幅SAR影像去重处理，假设最终得到N+1幅SAR影像；
    将N+1幅SAR影像，按成像时间排序，依次表示为：SAR影像I 0，SAR影像I 1，...，SAR影像I N，对应的成像时间依次为：t 0，t 1，...，t N；即：SAR影像I 0的成像时间为t 0，SAR影像I 1的成像时间为t 1，...，SAR影像I N的成像时间为t N；
    因此，将M幅解缠差分干涉相位图中的任意解缠差分干涉相位图，表示为解缠差分干涉相位图tu j(I c，I d)，其中，c＞d，c≠d，c＝1，2，...，N；d＝0，1，2，...，N-1，含义为：解缠差分干涉相位图tu j(I c，I d)，是由SAR影像I c和SAR影像I d共轭相乘后，经过预处理、掩膜处理和相位解缠处理后生成，解缠差分干涉相位图tu j(I c，I d)中，任意一个像元的解缠相位差，可直接通过解缠差分干涉相位图读取获得；
    步骤7，M幅解缠差分干涉相位图和N+1幅SAR影像的图像大小完全相同，具有相同数量的像元；
    对于任意解缠差分干涉相位图中的像元(l，p)，其中，l，p分别为像元的行号和列号，执行步骤7.1-步骤7.3，计算得到像元(l，p)对应的在地理坐标系下的InSAR雷达视线向形变量相对位移变化序列：
    步骤7.1，建立以下方程：
    Figure PCTCN2022106155-appb-100001
    其中：
    Figure PCTCN2022106155-appb-100002
    分别表示SAR影像I 1相对于SAR影像I 0在像元(l，p)位置的解缠相位差，SAR影像I 2相对于SAR影像I 0在像元(l，p)位置的解缠相位差，...，SAR影像I N相对于SAR影像I 0在像元(l，p)位置的解缠相位差；
    Figure PCTCN2022106155-appb-100003
    为未知待求值；
    Figure PCTCN2022106155-appb-100004
    分别表示解缠差分干涉相位图tu 1在像元(l，p)位置的解缠相位差，解缠差分干涉相位图tu 2在像元(l，p)位置的解缠相位差，...，解缠差分干涉相位图tu M在像元(l，p)位置的解缠相位差；
    Figure PCTCN2022106155-appb-100005
    Figure PCTCN2022106155-appb-100006
    可通过对应的解缠差分干涉相位图直接读取到；
    Figure PCTCN2022106155-appb-100007
    为M*N的系数矩阵A，为已知值，通过以下方式确定：
    假设对于系数矩阵A中的任意第j行，其中，j＝1，2，...，M，则有以下关系式：
    Figure PCTCN2022106155-appb-100008
    Figure PCTCN2022106155-appb-100009
    为解缠差分干涉相位图tu j(I c，I d)在像元(l，p)位置的解缠相位差；解缠差分干涉相位图tu j(I c，I d)对应SAR影像I c和SAR影像I d；
    如果d＝0，则：在
    Figure PCTCN2022106155-appb-100010
    中，令
    Figure PCTCN2022106155-appb-100011
    前面的系数a jc为1，其他N-1个系数均为0，由此确定了系数a j1，...，a jN的值；
    如果d≠0，则：在
    Figure PCTCN2022106155-appb-100012
    中，令
    Figure PCTCN2022106155-appb-100013
    前面的系数a jc为1，令
    Figure PCTCN2022106155-appb-100014
    前面的系数a jd为-1，其他N-2个系数均为0，由此确定了系数a j1，...，a jN的值；
    步骤7.2，求解步骤7.1建立的方程，得到桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序列
    Figure PCTCN2022106155-appb-100015
    Figure PCTCN2022106155-appb-100016
    即：桥梁在像元(l，p)位置，相对于t 0时刻，分别在t 1，...，t N时刻在图像坐标系下的InSAR雷达视线向形变量相对相位变化值；
    步骤7.3，依据SAR成像几何参数和波长信息，将桥梁在像元(l，p)位置在图像坐标系下的InSAR雷达视线向形变量相对相位变化序列
    Figure PCTCN2022106155-appb-100017
    转换为图像坐标系下InSAR雷达视线向形变量相对位移变化序列Δd(t 1)，Δd(t 2)，...，Δd(t N)，并将图像坐标系下InSAR雷达视线向形变量相对位移变化序列地理编码，得到地理坐标系下的InSAR雷达视线向形变量相对位移变化序列；表示为：
    Figure PCTCN2022106155-appb-100018
    通过以上方式，假设任意解缠差分干涉相位图中，桥梁监测对象覆盖N 0个像元，则t 1时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 1；
    t 2时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr 2；
    依此类推
    t N时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量相对位移变化值，绘制形成桥梁形变图gr N；
    步骤8，在桥梁监测对象布置num个GNSS测点；
    从初始化时刻开始，各个GNSS测点按采样间隔，可持续测得测点所在位置的三维形变位移值；
    因此，获得t i时刻各个GNSS测点测得的三维形变位移值，其中，i＝1，2，...，N，采用t i时刻各个GNSS测点测得的三维形变位移值，校正桥梁形变图gr i，得到校正后的桥梁形变图gr’ i，表示t i时刻相对于初始化时刻，各个像元在地理坐标系下的InSAR雷达视线向形变量位移变化值。
  2. 根据权利要求1所述的融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，步骤5具体为：
    步骤5.1，预处理：
    对每幅干涉相位图进行预处理，包括平地相位去除、地形相位去除和滤波处理，得到处理后的差分干涉相位图；
    步骤5.2，掩膜处理：
    步骤5.2.1，制作桥梁监测对象掩膜；其中，所述桥梁监测对象掩膜是指：在掩膜图像中，桥梁监测对象的像素点为1，其他所有位置的像素点为0；
    步骤5.2.2，采用桥梁监测对象掩膜，对每幅处理后的差分干涉相位图进行掩膜处理，得到掩膜处理后的差分干涉相位图；其中，在掩膜处理后的差分干涉相位图中，桥梁监测对象的像素点相位保持不变，其他位置的像素点相位统一变为0；
    步骤5.3，相位解缠处理：
    对于每幅掩膜处理后的差分干涉相位图，均采用多种相位解缠算法进行相位解缠运算，由此得到多幅解缠差分干涉相位图，并采用筛选算法，剔除存在解缠误差的差分干涉相位图，保留相位解缠结果正确的解缠差分 干涉相位图；
    将多种相位解缠算法获取的解缠结果正确的解缠差分干涉相位图合并，最终得到多个解缠差分干涉相位图。
  3. 根据权利要求1所述的融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，步骤5.1中，采用以下方法，进行地形相位去除：
    获取覆盖桥梁监测对象的高精度DSM数据；
    基于高精度DSM数据，进行地形相位去除。
  4. 根据权利要求1所述的融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，步骤7.2中，采用最小二乘或奇异值分解方法，求解步骤7.1建立的方程。
  5. 根据权利要求4所述的融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，当M≥N，且矩阵A的秩为N，通过最小二乘法求解；
    当解缠差分干涉相位图集分为上组，并且，矩阵A的秩小于N，为N-L+1，使用奇异值分解方法求解。
  6. 根据权利要求1所述的融合GNSS数据与InSAR技术的桥梁变形监测方法，其特征在于，步骤8中，具体采用以下方法，对桥梁形变图gr i进行校正：
    步骤8.1，获得t i时刻各个GNSS测点测得的三维形变位移值；其中，对于任意第z个GNSS测点，z＝1，2，...，num，其在t i时刻的三维形变位移值包括：d Xi、d Yi、d Ui，分别代表东西向、南北向和垂向上的位移值；
    步骤8.2，将第z个GNSS测点测得的三维形变位移值，采用以下公式，转化为雷达视线向形变测量值：
    Figure PCTCN2022106155-appb-100019
    θ为步骤2中SAR主影像对应的雷达波局部入射角；
    α为步骤2中SAR主影像的卫星轨道方位角；
    Figure PCTCN2022106155-appb-100020
    在t i时刻，第z个GNSS测点测得的雷达视线向形变测量值；
    步骤8.3，对于第z个GNSS测点，查找桥梁形变图gr i，获得第z个GNSS测点所在位置的像元的InSAR雷达视线向形变量相对位移变化值
    Figure PCTCN2022106155-appb-100021
    步骤8.4，采用下式，得到校正偏移量
    Figure PCTCN2022106155-appb-100022
    Figure PCTCN2022106155-appb-100023
    步骤8.5，使桥梁形变图gr i整体发生形变
    Figure PCTCN2022106155-appb-100024
    得到校正后的桥梁形变图gr’ i。
PCT/CN2022/106155 2021-08-10 2022-07-18 融合GNSS数据与InSAR技术的桥梁变形监测方法 WO2022214114A2 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
ZA2023/00146A ZA202300146B (en) 2021-08-10 2023-01-03 Bridge deformation monitoring method fusing gnss data and insar technology

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202110911658.8 2021-08-10
CN202110911658.8A CN113624122B (zh) 2021-08-10 2021-08-10 融合GNSS数据与InSAR技术的桥梁变形监测方法

Publications (2)

Publication Number Publication Date
WO2022214114A2 true WO2022214114A2 (zh) 2022-10-13
WO2022214114A3 WO2022214114A3 (zh) 2022-11-10

Family

ID=78383877

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2022/106155 WO2022214114A2 (zh) 2021-08-10 2022-07-18 融合GNSS数据与InSAR技术的桥梁变形监测方法

Country Status (3)

Country Link
CN (1) CN113624122B (zh)
WO (1) WO2022214114A2 (zh)
ZA (1) ZA202300146B (zh)

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115451874A (zh) * 2022-11-11 2022-12-09 深圳市科荣软件股份有限公司 一种基于InSAR图像的水坝形变监测方法及装置
CN115453520A (zh) * 2022-10-26 2022-12-09 中国电子科技集团公司第十四研究所 基于双频多极化差分干涉的地表形变测量方法及设备
CN116150301A (zh) * 2023-04-24 2023-05-23 北京中舆达信息科技有限公司 一种基于3DWebGIS的地表形变监测系统及方法
CN116224327A (zh) * 2023-02-20 2023-06-06 中国矿业大学 一种基于学习网络的矿区大梯度形变区相位解缠方法
CN116258005A (zh) * 2023-02-23 2023-06-13 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN116338690A (zh) * 2023-03-28 2023-06-27 中南林业科技大学 一种无模型约束的时序InSAR地形残差与地表形变估计方法
CN116736305A (zh) * 2023-08-14 2023-09-12 北京观微科技有限公司 地物形变确定方法、装置、电子设备及存储介质
CN117036965A (zh) * 2023-10-08 2023-11-10 四川正路建设工程检测咨询有限公司 桥梁维修设备控制方法、电子设备和计算机可读介质
CN117213443A (zh) * 2023-11-07 2023-12-12 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法
CN117274342A (zh) * 2023-11-21 2023-12-22 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法
CN117368916A (zh) * 2023-10-12 2024-01-09 云南大学 一种InSAR相位解缠方法、系统、设备及介质
CN117435902A (zh) * 2023-12-20 2024-01-23 武汉华威科智能技术有限公司 一种基于机器学习确定rfid标签运动行为的方法及装置
CN117724102A (zh) * 2024-02-18 2024-03-19 中国特种设备检测研究院 一种结合ekfpu的mcf相位解缠方法及系统

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113624122B (zh) * 2021-08-10 2022-09-20 中咨数据有限公司 融合GNSS数据与InSAR技术的桥梁变形监测方法

Family Cites Families (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101770027B (zh) * 2010-02-05 2012-05-16 河海大学 基于InSAR与GPS数据融合的地表三维形变监测方法
CN106772342B (zh) * 2017-01-11 2020-02-07 西南石油大学 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN107621636B (zh) * 2017-11-13 2021-08-06 河海大学 一种基于psi的大型铁路桥梁健康监测方法
CN207622735U (zh) * 2017-12-22 2018-07-17 交通运输部科学研究院 一种基于北斗gnss的桥梁变形监测系统
CN109212522B (zh) * 2018-05-28 2020-01-21 中国科学院电子学研究所 一种基于双基星载sar的高精度dem反演方法及装置
CN109375222B (zh) * 2018-12-17 2019-12-06 中国自然资源航空物探遥感中心 一种合成孔径雷达干涉测量电离层相位估计和补偿方法
US10852421B1 (en) * 2019-01-24 2020-12-01 Descartes Labs, Inc. Sparse phase unwrapping
CN109782282B (zh) * 2019-03-13 2023-04-18 武汉大学 一种集成对流层大气延迟改正的时间序列InSAR分析方法
US10895637B1 (en) * 2019-07-17 2021-01-19 BGA Technology LLC Systems and methods for mapping manmade objects buried in subterranean surfaces using an unmanned aerial vehicle integrated with radar sensor equipment
CN111174689A (zh) * 2020-03-04 2020-05-19 广东明源勘测设计有限公司 一种桥梁形变监测方法
CN111998766B (zh) * 2020-08-31 2021-10-15 同济大学 一种基于时序InSAR技术的地表形变反演方法
CN112050725A (zh) * 2020-09-14 2020-12-08 广东省核工业地质局测绘院 一种融合InSAR和GPS的地表形变监测方法
CN112540369A (zh) * 2020-11-27 2021-03-23 武汉大学 融合GNSS与升降轨时序InSAR的滑坡三维形变解算方法及系统
CN112540370A (zh) * 2020-12-08 2021-03-23 鞍钢集团矿业有限公司 InSAR和GNSS融合的露天矿边坡形变测量方法
CN112859077B (zh) * 2021-01-27 2023-03-07 中国测绘科学研究院 一种多级合成孔径雷达干涉相位解缠方法
CN113156456A (zh) * 2021-04-21 2021-07-23 中国公路工程咨询集团有限公司 一种路面、隧道一体化检测方法及检测设备、车辆
CN113624122B (zh) * 2021-08-10 2022-09-20 中咨数据有限公司 融合GNSS数据与InSAR技术的桥梁变形监测方法

Cited By (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115453520A (zh) * 2022-10-26 2022-12-09 中国电子科技集团公司第十四研究所 基于双频多极化差分干涉的地表形变测量方法及设备
CN115453520B (zh) * 2022-10-26 2023-01-03 中国电子科技集团公司第十四研究所 基于双频多极化差分干涉的地表形变测量方法及设备
CN115451874B (zh) * 2022-11-11 2023-02-28 深圳市科荣软件股份有限公司 一种基于InSAR图像的水坝形变监测方法及装置
CN115451874A (zh) * 2022-11-11 2022-12-09 深圳市科荣软件股份有限公司 一种基于InSAR图像的水坝形变监测方法及装置
CN116224327A (zh) * 2023-02-20 2023-06-06 中国矿业大学 一种基于学习网络的矿区大梯度形变区相位解缠方法
CN116258005A (zh) * 2023-02-23 2023-06-13 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN116258005B (zh) * 2023-02-23 2023-11-21 哈尔滨工业大学（深圳） 一种基于InSAR数据的桥梁形变预测方法及相关装置
CN116338690B (zh) * 2023-03-28 2024-01-16 中南林业科技大学 一种无模型约束的时序InSAR地形残差与地表形变估计方法
CN116338690A (zh) * 2023-03-28 2023-06-27 中南林业科技大学 一种无模型约束的时序InSAR地形残差与地表形变估计方法
CN116150301A (zh) * 2023-04-24 2023-05-23 北京中舆达信息科技有限公司 一种基于3DWebGIS的地表形变监测系统及方法
CN116736305A (zh) * 2023-08-14 2023-09-12 北京观微科技有限公司 地物形变确定方法、装置、电子设备及存储介质
CN116736305B (zh) * 2023-08-14 2023-10-27 北京观微科技有限公司 地物形变确定方法、装置、电子设备及存储介质
CN117036965A (zh) * 2023-10-08 2023-11-10 四川正路建设工程检测咨询有限公司 桥梁维修设备控制方法、电子设备和计算机可读介质
CN117036965B (zh) * 2023-10-08 2024-01-05 四川正路建设工程检测咨询有限公司 桥梁维修设备控制方法、电子设备和计算机可读介质
CN117368916A (zh) * 2023-10-12 2024-01-09 云南大学 一种InSAR相位解缠方法、系统、设备及介质
CN117213443A (zh) * 2023-11-07 2023-12-12 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法
CN117213443B (zh) * 2023-11-07 2024-03-19 江苏省地质调查研究院 一种天地深一体化地面沉降监测网建设与更新方法
CN117274342A (zh) * 2023-11-21 2023-12-22 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法
CN117274342B (zh) * 2023-11-21 2024-02-13 中铁水利水电规划设计集团有限公司 一种基于卫星数据的水利工程形变监测方法
CN117435902A (zh) * 2023-12-20 2024-01-23 武汉华威科智能技术有限公司 一种基于机器学习确定rfid标签运动行为的方法及装置
CN117435902B (zh) * 2023-12-20 2024-04-02 武汉华威科智能技术有限公司 一种基于机器学习确定rfid标签运动行为的方法及装置
CN117724102A (zh) * 2024-02-18 2024-03-19 中国特种设备检测研究院 一种结合ekfpu的mcf相位解缠方法及系统

Also Published As

Publication number Publication date
CN113624122B (zh) 2022-09-20
ZA202300146B (en) 2023-03-29
CN113624122A (zh) 2021-11-09
WO2022214114A3 (zh) 2022-11-10

Similar Documents

Publication Publication Date Title
WO2022214114A2 (zh) 融合GNSS数据与InSAR技术的桥梁变形监测方法
CN107389029B (zh) 一种基于多源监测技术融合的地面沉降集成监测方法
US20200103530A1 (en) Method for extracting elevation control point with assistance of satellite laser altimetry data
CN106526590B (zh) 一种融合多源sar影像工矿区三维地表形变监测及解算方法
US20140191894A1 (en) Three-dimensional positioning method
CN113340191B (zh) 时间序列干涉sar的形变量测量方法及sar系统
CN106526593B (zh) 基于sar严密成像模型的子像素级角反射器自动定位方法
Zhou et al. GPS/terrestrial 3D laser scanner combined monitoring technology for coal mining subsidence: a case study of a coal mining area in Hebei, China
Liu et al. Detecting subsidence in coastal areas by ultrashort-baseline TCPInSAR on the time series of high-resolution TerraSAR-X images
CN106960174A (zh) 高分影像激光雷达高程控制点提取及其辅助定位方法
CN100545676C (zh) 基于相关加权的干涉合成孔径雷达干涉相位估计方法
CN109100719B (zh) 基于星载sar影像与光学影像的地形图联合测图方法
CN103698764A (zh) 一种稀疏采样条件下的干涉合成孔径雷达成像方法
CN112051571A (zh) 一种新型差分InSAR的LOS向形变量估计方法
Cao et al. Bundle adjustment of satellite images based on an equivalent geometric sensor model with digital elevation model
CN108983231B (zh) 一种基于视频合成孔径雷达的干涉视频测量方法
Liu et al. Accurate mapping method for UAV photogrammetry without ground control points in the map projection frame
Xie et al. Analysis of deformation over permafrost regions of Qinghai-Tibet plateau based on permanent scatterers
CN102798851B (zh) 一种基于几何成像的modis lai产品验证方法
Wang et al. Layover compensation method for regional spaceborne SAR imagery without GCPs
Wang et al. A high precision DEM extraction method based on InSAR data
CN109696155A (zh) 光线共面约束的弱交会光学卫星影像联合平差方法及系统
CN114706074A (zh) 一种基于二维观测的地表三维变形InSAR测量方法2T3D-InSAR
Liu et al. Ultrashort-baseline persistent scatterer radar interferometry for subsidence detection
Falco et al. Infrastructure Stability Analysis by COSMO-SkyMed PSP SAR Interferometry: Spatio-Temporal Analysis and 3D Modeling

Legal Events

Date Code Title Description
NENP Non-entry into the national phase

Ref country code: DE