CN116449365A - Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology - Google Patents

Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology Download PDF

Info

Publication number
CN116449365A
CN116449365A CN202211456837.8A CN202211456837A CN116449365A CN 116449365 A CN116449365 A CN 116449365A CN 202211456837 A CN202211456837 A CN 202211456837A CN 116449365 A CN116449365 A CN 116449365A
Authority
CN
China
Prior art keywords
deformation
phase
image
images
dimensional
Prior art date
Legal status (The legal status 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 status listed.)
Pending
Application number
CN202211456837.8A
Other languages
Chinese (zh)
Inventor
姚建群
杨书仁
郭东华
于壮强
丁松
李超
马正源
丁力
乔炜
孟庆阳
苏新
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CCCC Infrastructure Maintenance Group Co Ltd
Original Assignee
CCCC Infrastructure Maintenance Group Co Ltd
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 CCCC Infrastructure Maintenance Group Co Ltd filed Critical CCCC Infrastructure Maintenance Group Co Ltd
Priority to CN202211456837.8A priority Critical patent/CN116449365A/en
Publication of CN116449365A publication Critical patent/CN116449365A/en
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR image post-processing techniques combined with interferometric techniques
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/38Registration of image sequences
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10032Satellite or aerial image; Remote sensing
    • G06T2207/10044Radar image

Landscapes

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

Abstract

The invention relates to a monitoring method of a subway along-line surface two-dimensional deformation field based on a time sequence InSAR technology, which comprises the following steps: acquiring SAR image data covering a monitoring range; optimizing an image pairing connection network; image registration and coherence computation; generating a differential interference pattern; image filtering and phase unwrapping; selecting GCP points; removing the atmospheric phase error; extracting a deformation field of the satellite sight line direction; calculating a true two-dimensional deformation field of the earth surface by combining monitoring results of different tracks; the method is characterized in that deformation results of subway regions along the subway line are extracted based on a subway network buffer zone, urban earth surface deformation monitoring based on a time sequence radar interferometry technology is carried out by utilizing SAR images of different satellites and different orbits covering the same region aiming at the condition that the InSAR technology can only monitor earth surface deformation in the satellite sight line direction, the deformation of the subway regions along the subway line is accurately extracted by using the subway network generation buffer zone, and data and theoretical support are provided for urban earth surface deformation monitoring and early warning of geological disasters.

Description

Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology
Technical Field
The invention relates to the technical field of synthetic aperture radar interferometry, in particular to a subway line surface two-dimensional deformation field monitoring method based on a time sequence InSAR technology.
Background
The synthetic aperture radar interferometry (InSAR, interferometric Synthetic Aperture Radar) is a large-range measurement technology based on a surface area, and has the advantages of high spatial resolution, wide coverage, low cost, rapidness and accuracy and large-scale continuous coverage capability. The method is mainly based on the technology that satellites are used for carrying out differential processing on phase data in SAR monitoring data of the same area in different periods, so that earth surface deformation information is extracted. The current means for monitoring the deformation of the ground surface along the subway mainly depend on traditional measuring means such as a level gauge, a theodolite, a total station, a GPS and the like. However, the methods need to be monitored on site by no exception, and the defects of long observation period, large consumption of manpower, material resources and financial resources, easiness in being influenced by natural conditions such as terrain, weather and the like, point measurement and the like seriously limit the general investigation and detection of the safety of the subway along the ground surface. Furthermore, these methods monitor objects based on point-based observations, rather than monitoring the entire area at once, even for a single geowire, the monitoring efficiency is low. The InSAR technology measures the tiny change of the earth surface through microwave imaging, adopts a synthetic aperture technology and a pulse compression technology in the azimuth direction and the distance direction respectively, has incomparable azimuth direction and distance direction high resolution characteristics compared with the conventional radar, and meanwhile, the active measurement mode also ensures that the radar is not interfered by weather. In other words, the InSAR technique has both high spatial resolution and high accuracy measurement scale, and with the reduction of satellite return period, the realization of high temporal resolution of the technique is also called japanese-stand. The InSAr technology is characterized by area imaging, all weather, all-day time, large range, high precision and high resolution, and is an optimal solution for solving the problem of monitoring the deformation of the ground surface along the large-range subway at present.
Under normal conditions, the InSAR technology can only acquire the surface deformation between two time points, but the time sequence InSAR technology can acquire deformation characteristics, namely the deformation time sequence of the surface in a whole period of time through joint calculation of a plurality of interference subsets, so that long-period dynamic monitoring of the surface is realized. The algorithms that are more typical in current time series InSAR technology are the small baseline set (Small Baseline Subset, SBAS) and permanent scatterer (Persistent Scatterer, PS) technologies. However, even in the above algorithm, the monitoring of the earth surface only stays in the LOS direction, and the monitoring result cannot reflect the actual deformation of the earth surface. And because the imaging angles are different, the solving results of the LOS direction deformation fields in the same area are also different by different tracks, and even the situation that the line outgoing is completely opposite can occur, the ground surface monitoring is unreliable by only using SAR data of a single track, and the misjudgment of the deformation situation can be caused.
As an important ring in urban traffic planning, subways are always an object of great interest in operation safety and monitoring. With the continuous expansion of economic zones in recent years, the scale of subways is also rapidly growing. However, the large scale extension of subways has a subtle and difficult-to-detect effect on urban surface deformation. Once the tiny deformation which is large in range and difficult to monitor is accumulated into geological disasters, serious troubles are brought to urban management, such as urban waterlogging, underground pipeline rupture and even surface collapse, and life and property security of people are threatened. Therefore, the method monitors the surface deformation of the subway along the line area, acquires the real deformation condition of the surface, evaluates the influence of the surface on surrounding buildings, and has important significance for urban planning management and people safety.
Disclosure of Invention
The invention aims to provide a subway line earth surface two-dimensional deformation field monitoring method based on a time sequence InSAR technology, which has the advantages that after SAR images of the same imaging area of two groups of different tracks are acquired, overlapping interference processing is respectively carried out on each track image set, after phases irrelevant to deformation are removed, deformation fields in LOS directions of different tracks are acquired, the results are jointly solved, the deformation fields in the vertical direction and the east-west direction of the earth surface are obtained, the real earth surface deformation condition is presented, and misjudgment caused by monitoring earth surface deformation in the LOS direction by a single track is reduced, so that the problems in the prior art are solved.
In order to achieve the above purpose, the present invention provides the following technical solutions: a subway line surface two-dimensional deformation field monitoring method based on a time sequence InSAR technology specifically comprises the following steps:
step one: acquiring SAR image data covering a monitoring range
Step two: image pairing connection network optimization
Step three: image registration and coherence computation
Step four: generating differential interferograms
Step five: image filtering and phase unwrapping
Step six: GCP point selection
Step seven: eliminating atmospheric phase error
Step eight: extracting satellite sight direction deformation field
Step nine: calculating real two-dimensional deformation field of earth surface by combining monitoring results of different tracks
Step ten: and extracting deformation results of the subway region along the subway line based on the subway network buffer area.
In the first step, taking one set of SAR images as an example, N images acquired in the same area need to be selected for subsequent interferometry, the N SAR image images are paired in a free combination form to form a connection network, but because of the fact that the images of image time and space are incoherent, if any images are freely paired, a large number of invalid low coherence interferograms are generated, and an error which cannot be estimated is brought to deformation calculation results, therefore, the condition of image pairing needs to be restrained, the connection network formed by image sets reaches the optimal value of coherence, the baseline time threshold is set to be 4 times of the interval of single images, the baseline space threshold is set to be 45% of the longest space baseline of the whole image set, the coherence of the whole connection network can be effectively provided, the success rate of image interference is improved, the reality and reliability of deformation measurement results are ensured, at this moment, the N scene images of a research area are acquired in a monitoring period, and after the threshold is screened, M differential interferograms can be obtained, and the relationship satisfies:
further, after the interference sequence between the master image and the slave image is determined, considering the offset and the stretching deformation of each SAR image pair between the master image and the slave image caused by different observation geometries, in order to realize the phase subtraction of the same name points of the master image and the slave image in the subsequent interference phase calculation process, data registration needs to be performed on each SAR image pair, and the step is completed by adopting a maximum spectrum method, and the implementation process is as follows: extracting an i×i window from the basically same area in the paired images, obtaining the interferogram of the group of images by common-short multiplication, carrying out two-dimensional fast Fourier transform on the interferogram to obtain a spectrogram of interference data, calculating the signal-to-noise ratio of the frequency spectrum, fixing a source image data block, shifting a target image data block by one pixel in azimuth or distance direction, repeating the steps, scanning n pixels in azimuth and distance direction, obtaining a (n-1) x (n-1) frequency spectrum maximum matrix, finding the extremum of the frequency spectrum maximum matrix, wherein the optimal offset of the images is located at the extremum, coherence is an important basis for evaluating the echo similarity degree between two images and the phase quality of the interferogram, the coherence between the images largely determines the effect of interference processing, and the coherence is also required to be used as a processing of a filtering threshold when the interferogram is subjected to subsequent filtering processing, so that after the images are registered, the coherence calculation is required for the paired images, the value of the coherence can be calculated by determining a sliding window, and the SAR (synthetic aperture) is defined as follows:
wherein, c 1 And c 2 Complex value of SAR image, c 2 * Is c 2 Complex conjugate of phi det The method is an interference phase value superimposed by various factors such as baseline errors, topography fluctuation, surface deformation and the like in a coherence estimation window, and the coherence value of each pixel can be calculated through a preset moving window to obtain a coherence map of an interference pair.
Preferably, in the fourth step, the phase map in the paired images is subjected to differential processing, and the differential phase map includes deformation phase changes of the earth surface of the images imaged twice, and a specific differential formula is as follows:
δφ x =φ x (t b )-φ x (t a )=φ flat,xele,xdisp,xatm,xnoise
wherein phi is x (t b ) And phi x (t a ) At t a And t b The images acquired by the radar images at two different moments in time in the area, x represents the value of a pixel in the image, delta phi x The phase difference obtained after the interference of the two images. Phi is flat,x 、φ ele,x 、φ disp,x 、φ atm,x 、φ noise The InSAR technology is mainly aimed at inverting the surface deformation value by extracting deformation phases in the interferogram, so that after the paired image interferogram is obtained, the phase irrelevant to deformation is needed to be removed, and the deformation phases are extracted, which specifically comprises the step of introducing external DEM eliminates the ground elevation phase (phi) ele,x ) The method comprises the steps of carrying out a first treatment on the surface of the Cancellation of phase (phi) generated by ground leveling effect by track data flat,x ) The method comprises the steps of carrying out a first treatment on the surface of the Attenuation of measured noise corresponding phase (phi) by image Goldstein filtering algorithm noise ) Is a function of (a) and (b). The Goldstein filtering algorithm mainly uses Fourier transformation to convert an interferogram from a space domain to a frequency domain, divides the interferogram into overlapped sliding windows, and then performs smoothing processing on a power spectrum of each window, wherein a filtering formula is as follows:
H(u,v)=S{|Z(u,v)|} α *Z(u,v)
wherein a is a filter factor; s is a smoothing operator; z (u, v) is the image to be processed; h (uv) is the filtered image; u and v are spatial frequencies. When the InSAR interferogram is generated, the interference phase is in a winding state and the phase value is between [ -pi, pi), so that the real phase information can be obtained by disentangling the phase after filtering, under ideal conditions, the relative phase between any point and a certain principal value can be obtained through simple integration, but under actual conditions, the phase disentanglement becomes extremely complex due to the difference of complex ground fluctuation and the data quality of the interference image and the influence of noise and the like, the minimum cost flow algorithm is an algorithm with stronger adaptability and higher precision in a plurality of disentanglement algorithms, the main idea is to minimize the difference between the winding phase and the disentanglement phase, then convert the minimization problem into the minimum cost flow problem, greatly reduce the difficulty of knowing the winding algorithm from the network optimization angle, and although the ground surface elevation phase and the ground surface phase are removed to a certain extent in the previous step, a part of the phase can not be removed, therefore, the ground control point (Ground Control Point, GCP) is required to be selected, the residual phase is subjected to interpretation, and the phase is completely removed, the effect of the phase is very important in the selection of the phase, the phase is not required to be manually selected, and the phase is not subjected to the distortion control in the conventional method.
Preferably, in the seventh step, since the particles in the atmosphere scatter radar waves to a certain extent, and the atmospheric delay effects caused by different environments and different times are different, in the differential interference technology, the atmospheric effect is always an error source that affects the accuracy quite but is not easy to remove, although the atmospheric effect can be eliminated to a certain extent by using an image-to-image superposition method, the atmospheric effect error can only be dispersed in opposition by the image-to-image superposition method, but the atmospheric effect cannot be estimated and eliminated truly, and the current main method for removing the atmospheric phase is to use a time-space filter, which assumes that the atmospheric phase and the track error phase are spatially related in a certain area range, and uses low-pass filtering on the time domain, so that the atmospheric effect error of the main image and each sub-image can be effectively estimated, thereby providing the atmospheric effect error in the real phase.
Preferably, in the step eight, since the baseline threshold is set in the step two, the SAR image is divided into a plurality of independent data subsets, and the rank of the coefficient matrix is not the full rank, so that the deformation solution equation of the time sequence InSAR is rank deficient, and the solution of the equation set is not unique, for this purpose, singular value decomposition is performed on the coefficient matrix, the equation reaches the state of full rank by adding constraint, and the least square solution under the meaning of the minimum norm is obtained, and then the surface deformation is obtained, and the specific process is as follows: after eliminating the phase irrelevant to deformation in the fifth step and the seventh step, the information represented by the phase in the interferogram is:
δφ x =φ disp,xnoise
the expression of the surface deformation phase and deformation rate is as follows:
in the formula, lambda is the radar wavelength, V i For deformation rate, and meanwhile, in order to eliminate errors caused by measurement noise, a least square method is needed to carry out joint solution on a plurality of groups of interference pairs, and all interference pairs are combined into a matrix form, which can be expressed as:
δφ=BV
B MxN the method is characterized in that the method is a coefficient matrix, a minimum norm solution of a velocity vector can be obtained by utilizing a singular value decomposition method, and finally, the time of each velocity is integrated to obtain an accumulated deformation value in a monitoring period, although the time sequence InSAR technology can monitor the surface deformation of millimeter level, the deformation value result is not in a vertical direction, but is a one-dimensional projection of a radar Sight Line direction (LOS) in the LOS direction, wherein the deformation value is actually in a three-dimensional direction (east direction, north direction and vertical direction), in other words, the result is composed of the comprehensive contribution of the three-dimensional deformation, and the result is not the real deformation of the surface. To calculate the real earth surface three-dimensional deformation field of the monitoring area, the LOS deformation field after the SAR data interference of the lifting rail is needed to be used, a plurality of equations are constructed through the imaging principle of the lifting rail to solve the eastern deformation (D E ) Deformation in the north direction (D N ) Deformation in the vertical direction (D U ) These three unknowns. For the track lifting image, the three-dimensional deformation of the earth surface and the deformation in the LOS direction are as follows:
wherein the method comprises the steps ofIs deformation in LOS direction, +.>For the angle of incidence of an orbiting satellite->For azimuth angles of the orbiting satellites, the imaging geometry equation for the same thing as the orbit is:
since the sight direction of the lifting rail is basically consistent with the east-west direction, the lifting rail is insensitive to the deformation value in the north-south direction, in other words, the lifting rail is formedUnknown quantity D in image geometry equation N The corresponding coefficient is very small, we consider D N The contribution to the LOS directional deformation is negligible, so D in the formula can be used in inverting the two-dimensional deformation field N The corresponding terms are deleted, and the formula for solving the two-dimensional deformation field by the lifting rail can be written into a matrix form:
L=BX
wherein the method comprises the steps ofX=[D U D E ] T The method comprises the steps of carrying out a first treatment on the surface of the B is a coefficient array:
the calculation unit for inversion of the earth surface two-dimensional deformation by the lifting rail time sequence InSAR is a single pixel point, the two deformation graphs can be converted into the scale of the same resolution and the same matrix size after resampling, at the moment, the pixels of the same row and column numbers of the two deformation graphs represent the same ground object, and finally, the two-dimensional deformation value of each pixel point is calculated by using a matrix operation cycle, so that the deformation graph in the east-west and vertical directions can be obtained.
Preferably, in the step ten, after the deformation fields in the vertical and east-west directions of the monitored area are obtained, the subway net vector data is used to generate a buffer area of 500m, and the buffer area element is used as a clipping area to clip the two-dimensional velocity field and the accumulated two-dimensional deformation field of the SAR data, so that the deformation monitoring result only aiming at the area along the subway is extracted.
Compared with the prior art, the invention has the following beneficial effects:
according to the subway line surface two-dimensional deformation field monitoring method based on the time sequence InSAR technology, after SAR images of the same imaging area of two groups of different tracks are acquired, overlapping interference processing is carried out on each track image set, after phases irrelevant to deformation are removed, deformation fields in LOS directions of different tracks are acquired, the results are combined and calculated, deformation fields in the vertical direction and the east-west direction of the surface are obtained, a real surface deformation condition is presented, and misjudgment caused by single track monitoring LOS direction surface deformation is reduced.
Drawings
FIG. 1 is a process flow diagram of the present invention;
FIG. 2 is a graph of the deformation rate of the elevated rail surface in the LOS direction of the present invention;
FIG. 3 is a graph of the deformation rate of the derailment earth surface in the LOS direction of the present invention;
FIG. 4 is a plot of the deformation rate of the earth's surface in the vertical direction according to the present invention;
FIG. 5 is an illustration of the east-west deformation rate of the earth's surface in accordance with the present invention;
FIG. 6 is a graph of deformation rate in the vertical direction of a subway according to the invention;
fig. 7 is a graph of sampling point LOS direction and two-dimensional direction deformation rate according to the present invention.
Detailed Description
The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention, and it is apparent that the described embodiments are only some embodiments of the present invention, but not all embodiments, and all other embodiments obtained by those skilled in the art without making creative efforts based on the embodiments of the present invention are included in the protection scope of the present invention.
In the description of the present invention, it should be understood that the terms "upper," "lower," "front," "rear," "left," "right," "top," "bottom," "inner," "outer," and the like indicate or are based on the orientation or positional relationship shown in the drawings, merely to facilitate description of the present invention and to simplify the description, and do not indicate or imply that the devices or elements referred to must have a specific orientation, be configured and operated in a specific orientation, and thus should not be construed as limiting the present invention.
Example 1:
referring to fig. 1-7, a method for monitoring a two-dimensional deformation field of a subway along a line surface based on a time sequence InSAR technology comprises the following steps: acquiring SAR image data covering a monitoring range
The method comprises the following steps:
the adult is one of the most rapidly developed cities in China in recent years, and as a city of the world known by the tourism industry and the catering industry, the adult government has carried out large-scale extension on urban areas and traffic routes in order to increase the population capacity of the city. By 2021, 1 month, 12 subway lines of the capital have been built, and the total mileage reaches 518.96 km. Aiming at cities with a large number of high-rise building groups and complex underground engineering structures, the method is indispensable to large-scale and high-frequency settlement observation of subways along the lines. The experiment utilizes 37 Jing Shenggui and a C-band sentinel satellite for monitoring a city district of a capital center, the monitoring period is from 2017 month 3 to 2021 month 3, and the time sequence InSAR technology is adopted to acquire the earth surface deformation data of a research area from SAR images. The basic information of the two sets of data is shown in Table 1.
TABLE 1
Step two: image pairing connection network optimization
The method comprises the following steps:
the 37 SAR image images of the ascending rail and the descending rail are paired in a free combination mode, so that a connecting network is formed. However, since images of which the time and space are incoherent are freely paired, a large number of invalid low-coherence interferograms are generated and an error which cannot be estimated is brought to a deformation calculation result, the conditions of image pairing are required to be restrained, and a connection network formed by image sets reaches an optimal value of coherence. In the experiment, the baseline time threshold is set to 150 days, the baseline space threshold is set to 45% of the longest space baseline of the whole image set, and the image with the smallest combination of three parameters of the time baseline, the space baseline and the Doppler centroid frequency difference is comprehensively considered as the main image. The common main image of the track lifting is an image of 2018, 09 and 11 months, and the common main image of the track lowering is an image of 2018, 10 and 12 months.
Step three: image registration and coherence computation
The method comprises the following steps:
after the interference sequence between the main image and the auxiliary image is determined, considering the offset and the stretching deformation of each SAR image pair between the main image and the auxiliary image caused by different observation geometries, in order to realize the phase subtraction of the same-name points of the main image and the auxiliary image in the subsequent interference phase calculation process, data registration needs to be carried out on each SAR image pair. This step will be done using maximum spectroscopy. The implementation flow is as follows: and respectively extracting an i multiplied by i window from the basically same region in the paired images, obtaining an interference pattern of the group of images through common multiplication, carrying out two-dimensional fast Fourier transform on the interference pattern, obtaining a spectrogram of interference data, and calculating the signal to noise ratio of the frequency spectrum. The source image data block is fixed and the target image data block is shifted by one pixel in azimuth or distance, repeating the above steps. Thus, n pixels are scanned upwards in azimuth and distance, a (n-1) x (n-1) spectrum maximum matrix is obtained, the extreme value of the spectrum maximum matrix can be found, and the optimal offset of the image is located at the position of the extreme value. Coherence is an important basis for evaluating the degree of echo similarity and the phase quality of interferograms between two images. The coherence between images largely determines the effect of the interference process, and the correlation is also required as a filter threshold when the interference pattern is subjected to subsequent filter processes, so that the coherence calculation is required for the paired images after image registration is performed. The value of the coherence can be obtained by determining a sliding window and calculating the cross-correlation of two SAR images in the window, and the definition formula is as follows:
wherein, c 1 And c 2 Complex value of SAR image, c 2 * Is c 2 Complex conjugate of phi det The method is an interference phase value superimposed by various factors such as baseline errors, topography fluctuation, surface deformation and the like in a coherence estimation window, and the coherence value of each pixel can be calculated through a preset moving window to obtain a coherence map of an interference pair.
Step four: generating differential interferograms
The method comprises the following steps:
after the image pairing is completed, carrying out differential processing on the paired phase diagram, wherein the differential phase diagram comprises deformation phase change of the earth surface of the image formed by two times of imaging, and the specific differential formula is as follows:
δφ x =φ x (tb)-φ x (t a )=φ flat,xele,xdiap,xatm,xnoise
wherein phi is x (t b ) And phi x (t a ) At t a And t b The images acquired by the radar images at two different moments in time in the area, x represents the value of a pixel in the image, delta phi x The phase difference obtained after the interference of the two images. Phi is flat,x 、φ ele,x 、φ disp,x 、φ atm,x 、φ noise The phase caused by the earth effect, the phase caused by the earth surface elevation, the phase caused by the earth surface deformation, the phase caused by the atmospheric effect and the phase caused by the measurement noise are respectively shown.
Step five: image filtering and phase unwrapping
The method comprises the following steps:
the InSAR technology is mainly used for inverting the surface deformation value through extracting the deformation phase in the interferogram, so that after the paired image interferogram is obtained, the phase irrelevant to deformation is needed to be removed, and the deformation phase is extracted. Specifically includes introducing an external DEM to eliminate the ground surface elevation phase (phi) ele,x ) The method comprises the steps of carrying out a first treatment on the surface of the Cancellation of phase (phi) generated by ground leveling effect by track data flat,x ) The method comprises the steps of carrying out a first treatment on the surface of the Attenuation of measured noise corresponding phase (phi) by image Goldstein filtering algorithm noise ) Is a function of (a) and (b). The Goldstein filtering algorithm mainly uses Fourier transformation to convert an interferogram from a space domain to a frequency domain, divides the interferogram into overlapped sliding windows, and then performs smoothing processing on a power spectrum of each window, wherein a filtering formula is as follows:
H(u,v)=S{|Z(u,v)|} α *Z(u,v)
wherein a is a filter factor; s is a smoothing operator; z (u, v) is the image to be processed; h (uv) is the filtered image; u and v are spatial frequencies. When the InSAR interferogram is generated, the interference phase is in a winding state, and the phase value is between [ -pi, pi), so that the actual phase information can be obtained by disentangling the phase after filtering, and the relative phase between any point and a certain principal value can be obtained through simple integration under ideal conditions, but under actual conditions, the disentanglement of the phase is complicated due to the difference of complex ground fluctuation and the data quality of the interference image and the influence of noise and the like. The minimum cost flow algorithm is an algorithm with stronger adaptability and higher precision in a plurality of unwrapping algorithms. The main idea is to minimize the difference between winding phase and unwinding phase, and then to convert the minimization problem into the minimum cost flow problem, so as to greatly reduce the difficulty of knowing the winding algorithm from the network optimization perspective.
Step six: GCP point selection
The method comprises the following steps:
although the ground surface elevation phase and the ground phase are removed to a certain extent in the previous step, a part of the phases still remain and cannot be removed, so GCP is selected, and the residual phases are fitted, so that the phases are completely removed. GCP is very important in the selection of InSAR, where control points are selected generally requiring locations without residual terrain phase, deformation phase and phase jumps. The experiment is mainly based on visual interpretation and a coherence map to manually select GCP.
Step seven: eliminating atmospheric phase error
The method comprises the following steps:
because particles in the atmosphere can scatter radar waves to a certain extent, and the atmospheric delay effects caused by different environments and different times are different, in the differential interference technology, the atmospheric effect is always an error source which has quite large influence on accuracy but is not easy to remove. Although the method of overlapping the images can eliminate the atmospheric effect to a certain extent, the processing of overlapping the images can only disperse the atmospheric effect error evenly on the opposite of each image, and cannot actually estimate and eliminate the atmospheric effect. The experiment was rejected using a time-space filter. The method assumes that the atmospheric phase and the orbit error phase are spatially related in a certain region range, carries out high-pass filtering on the unwrapping phase in a time domain, and uses low-pass filtering on the space domain, so that the atmospheric effects of the main image and each auxiliary image can be effectively estimated, and the atmospheric effect errors are proposed in the real phase.
Step eight: extracting satellite sight direction deformation field
The method comprises the following steps:
because the baseline threshold is set in the second step, the SAR image is divided into a plurality of independent data subsets, the rank of the coefficient matrix is not full rank, therefore, the deformation solution equation of the time sequence InSAR is rank deficient, and the solution of the equation set is not unique. Therefore, singular value decomposition is carried out on the coefficient array, constraint is added to enable the equation to reach a state of full rank, a least square solution in the meaning of the minimum norm is obtained, and then the surface deformation is obtained. The specific process is as follows: after eliminating the phase irrelevant to deformation in the fifth step and the seventh step, the information represented by the phase in the interferogram is:
δφ x =φ disp,xnoise
the expression of the surface deformation phase and deformation rate is as follows:
in the formula, lambda is the radar wavelength, V i For deformation rate, and meanwhile, in order to eliminate errors caused by measurement noise, a least square method is needed to carry out joint solution on a plurality of groups of interference pairs, and all interference pairs are combined into a matrix form, which can be expressed as:
δφ=BV
B MxN at this time, a minimum norm solution of the velocity vector can be obtained by using a singular value decomposition method, and finally, the time of each velocity is integrated to obtain the accumulated deformation value in the monitoring period.
Step nine: calculating real two-dimensional deformation field of earth surface by combining monitoring results of different tracks
The method comprises the following steps:
although the time series InSAR technique can monitor the surface deformation in millimeter scale, the deformation value result is not in the vertical direction, but in the Line of Sight (LOS) direction of the radar. The deformation value is actually a one-dimensional projection of the three-dimensional direction (east direction, north direction, vertical direction) in the LOS direction. In other words, the result consists of a comprehensive contribution of the three-dimensional directional deformation, which results exhibit a true deformation not of the earth's surface. To calculate the real earth surface three-dimensional deformation field of the monitoring area, the LOS deformation field after the SAR data interference of the lifting rail is needed to be used, a plurality of equations are constructed through the imaging principle of the lifting rail to solve the eastern deformation (D E ) Deformation in the north direction (D N ) Deformation in the vertical direction (D U ) These three unknowns. For the track lifting image, the three-dimensional deformation of the earth surface and the deformation in the LOS direction are as follows:
(6)
wherein the method comprises the steps ofIs deformation in LOS direction, +.>For the angle of incidence of an orbiting satellite->For azimuth angles of the orbiting satellites, the imaging geometry equation for the same thing as the orbit is:
since the sight direction of the lifting rail is basically consistent with the east-west direction, the lifting rail is insensitive to the deformation value in the north-south direction, in other words, the unknown quantity D in the imaging geometric equation N The corresponding coefficient is very small, we consider D N For LOS directionThe deformation contribution is negligible, so D in the formula can be used in inverting the two-dimensional deformation field N The corresponding terms are deleted, and the formula for solving the two-dimensional deformation field by the lifting rail can be written into a matrix form:
L=BX
wherein the method comprises the steps ofX=[D U D E ] T The method comprises the steps of carrying out a first treatment on the surface of the B is a coefficient array:
the calculation unit for inversion of the earth surface two-dimensional deformation by the lifting rail time sequence InSAR is a single pixel point, the two deformation graphs can be converted into the scale of the same resolution and the same matrix size after resampling, and at the moment, the pixels of the same row and column numbers of the two deformation graphs represent the same ground object. And finally, circularly calculating the two-dimensional deformation value of each pixel point by using matrix operation, and obtaining deformation graphs in the east-west and vertical directions.
Step ten: extracting deformation results of subway line area based on subway network buffer area
The method comprises the following steps:
after the deformation fields in the vertical and east-west directions of the monitoring area are obtained, the subway net vector data is utilized to generate a buffer area of 500m, the buffer area elements are utilized as the cutting area, and the two-dimensional velocity field and the accumulated two-dimensional deformation field of SAR data are cut, so that the deformation monitoring result only aiming at the subway area along the subway line is extracted.
While the fundamental and principal features of the invention and advantages of the invention have been shown and described, it will be apparent to those skilled in the art that the invention is not limited to the details of the foregoing exemplary embodiments, but may be embodied in other specific forms without departing from the spirit or essential characteristics thereof; the present embodiments are therefore to be considered in all respects as illustrative and not restrictive, the scope of the invention being indicated by the appended claims rather than by the foregoing description, and all changes which come within the meaning and range of equivalency of the claims are therefore intended to be embraced therein.
Furthermore, it should be understood that although the present disclosure describes embodiments, not every embodiment is provided with a separate embodiment, and that this description is provided for clarity only, and that the disclosure is not limited to the embodiments described in detail below, and that the embodiments described in the examples may be combined as appropriate to form other embodiments that will be apparent to those skilled in the art.

Claims (6)

1. A subway line surface two-dimensional deformation field monitoring method based on a time sequence InSAR technology is characterized by comprising the following steps of: the method specifically comprises the following steps:
step one: acquiring SAR image data covering a monitoring range
Step two: image pairing connection network optimization
Step three: image registration and coherence computation
Step four: generating differential interferograms
Step five: image filtering and phase unwrapping
Step six: GCP point selection
Step seven: eliminating atmospheric phase error
Step eight: extracting satellite sight direction deformation field
Step nine: calculating real two-dimensional deformation field of earth surface by combining monitoring results of different tracks
Step ten: and extracting deformation results of the subway region along the subway line based on the subway network buffer area.
2. The subway line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology according to claim 1, wherein the method is characterized by comprising the following steps: taking one set of SAR images as an example, N images acquired in the same area need to be selected for subsequent interferometry, the N SAR image images are paired in a free combination form to form a connection network, but if any images are freely paired, a large number of invalid low-coherence interferograms are generated, which brings an estimation error to deformation calculation results, so that conditions of image pairing need to be restrained, the connection network formed by image sets reaches an optimal value of coherence, a baseline time threshold is set to be 4 times of an interval of a single image, a baseline space threshold is set to be 45% of a longest spatial baseline of the whole image set, the coherence of the whole connection network can be effectively provided, the success rate of image interference is improved, the reality and reliability of deformation measurement results are ensured, at this time, N scene images of a research area are acquired in a monitoring period, and after threshold screening, M differential interferograms can be obtained, and the relationship satisfies:
further, after the interference sequence between the master image and the slave image is determined, considering the offset and the stretching deformation of each SAR image pair between the master image and the slave image caused by different observation geometries, in order to realize the phase subtraction of the same name points of the master image and the slave image in the subsequent interference phase calculation process, data registration needs to be performed on each SAR image pair, and the step is completed by adopting a maximum spectrum method, and the implementation process is as follows: extracting an i×i window from the basically same area in the paired images, obtaining the interferogram of the group of images by common-short multiplication, carrying out two-dimensional fast Fourier transform on the interferogram to obtain a spectrogram of interference data, calculating the signal-to-noise ratio of the frequency spectrum, fixing a source image data block, shifting a target image data block by one pixel in azimuth or distance direction, repeating the steps, scanning n pixels in azimuth and distance direction, obtaining a (n-1) x (n-1) frequency spectrum maximum matrix, finding the extremum of the frequency spectrum maximum matrix, wherein the optimal offset of the images is located at the extremum, coherence is an important basis for evaluating the echo similarity degree between two images and the phase quality of the interferogram, the coherence between the images largely determines the effect of interference processing, and the coherence is also required to be used as a processing of a filtering threshold when the interferogram is subjected to subsequent filtering processing, so that after the images are registered, the coherence calculation is required for the paired images, the value of the coherence can be calculated by determining a sliding window, and the SAR (synthetic aperture) is defined as follows:
wherein, c 1 And c 2 Complex value of SAR image, c 2 * Is c 2 Complex conjugate of phi det The method is an interference phase value superimposed by various factors such as baseline errors, topography fluctuation, surface deformation and the like in a coherence estimation window, and the coherence value of each pixel can be calculated through a preset moving window to obtain a coherence map of an interference pair.
3. The subway line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology according to claim 1, wherein the method is characterized by comprising the following steps: in the fourth step, the phase map in the paired images is subjected to differential processing, and the differential phase map contains deformation phase changes of the earth surface of the images imaged twice, wherein a specific differential formula is as follows:
δφ x =φ x (t b )-φ x (t a )=φ flat,xele,xdisp,xatm,xnoise
wherein phi is x (t b ) And phi x (t a ) At t a And t b The images acquired by the radar images at two different moments in time in the area, x represents the value of a pixel in the image, delta phi x The phase difference obtained after the interference of the two images. Phi is flat,x 、φ ele,x 、φ disp,x 、φ atm,x 、φ noise Respectively representing the phase generated by the earth effect, the phase generated by the earth surface elevation, the phase generated by the earth surface deformation, the phase generated by the atmospheric effect and the phase generated by the measurement noise, inSThe main purpose of the AR technique is to invert the surface deformation value by extracting the deformation phase in the interferogram, so that after the paired image interferogram is obtained, the phase irrelevant to deformation needs to be removed, thereby extracting the deformation phase, specifically including introducing an external DEM to eliminate the surface elevation phase (phi) ele,x ) The method comprises the steps of carrying out a first treatment on the surface of the Cancellation of phase (phi) generated by ground leveling effect by track data flat,x ) The method comprises the steps of carrying out a first treatment on the surface of the Attenuation of measured noise corresponding phase (phi) by image Goldstein filtering algorithm noise ) Is a function of (a) and (b). The Goldstein filtering algorithm mainly uses Fourier transformation to convert an interferogram from a space domain to a frequency domain, divides the interferogram into overlapped sliding windows, and then performs smoothing processing on a power spectrum of each window, wherein a filtering formula is as follows:
H(u,v)=S{|Z(u,v)|} α *Z(u,v)
wherein a is a filter factor; s is a smoothing operator; z (u, v) is the image to be processed; h (uv) is the filtered image; u and v are spatial frequencies. When the InSAR interferogram is generated, the interference phase is in a winding state and the phase value is between [ -pi, pi), so that the real phase information can be obtained by disentangling the phase after filtering, under ideal conditions, the relative phase between any point and a certain principal value can be obtained through simple integration, but under actual conditions, the phase disentanglement becomes extremely complex due to the difference of complex ground fluctuation and the data quality of the interference image and the influence of noise and the like, the minimum cost flow algorithm is an algorithm with stronger adaptability and higher precision in a plurality of disentanglement algorithms, the main idea is to minimize the difference between the winding phase and the disentanglement phase, then convert the minimization problem into the minimum cost flow problem, greatly reduce the difficulty of knowing the winding algorithm from the network optimization angle, and although the ground surface elevation phase and the ground surface phase are removed to a certain extent in the previous step, a part of the phase can not be removed, therefore, the ground control point (Ground Control Point, GCP) is required to be selected, the residual phase is subjected to interpretation, and the phase is completely removed, the effect of the phase is very important in the selection of the phase, the phase is not required to be manually selected, and the phase is not subjected to the distortion control in the conventional method.
4. The subway line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology according to claim 1, wherein the method is characterized by comprising the following steps: in the seventh step, since the particles in the atmosphere scatter radar waves to a certain extent, and the atmospheric delay effects caused by different environments and different times are different, in the differential interference technology, the atmospheric effect is always an error source with quite large influence precision but not easy to remove, although the atmospheric effect can be eliminated to a certain extent by using an image-to-image superposition method, the atmospheric effect error can only be averagely dispersed in each image by the image-to-image superposition method, but the atmospheric effect cannot be really estimated and eliminated, and the mainstream method for removing the atmospheric phase is to use a time-space filter.
5. The subway line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology according to claim 1, wherein the method is characterized by comprising the following steps: in the eighth step, since the baseline threshold is set in the second step, the SAR image is divided into a plurality of independent data subsets, the rank of the coefficient matrix is not the full rank, so the deformation solution equation of the time sequence InSAR is rank deficient, the solution of the equation set is not unique, for this purpose, singular value decomposition is performed on the coefficient matrix, the equation is brought to the state of full rank by adding constraint, the least square solution in the sense of the minimum norm is obtained, and the surface deformation is obtained, and the specific process is as follows: after eliminating the phase irrelevant to deformation in the fifth step and the seventh step, the information represented by the phase in the interferogram is:
δφ x =φ disp,xnoise
the expression of the surface deformation phase and deformation rate is as follows:
in the formula, lambda is the radar wavelength, V i For deformation rate, and meanwhile, in order to eliminate errors caused by measurement noise, a least square method is needed to carry out joint solution on a plurality of groups of interference pairs, and all interference pairs are combined into a matrix form, which can be expressed as:
δφ=BV
B MxN the method is characterized in that the method is a coefficient matrix, a minimum norm solution of a velocity vector can be obtained by utilizing a singular value decomposition method, and finally, the time of each velocity is integrated to obtain an accumulated deformation value in a monitoring period, although the time sequence InSAR technology can monitor the surface deformation of millimeter level, the deformation value result is not in a vertical direction, but is a one-dimensional projection of a radar Sight Line direction (LOS) in the LOS direction, wherein the deformation value is actually in a three-dimensional direction (east direction, north direction and vertical direction), in other words, the result is composed of the comprehensive contribution of the three-dimensional deformation, and the result is not the real deformation of the surface. To calculate the real earth surface three-dimensional deformation field of the monitoring area, the LOS deformation field after the SAR data interference of the lifting rail is needed to be used, a plurality of equations are constructed through the imaging principle of the lifting rail to solve the eastern deformation (D E ) Deformation in the north direction (D N ) Deformation in the vertical direction (D U ) These three unknowns. For the track lifting image, the three-dimensional deformation of the earth surface and the deformation in the LOS direction are as follows:
wherein the method comprises the steps ofIs deformation in LOS direction, +.>For the angle of incidence of an orbiting satellite->For azimuth angles of the orbiting satellites, the imaging geometry equation for the same thing as the orbit is:
since the sight direction of the lifting rail is basically consistent with the east-west direction, the lifting rail is insensitive to the deformation value in the north-south direction, in other words, the unknown quantity D in the imaging geometric equation N The corresponding coefficient is very small, we consider D N The contribution to the LOS directional deformation is negligible, so D in the formula can be used in inverting the two-dimensional deformation field N The corresponding terms are deleted, and the formula for solving the two-dimensional deformation field by the lifting rail can be written into a matrix form:
L=BX
wherein the method comprises the steps ofX=[D U D E ] T The method comprises the steps of carrying out a first treatment on the surface of the B is a coefficient array:
the calculation unit for inversion of the earth surface two-dimensional deformation by the lifting rail time sequence InSAR is a single pixel point, the two deformation graphs can be converted into the scale of the same resolution and the same matrix size after resampling, at the moment, the pixels of the same row and column numbers of the two deformation graphs represent the same ground object, and finally, the two-dimensional deformation value of each pixel point is calculated by using a matrix operation cycle, so that the deformation graph in the east-west and vertical directions can be obtained.
6. The subway line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology according to claim 1, wherein the method is characterized by comprising the following steps: and step ten, after the deformation fields in the vertical and east-west directions of the monitoring area are obtained, generating a buffer area of 500m by utilizing subway network vector data, and utilizing the buffer area element as a cutting area to cut the two-dimensional velocity field and the accumulated two-dimensional deformation field of SAR data, thereby extracting the deformation monitoring result only aiming at the area along the subway.
CN202211456837.8A 2022-11-21 2022-11-21 Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology Pending CN116449365A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211456837.8A CN116449365A (en) 2022-11-21 2022-11-21 Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211456837.8A CN116449365A (en) 2022-11-21 2022-11-21 Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology

Publications (1)

Publication Number Publication Date
CN116449365A true CN116449365A (en) 2023-07-18

Family

ID=87129051

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211456837.8A Pending CN116449365A (en) 2022-11-21 2022-11-21 Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology

Country Status (1)

Country Link
CN (1) CN116449365A (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117214898A (en) * 2023-11-09 2023-12-12 中国电子科技集团公司第十四研究所 Wide-scale earth surface deformation refined remote sensing method based on multidimensional electromagnetic information depth fusion
CN117826148A (en) * 2023-11-29 2024-04-05 北京市市政工程研究院 Method and system for identifying coherent point
CN118089611A (en) * 2024-04-17 2024-05-28 东南大学 Building three-way displacement monitoring method and system integrating InSAR data and physical knowledge

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117214898A (en) * 2023-11-09 2023-12-12 中国电子科技集团公司第十四研究所 Wide-scale earth surface deformation refined remote sensing method based on multidimensional electromagnetic information depth fusion
CN117214898B (en) * 2023-11-09 2024-01-23 中国电子科技集团公司第十四研究所 Wide-scale earth surface deformation refined remote sensing method based on multidimensional electromagnetic information depth fusion
CN117826148A (en) * 2023-11-29 2024-04-05 北京市市政工程研究院 Method and system for identifying coherent point
CN118089611A (en) * 2024-04-17 2024-05-28 东南大学 Building three-way displacement monitoring method and system integrating InSAR data and physical knowledge

Similar Documents

Publication Publication Date Title
Colesanti et al. Monitoring landslides and tectonic motions with the Permanent Scatterers Technique
CN116449365A (en) Metro along-line surface two-dimensional deformation field monitoring method based on time sequence InSAR technology
Colesanti et al. SAR monitoring of progressive and seasonal ground deformation using the permanent scatterers technique
Crosetto et al. Persistent scatterer interferometry
Noferini et al. Using GB-SAR technique to monitor slow moving landslide
Price et al. Small‐scale deformations associated with the 1992 Landers, California, earthquake mapped by synthetic aperture radar interferometry phase gradients
Guéguen et al. Monitoring residual mining subsidence of Nord/Pas-de-Calais coal basin from differential and Persistent Scatterer Interferometry (Northern France)
CN108957456A (en) Landslide monitoring and EARLY RECOGNITION method based on multi-data source SBAS technology
Guo et al. Land subsidence in Tianjin for 2015 to 2016 revealed by the analysis of Sentinel-1A with SBAS-InSAR
Yao et al. Research on Surface Deformation of Ordos Coal Mining Area by Integrating Multitemporal D‐InSAR and Offset Tracking Technology
CN107064933A (en) The method that SAR based on circulation Power estimation chromatographs depth of building
Hu et al. Time-series InSAR technology for ascending and descending orbital images to monitor surface deformation of the metro network in Chengdu
CN113238228B (en) Three-dimensional earth surface deformation obtaining method, system and device based on level constraint
CN114091274A (en) Landslide susceptibility evaluation method and system
Yang et al. Monitoring urban subsidence with multi-master radar interferometry based on coherent targets
CN117541929A (en) Deformation risk assessment method for large-area power transmission channel of InSAR in complex environment
Berardino et al. Small baseline DIFSAR techniques for earth surface deformation analysis
Vadon et al. Earthquake displacement fields mapped by very precise correlation. Complementarity with radar interferometry
Azeriansyah et al. Land Subsidence Monitoring in Semarang and Demak Coastal Areas 2016-2017 Using Persistent Scatterer Interferometric Synthetic Aperture Radar
Hammad et al. Landslide investigation using differential synthetic aperture radar interferometry: a case study of Balloran dam area in Syria
Haque et al. Time series analysis of subsidence in Dhaka City, Bangladesh using Insar
Zhao et al. Different scale land subsidence and ground fissure monitoring with multiple InSAR techniques over Fenwei basin, China
Suhadha et al. Precise coseismic displacement related to the 2022 pasaman earthquake using multi-geometry of Sentinel-1 InSAR
Atanasova et al. Co-seismic surface displacements after the earthquakes in Larissa, 3 March 2021, derived by DInSAR
Dimitrov et al. Monitoring of the landslide processes at the" Dalgiya Yar" landslide

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination