Method for separating InSAR (Interferometric Synthetic Aperture Radar) earth surface tiny deformation phase

Abstract

The invention discloses a method for separating InSAR (Interferometric Synthetic Aperture Radar) earth surface tiny deformation phase, which comprises the following steps of: obtaining a radar remote sensing image related to a target area, and registering image data by combining precise orbit data and ground digital elevation model data of the target area to obtain full-time-domain image information; an InSAR differential interferogram is extracted according to the full-time-domain image information, then phase feature extraction and optimization processing are carried out on the differential interferogram, so that a regional deformation graph is solved, and each pixel in the regional deformation graph comprises deformation time sequence phase features arranged according to an image acquisition sequence; and performing independent component analysis and inspection on the regional deformation graph to obtain a plurality of mutually independent signal source matrixes, thereby separating a deformation characteristic signal source by interpreting the plurality of signal source matrixes. According to the method, time-space separation can be carried out on the InSAR time sequence under the condition that no priori information exists.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing, and in particular to a method for separating the micro surface deformation phase of InSAR. Background Art

[0002] The micro surface deformation is often the primary manifestation of land subsidence. It is an intuitive perturbation form of the change of the subsurface tectonic physical mechanism and is the comprehensive result of various inducements, including over-exploitation of groundwater, tectonic activities, change of surface load, development of underground oil, gas and other resources, and activities of shallow soil layers such as frozen soil and expansive soil. The consequences related to land subsidence include: damage to infrastructure, roads and buildings; formation of ground fissures; and, coastal areas are more vulnerable to floods and land salinization. The more obvious land subsidence reflects the highly developed change of the subsurface physical mechanism, and the optimal treatment time for the factors caused by human activities has been missed. The accurate capture of micro surface deformation and the scientific inversion of human activity inducements are particularly crucial for the treatment of this environmental geological disaster of land subsidence.

[0003] At present, the conventional monitoring methods for regional surface deformation problems include GPS (global positioning system) measurement, leveling measurement, layer calibration technology, etc. However, the above monitoring methods have relatively high artificial equipment costs, can only obtain relatively discrete point deformation data, may miss deformation information, and it is difficult to conduct overall deformation monitoring on the research area. Synthetic aperture radar interferometry (InSAR) technology, as a new type of earth observation means, can overcome the deficiencies of traditional observation means, and obtain high-resolution surface deformation information through two observation values of intensity and phase, so as to achieve large-scale spatial detection and long-term precise monitoring.

[0004] In the study of large-area synthetic aperture radar interferometry (InSAR), separating different signal sources is a challenge, especially in areas with obvious changes in atmospheric conditions, active tectonic deformation and frequent human activities, where: (1) The change of atmospheric conditions mainly causes relatively serious tropospheric delay when the InSAR radar penetrates the troposphere. The tropospheric phase delay can reach up to 20 cm of apparent ground displacement and has a complex time evolution; (2) Active tectonic deformation is mainly long-wave deformation, which reflects micro surface deformation, but will interfere with obtaining the surface deformation component caused by human activities.

[0005] In summary, there are still great challenges in how to separate the surface deformation caused by different inducements, especially the spatio-temporal scale characteristics of the surface deformation caused by human activities. Summary of the Invention

[0006] The object of the present invention is to provide an InSAR surface micro-deformation phase separation technology based on independent component analysis to separate surface deformations caused by different incentives (in particular, the spatio-temporal scale characteristics of surface deformations caused by human activities).

[0007] To solve the above technical problems, an embodiment of the present invention provides a method for separating InSAR surface micro-deformation phases, including: obtaining radar remote sensing images of a target area, registering the image data by combining precise orbit data and ground digital elevation model data of the target area to obtain full-time domain image information; extracting an InSAR differential interferogram according to the full-time domain image information, and then performing phase feature extraction and optimization processing on the differential interferogram to calculate a regional deformation map, wherein each pixel in the regional deformation map contains deformation time-series phase features arranged in the order of image acquisition; performing independent component analysis and testing on the regional deformation map to obtain several mutually independent signal source matrices, and thus separating deformation feature signal sources by interpreting the several signal source matrices.

[0008] Preferably, in the process of registering radar image data, it includes: referring to the precise orbit data and the ground digital elevation model data, registering radar remote sensing images collected at different times to the radar coordinates of a preset reference image; resampling the images under the same reference coordinates according to the range direction and azimuth direction through a registration offset polynomial, and then performing registration in the overlapping area based on spectral diversity to eliminate satellite orbit errors, thereby forming the full-time domain image information.

[0009] Preferably, in the step of performing phase feature extraction and optimization processing on the differential interferogram to calculate a regional deformation map, it includes: preprocessing the InSAR differential interferogram; estimating the baseline information of the interferometric pair and performing baseline refinement processing on the preprocessed differential interferometric pair; performing phase unwrapping and error rejection processing on the differentially interferometric pair after baseline refinement; performing quality control on the processed interferometric pair; and obtaining the displacement time series of the interferometric pair after quality control based on the least squares criterion to form a regional deformation map arranged in a time series.

[0010] Preferably, in the step of preprocessing the InSAR differential interferogram, it includes: using the ground digital elevation model data to remove the topographic phase error in the InSAR differential interferometric pair; and performing filtering on the differential interferometric pair after topographic phase error processing by using Goldstein filtering.

[0011] Preferably, in the step of performing phase unwrapping and error rejection processing on the differentially interfered pair after baseline refinement, it includes: re-performing differential interference on the differentially interfered pair after baseline refinement, and estimating the average coherence of the current differentially interfered pair; determining a mask threshold according to the average coherence, so as to create a corresponding mask file based on the mask threshold; using the phase unwrapping method based on the minimum cost flow to unwrap the phase wrapped by the current interfered pair stack according to the mask file based on the mask threshold, so as to avoid the unwrapping error caused by low coherence.

[0012] Preferably, in the step of performing quality control on the processed interfered pair, it includes: removing the tropospheric delay error and trend error of the interfered pair after phase unwrapping; positioning and rejecting the interfered pairs with poor quality based on the closed phase difference of the current interfered pair loop, so as to complete the quality control of the interfered pair.

[0013] Preferably, in the process of performing independent component analysis and inspection on the regional deformation map to obtain a number of independent signal sources, it includes: writing the regional deformation map stack into a three-dimensional matrix and unfolding it two-dimensionally in space; performing feature component extraction and independent separation processing of feature components on the deformation time series phase matrix of each pixel point in the unfolded deformation map; using the projection matrix of a number of independent signal source vectors corresponding to each pixel point to separate the regional deformation map, so as to obtain signal sources respectively based on each independent component.

[0014] Preferably, in the process of performing feature component extraction and independent separation processing of feature components on the deformation time series phase matrix of each pixel point, it includes: the first step, centering the deformation time series phase matrix of each pixel point; the second step, determining the number of independent components of the regional deformation map through the eigenvalues and eigenvectors of each phase matrix, and determining the number of eigenvectors corresponding to the phase matrix of each pixel point according to the current number of independent components; the third step, performing non-Gaussianity measurement on the projection matrix of the number of eigenvectors corresponding to each pixel point, and updating the projection matrix of the number of eigenvectors by maximizing the non-Gaussianity measurement; the fourth step, orthogonalizing the projection matrix of the number of eigenvectors updated above, so as to ensure the mutual independence between the eigenvectors; the fifth step, checking whether the current independent component decomposition algorithm converges according to the projection matrix of the number of orthogonalized eigenvectors corresponding to each pixel point. If it converges, taking the projection matrix of the number of eigenvectors corresponding to each pixel point as the independent signal source vectors corresponding to the corresponding pixel points. Among them, if it does not converge, returning to the second step to re-perform signal decomposition processing by adjusting the number of independent components until the decomposition result converges.

[0015] Preferably, in the second step, it includes: Step S1, calculating the eigenvalues and eigenvectors of the covariance matrix corresponding to each phase matrix, and based on this, combining the number of eigenvalues and the number of important eigenvalues of each pixel point to determine the current component number threshold and the preset difference threshold; Step S2, respectively calculating, for each phase matrix, the difference between the sum of the eigenvalues before the current component number threshold and the sum of all the subsequent eigenvalues, so as to determine whether the difference corresponding to each pixel point is greater than the preset difference threshold. If all exceed, then use the current component number threshold as the current independent component number. Among them, if there are pixel points that do not exceed, return to Step S1 to readjust the component number threshold.

[0016] Preferably, in the step of separating the deformation feature signal source by interpreting the several signal source matrices, it includes: respectively performing DEM error source analysis on each separated signal source matrix; according to the prior information about the target area, estimating the fitting accuracy of the DEM error of each signal source matrix with the actual allowable error of the target area, and diagnosing whether the current independent component analysis is correct. Among them, if the decomposition result does not meet the accuracy requirements, perform independent component analysis processing again.

[0017] Compared with the prior art, one or more embodiments in the above solution may have the following advantages or beneficial effects:

[0018] The present invention proposes a method for separating the InSAR surface micro-deformation phase. This method performs spatio-temporal separation on the InSAR time series without any prior information to obtain m independent components, including m - 1 deformation components and one noise component. After removing the noise, the application of independent component analysis can improve the deformation accuracy. In addition, different deformation components also reveal different geomechanical processes in the study area, thereby helping researchers understand and interpret the geophysical change process below the surface.

[0019] Other features and advantages of the present invention will be described in the subsequent specification, and, in part, will become apparent from the specification, or will be understood by implementing the present invention. The objectives and other advantages of the present invention can be achieved and obtained through the structures specifically pointed out in the specification, claims, and drawings. Description of the Drawings

[0020] The drawings are used to provide a further understanding of the present invention, and constitute a part of the specification. They are used together with the embodiments of the present invention to explain the present invention, and do not constitute a limitation to the present invention. In the drawings:

[0021] Figure 1 It is a step diagram of the method for separating the InSAR surface micro-deformation phase according to the embodiment of the present application.

[0022] Figure 2 It is a schematic diagram of the specific process in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0023] Figure 3 It is a schematic diagram of the SAR image coverage range of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0024] Figure 4 It is a schematic diagram of the terrain features and geographical location of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0025] Figure 5 It is an example diagram of the spatio-temporal baseline distribution of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0026] Figure 6 It is an example diagram of the regional deformation map of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0027] Figure 7 It is an example diagram of the first independent component signal source in the regional deformation map of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0028] Figure 8 It is an example diagram of the second independent component signal source in the regional deformation map of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application.

[0029] Figure 9 It is an example diagram of the third independent component signal source in the regional deformation map of the target area to be evaluated in the method for separating the InSAR surface micro-deformation phase in the embodiment of the present application. Detailed implementation manners

[0030] The following will combine the accompanying drawings and embodiments to detail the implementation manners of the present invention, so as to fully understand how the present invention uses technical means to solve technical problems and the implementation process of achieving technical effects and implement accordingly. It should be noted that as long as there is no conflict, the various embodiments in the present invention and the various features in each embodiment can be combined with each other, and the formed technical solutions are all within the protection scope of the present invention.

[0031] In addition, the steps illustrated in the process flow diagrams of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions. And although a logical order is illustrated in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that herein.

[0032] The terms used herein are merely for the purpose of describing particular embodiments and are not intended to limit the exemplary embodiments. Unless the context clearly dictates otherwise, the singular forms "a", "an" used herein are also intended to include the plural. It should also be understood that the terms "comprises" and / or "comprising" specify the presence of the stated features, integers, steps, operations, units and / or components, and do not preclude the presence or addition of one or more other features, integers, steps, operations, units, components and / or combinations thereof.

[0033] The minute surface deformation is often the primary manifestation of land subsidence, an intuitive perturbation form of the physical mechanism change under the surface, and the comprehensive result of multiple inducing factors, including over-exploitation of groundwater, tectonic activities, surface load changes, development of underground oil, gas and other resources, and activities of shallow soil layers such as frozen soil and expansive soil. The consequences related to land subsidence include: damage to infrastructure, roads and buildings; formation of ground fissures; and, coastal areas are more vulnerable to floods and land salinization. The more obvious land subsidence reflects the highly developed physical mechanism change under the surface, and the optimal treatment time for the factors caused by human activities has been missed. The accurate capture of minute surface deformation and the scientific inversion of the inducing factors of human activities are particularly crucial for the treatment of this environmental geological disaster of land subsidence.

[0034] Currently, the conventional monitoring methods for regional surface deformation problems include GPS (global positioning system) measurement, leveling measurement, layer calibration technology, etc. However, the above monitoring methods have relatively high costs for artificial equipment, can only obtain relatively discrete point deformation data, may miss deformation information, and it is difficult to conduct overall deformation monitoring of the research area. The synthetic aperture radar interferometry (InSAR) technology, as a new means of earth observation, can overcome the deficiencies of traditional observation means, obtain high-resolution surface deformation information through two observation values of intensity and phase, and achieve large-scale spatial detection and long-term precise monitoring.

[0035] In the study of large - area interferometric synthetic aperture radar (InSAR), separating different signal sources is a challenge, especially in areas with significant atmospheric condition changes, active tectonic deformation, and frequent human activities. Among them: (1) Atmospheric condition changes mainly cause severe tropospheric delays when the InSAR radar penetrates the troposphere. The tropospheric phase delay can reach up to 20 cm of apparent ground displacement and has a complex temporal evolution. (2) Active tectonic deformation is mainly long - wave deformation, which reflects small surface deformations but interferes with obtaining the component of surface deformation caused by human activities.

[0036] In summary, there are still great challenges in separating the spatio - temporal scale characteristics of surface deformations caused by different inducements, especially those caused by human activities.

[0037] To solve the problems in the above - mentioned background technology, an embodiment of the present application proposes a method for separating the phase of small surface deformations in InSAR. The method includes: performing SAR data pre - processing by registering, cropping, and removing orbital errors on SAR satellite image data; obtaining and quality - controlling the time series of regional InSAR deformations; two - dimensional unfolding of the InSAR deformation map stack; calculating the eigenvalues and eigenvectors of the covariance matrix of the mixed observation matrix to determine the number of independent signal sources of the mixed matrix; non - Gaussianizing the projected mixed observation matrix and using the negative entropy function to measure the non - Gaussianity of the data, so as to perform independent component analysis on the deformation observation results of the research area obtained by the InSAR technology, and interpreting the physical mechanisms reflected by the independent signal sources according to the regional prior information.

[0038] Example 1

[0039] Figure 1 It is a step diagram of the method for separating the phase of small surface deformations in InSAR of the embodiment of the present application. Figure 2 It is a schematic diagram of the specific process in the method for separating the phase of small surface deformations in InSAR of the embodiment of the present application. The following combines Figure 1 and Figure 2 to illustrate the specific step - by - step process of the method for separating the phase of small surface deformations in InSAR (hereinafter referred to as the "deformation phase separation method") described in the embodiment of the present invention.

[0040] Step S110: Obtain radar remote - sensing images collected at different times for the target area. Combine the precise orbit data and ground digital elevation model data of the target area to register the image data and obtain full - time - domain image information.

[0041] In step S110, Sentinel-1A radar remote sensing images, precise orbit data, and ground digital elevation model data (Digital Elevation Model, DEM) of the target area to be studied are obtained. In the GAMMA software, registration, cropping, and orbit error removal are performed on the SAR satellite radar remote sensing image data (see Figure 3 ).

[0042] It should be noted that step S110 requires obtaining radar remote sensing images collected at different times, and each radar remote sensing image can cover the complete topographic features and geographical locations of the target area (see Figure 4 ). In addition, during the preprocessing of InSAR data in step S110, the precise orbit data of the target area used is the Sentinel-1A precise orbit data provided by the European Space Agency (ESA), and the ground digital elevation model data of the target area used is the SRTM 30-meter resolution digital elevation model provided by the National Aeronautics and Space Administration (NASA) of the United States.

[0043] During the registration of radar image data collected at different times, first, referring to the precise orbit data and ground digital elevation model data of the target area, the radar remote sensing images collected at different times are registered to the radar coordinates of a preset reference image; then, according to the range direction and azimuth direction, resampling is performed on the images under the same reference coordinate through a registration offset polynomial, and registration in the overlapping area is performed based on spectral diversity to eliminate satellite orbit errors (phase discontinuity), thereby forming the full-time domain image information.

[0044] Generally, the Sentinel-1A image data in the IW imaging mode realizes a 250-km-wide scanning imaging through TOPS imaging, but in the image assembly process, the TOPS mode is used. In the TOPS mode, there is a relatively large change in the Doppler center frequency between different bursts, and the registration error will cause a shift in the interference phase, which leads to an obvious phase jump phenomenon easily occurring between different bursts. Therefore, the embodiment of the present invention needs to use step S110 to carry out data preprocessing, add precise orbit parameters to all images, and accurately register them to the radar coordinates of a preset reference image.

[0045] To achieve precise registration between SLC images, the initial lookup table is calculated using the image parameter file of the radar image acquisition device and the DEM, and a series of images are resampled into the geometric framework of the reference image. Then, combined with the range direction and azimuth direction, the registration offset polynomial is used to resample the images. On this basis, a small amount of data still needs to be further registered in the overlapping area based on spectral diversity, so as to eliminate the error caused by the phase discontinuity between bursts in the differential interferogram, and thus obtain the full-time-domain image information and enter step S120.

[0046] In step S120, according to the full-time-domain image information obtained in step S110, the InSAR differential interferogram is extracted, and then the phase characteristics of the extracted differential interferogram are extracted and optimized, so as to calculate the regional deformation map arranged in time series. Among them, each pixel in the regional deformation map contains the deformation time-series phase characteristics arranged in the order of image acquisition (time).

[0047] In step S120, first, the full-time-domain image information is converted into an InSAR differential interferogram (InSAR differential interferogram pair). Specifically, according to the baseline data of the SAR satellite radar remote sensing image data, an interference network of the SAR data is constructed, and based on the constructed interference network, differential interference processing is performed on the full-time-domain image information to obtain the InSAR differential interferogram pair.

[0048] Next, step S120 also extracts the phase characteristics of the differential interferogram (InSAR differential interferogram pair) and performs optimization processing, so as to calculate the regional deformation map.

[0049] Step 1: Preprocess the InSAR differential interferogram. In step 1, first, the ground digital elevation model data is used to remove the topographic phase error in the InSAR differential interferogram pair; then, the Goldstein filter is used to filter the differential interferogram pair after the topographic phase error processing. Use the DEM data to remove the topographic phase error in the InSAR differential interferogram pair; use the Goldstein filter to filter the stack of interferogram pairs after the topographic phase error removal processing, so as to improve the coherence of the interferogram pair.

[0050] Step 2: Estimate the baseline information of the interferogram pair and perform baseline refinement processing on the preprocessed differential interferogram pair. In step 2, first, estimate the baseline information of the interferogram pair in the current target area to form an InSAR baseline map, see Figure 5 . Then, perform windowing processing according to the estimated InSAR baseline map to achieve the purpose of baseline refinement, so as to ensure that the interferogram pair is not affected by baseline errors.

[0051] Step 3: Perform phase unwrapping and error elimination on the differentially interfered pairs after baseline refinement. In Step 3, it includes: re-performing differential interference on the differentially interfered pairs after baseline refinement, and estimating the average coherence of the current differentially interfered pair; determining a mask threshold according to the average coherence, so as to create a corresponding mask file based on the current mask threshold; finally, using the MCF phase unwrapping method based on the minimum cost flow to unwrap the phase wrapped by the current stack of interfered pairs after double differential interference, so as to avoid the unwrapping error caused by low coherence.

[0052] Specifically, re-perform differential interference on the SAR images after baseline refinement, and estimate the average coherence of the interfered pairs; determine the mask threshold based on the coherence mean file of the interfered pair stack, so as to create a mask file, and the threshold generally refers to the range of 0.4 - 0.7. Then, use the phase unwrapping method based on the minimum cost flow to unwrap the phase wrapped by the interfered pair stack. Among them, during the unwrapping process, use the mask file based on the coherence threshold, so as to avoid the unwrapping error caused by low coherence to a certain extent.

[0053] Coherence is an important indicator to measure the quality of the InSAR interferogram. In the embodiment of the present invention, through the calculation of coherence, without relying on external measurement data and combining on-site investigation data, a relatively reliable InSAR deformation monitoring result can be obtained. In addition, in order to be able to select unwrapping points with stable and good ground coherence, the phase unwrapping method based on the minimum cost flow (MCF) is used to unwrap the phase wrapped by the current stack of interfered pairs, so as to avoid the interference of low coherence on unwrapping.

[0054] Next, Step 4: Perform quality control on the processed interfered pairs. In Step 4, it includes: removing the tropospheric delay error and (stratospheric) trend error of the interfered pairs after phase unwrapping; then, based on the closure phase difference of the current interfered pair loop, locate and eliminate the interfered pairs with poor quality, so as to complete the quality control of the interfered pairs. Specifically, after removing the tropospheric delay error and trend error of the interfered pairs and completing the error removal; based on the closure phase difference of the current interfered pair loop, locate and eliminate the interfered pairs with a phase closure difference greater than 3 times the mean square error.

[0055] After completing the quality control of the interfered pairs, Step 5: Based on the least squares criterion, obtain the displacement time series of the interfered pairs after quality control, and form a regional deformation map arranged in time series. According to the interfered pairs after quality control, based on the least squares criterion, obtain the displacement time series of the current interfered pair, and obtain the mixed deformation matrix (i.e., the regional deformation map), as Figure 6 shown.

[0056] In this way, after the InSAR time series deformation solution is performed on the preprocessed full-time domain image information in step S120 in the embodiment of the present invention, a regional deformation map represented by a mixed deformation matrix is generated, and thus step S130 is entered.

[0057] In step S130, independent component analysis and inspection are performed on the regional deformation map obtained in step S120 to obtain a number of mutually independent signal source matrices, and thus the deformation feature signal sources are separated by interpreting the number of signal source matrices.

[0058] In step S130, the stack of the regional deformation maps represented in the form of a mixed deformation matrix obtained in step S120 is written into a three-dimensional matrix and two-dimensionally expanded in space. Specifically, the stack of deformation maps arranged in time series is written into a three-dimensional matrix and two-dimensionally expanded in space, where the row vector of each pixel point in the regional deformation map represents the deformation time series (phase matrix) of the current single pixel.

[0059] Then, feature component extraction and feature component independence separation processing are performed on the deformation time series phase matrix of each pixel point in the two-dimensionally expanded deformation map. ICA decomposition is respectively performed on the deformation time series matrix of each pixel point, so that after the decomposition results of each deformation time series matrix meet the preset convergence conditions, a corresponding number of mutually independent signal source vectors can be obtained for each pixel point.

[0060] Specifically, in the first step, the deformation time series phase matrix of each pixel point is centralized, and the eigenvalues and eigenvectors of the covariance matrix of each deformation time series phase matrix are calculated.

[0061] In the second step, the number of independent components of the regional deformation map is determined through the calculated eigenvalues and eigenvectors of each phase matrix, and the number of eigenvectors (matrices) corresponding to the independent components of the phase matrix of each pixel point is determined according to the current number of independent components. It should be noted that the number of independent components of the regional deformation map described in the embodiment of the present invention is the number of independent components jointly possessed by each pixel point in the regional deformation map.

[0062] In the second step, the number of independent components of the current regional deformation map is solved according to the following process, and the eigenvector matrix of the current number of independent components is generated for each pixel point, where the eigenvector matrix corresponding to each pixel point represents the deformation time series phase feature of the pixel point.

[0063] Specifically, in the second step, in step S1, calculate the eigenvalues and eigenvectors of the covariance matrix corresponding to each phase matrix, and determine the current component number threshold and the preset difference threshold according to the eigenvalues and eigenvectors corresponding to each phase matrix, in combination with the number of eigenvalues of each pixel point and the number of important eigenvalues. According to the eigenvalues and eigenvectors of the phase covariance matrix of each pixel point (one pixel point as a group), each group retains the most important k eigenvectors (usually, k is less than n, where n represents the minimum value of the number of eigenvalues of all pixel points), determine a component number threshold k, and at the same time, set a difference threshold.

[0064] Step S2: For each phase matrix, calculate the difference between the sum of the eigenvalues before the current component number threshold and the sum of all the eigenvalues after (after the current component number threshold), so as to judge whether the difference corresponding to each pixel point is greater than the current preset difference threshold. If all exceed (that is, if the differences corresponding to all pixel points exceed the current preset difference threshold), then use the current component number threshold as the current independent component number. In addition, if there are pixel points that do not exceed (that is, if there is at least one pixel point whose corresponding difference does not exceed the current preset difference threshold), then return to the above step S1 to readjust the component number threshold until the above preset difference threshold condition is met and the corresponding current independent component number is determined.

[0065] For each phase matrix, determine the number of independent signal sources of the current mixing matrix according to the following process: if the sum of the first k eigenvalues minus the sum of the subsequent n - k eigenvalues is greater than the preset difference threshold, then select this k as the current independent component number. Thus, the eigenvectors corresponding to the k eigenvalues of the current pixel point are obtained. In this way, it is determined that the current deformation mixing signal matrix is composed of k independent components. Therefore, for each pixel point, the corresponding k eigenvectors (matrices) are obtained, and then enter the third step.

[0066] In the third step, perform non-Gaussianity measurement on the projection matrices (for example: the unit vectors of each eigenvector) of the number of eigenvectors corresponding to each pixel point, and update the projection matrices of the number of eigenvectors by maximizing the non-Gaussianity measurement.

[0067] According to information theory, among all random variables with equal variance, the entropy of Gaussian variables is the largest, so negative entropy can be used to measure non-Gaussianity. According to the central limit theorem, if a random variable X is composed of the sum of many independent random variables Si (i = 1, 2, 3..., N), as long as Si has a finite mean and variance, then no matter what kind of distribution it is, the random variable X is closer to the Gaussian distribution than Si. Therefore, when the Gaussianity measure reaches the maximum, it means that the separation of each independent component is completed. In this way, by maximizing the non-Gaussianity measure of the projection matrix of the number of independent component eigenvectors, the number of independent component eigenvectors of each pixel point are separated and updated, so that the separated (independent) number of independent component eigenvectors can be obtained for each pixel point.

[0068] The fourth step is to orthogonalize the projection matrix of the updated independent component number eigenvectors to ensure the independence of each eigenvector. For each pixel point, the projection matrix of the updated independent component number eigenvectors is orthogonalized, and for each pixel point, an independent component number eigenvector after orthogonalization is obtained.

[0069] The fifth step is to check whether the current independent component decomposition algorithm has converged according to the projection matrix of the orthogonalized independent component number eigenvectors corresponding to each pixel point. If converged, the projection matrix of the current independent component number eigenvectors corresponding to each pixel point is used as the (k) mutually independent signal source vectors corresponding to the corresponding pixel point. In addition, if it does not converge, return to the above second step and re-carry out the signal decomposition process by re-adjusting the current number of independent components until the decomposition result converges.

[0070] In this way, after completing the ICA decomposition of the deformation time series phase matrix of each pixel point, a projection matrix of mutually independent signal source vectors corresponding to the number of independent components is obtained for each pixel point.

[0071] Finally, the projection matrix of several (the number of current independent components) mutually independent signal source vectors corresponding to each pixel is used to separate the regional deformation map solved in the current step S120 to obtain the signal sources based on each independent component, see Figure 7 , Figure 8 and Figure 9 .

[0072] Furthermore, step S130 also separates the deformation feature signal source by interpreting a number (the number of current independent components) of signal source matrices.

[0073] Specifically, first, perform DEM error source analysis on each separated signal source matrix respectively. Then, according to the prior information about the target area (for example: actual precision requirements), estimate the fitting precision between the DEM error of each signal source matrix and the actual allowable error of the target area respectively, and diagnose whether the current independent component analysis is correct (that is, perform the goodness-of-fit F-test on each signal source matrix). Among them, if there is a situation where the precision requirements are not met in each separated signal source matrix, re-perform the independent component analysis process.

[0074] If the precision requirements are met in each separated signal source matrix, directly use each currently separated signal source matrix as the number of deformation feature signal sources corresponding to the current independent components. Among them, among several (m) deformation feature signal sources, it includes: m - 1 deformation components and 1 noise component.

[0075] Example 2

[0076] Based on the deformation phase separation method formed in the above-mentioned Embodiment 1, the embodiment of the present invention applies the above deformation phase separation method to the surface deformation monitoring and signal separation of a certain area. The topographic features and geographical location of the target area, such as Figure 4 shown, to illustrate the InSAR surface micro-deformation separation technology provided by the present invention based on independent component analysis, specifically including the following steps;

[0077] Step T1, preprocessing of SAR data in the target area:

[0078] Obtain the Sentinel-1A radar remote sensing image, precise orbit data and ground digital elevation model data (Digital Elevation Model, DEM) of the target area. In the GAMMA software, perform steps such as registration, cropping, and orbit error removal on the SAR satellite image data.

[0079] Specifically, the data mainly used in this research on emergency mapping of earthquake disaster events are: Sentinel-1A radar remote sensing images before and after the earthquake, and the image coverage range is as Figure 3 shown. In the process of InSAR data preprocessing, auxiliary data such as Sentinel-1A precise orbit data provided by the European Space Agency (ESA) and SRTM 30-meter resolution digital elevation model provided by the National Aeronautics and Space Administration (NASA) of the United States are used.

[0080] The Sentinel-1A image data in IW imaging mode realizes scanning imaging with a swath width of 250 km through TOPS imaging. However, in the TOPS mode, there is a relatively large change in the Doppler center frequency between different bursts, and registration errors will cause shifts in the interferometric phase, which leads to obvious phase jumps between different bursts. Therefore, data preprocessing is required. Precise orbit parameters need to be added to all images, and they need to be accurately registered to the radar coordinates of a reference image. To achieve precise registration between SLC images, an initial lookup table is calculated using the master and slave image parameter files and the DEM, and the slave image is resampled to the geometric framework of the reference image. The image is resampled by combining the range and azimuth registration offset polynomials. On this basis, a small amount of data needs to be further registered in the overlapping area based on spectral diversity to eliminate the errors caused by the phase discontinuity between bursts in the differential interferogram.

[0081] Step T2, obtaining regional deformation information by time-series InSAR:

[0082] To better reduce the influence of spatial decorrelation, the small baseline subset InSAR technique is used, and image pairs are selected according to the short baseline principle. When solving for deformation, the singular value decomposition (SVD) method is used to connect the independent data sets separated due to baseline limitations, and the surface deformation solution based on the least squares criterion is solved. The InSAR baseline map is as Figure 5 shown. The main processing steps are as follows: (1) Generate interferograms. In this paper, the spatial baseline threshold is set to ±150 m and the temporal baseline threshold is set to 36 d, and a total of 156 interferogram pairs are generated; (2) Calculate coherence. Coherence is an important indicator for measuring the quality of InSAR interferograms. By calculating the coherence, reliable InSAR deformation monitoring results can be obtained without relying on external measurement data and combined with on-site investigation data. (3) Select unwrapping points with stable and good ground coherence and perform phase unwrapping based on the minimum cost flow (MCF); (4) Remove the trend and atmospheric errors; (5) Solve the weighted least squares solution of the deformation time-series phase through the SVD method. The deformation solution results are as Figure 6 shown.

[0083] Step T3, phase separation technology based on independent component analysis:

[0084] Write the stack of deformation maps arranged in time series into a three-dimensional matrix and expand it two-dimensionally in space. Each row vector represents the deformation time series of a single pixel. Determine the number of independent signal sources in the mixing matrix: For the two-dimensional matrix, that is, remove the mean value row by row for the implementation sequence of each pixel, perform centering processing, calculate the eigenvalues and eigenvectors of the covariance matrix, and retain the most important 3 features (usually k should be less than n). Select a threshold. By subtracting the sum of the last n - 3 eigenvalues from the sum of the first 3 eigenvalues and being greater than this threshold, it is determined that the deformation mixing signal matrix is composed of 3 independent components.

[0085] Calculate the non-Gaussianity of the projected data: Use the negative entropy function to measure the non-Gaussianity of the projected data. Update the value of the projection matrix W by maximizing the non-Gaussianity metric function. Orthogonalize the updated projection matrix W to ensure that each component is independent of each other. Convergence check: Check whether the algorithm converges. If the convergence condition has not been reached, return to the non-Gaussianity metric step to continue the iteration. Separate the mixed signals using the estimated projection matrix to obtain 3 independent signal sources. According to the regional prior information, interpret the physical mechanisms reflected by the 3 signal sources. The 3 separated independent deformation signal sources are respectively as Figure 7 、 8 、9 shown. It can be seen that Figure 7 relatively obvious periodic deformations are accurately extracted.

[0086] The present invention discloses a method for separating the InSAR surface micro-deformation phase. This method performs spatio-temporal separation on the InSAR time series without any prior information to obtain m independent components, including m - 1 deformation components and one noise component. After removing the noise, the application of independent component analysis can improve the deformation accuracy. In addition, different deformation components also reveal different geomechanical processes in the study area, thereby helping researchers understand and interpret the geophysical change process below the surface.

[0087] The above is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those familiar with the technology within the technical scope disclosed by the present invention should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

[0088] In the description of the present invention, unless otherwise specified, the meaning of "a plurality of" is two or more; the orientation or positional relationship indicated by terms such as "upper", "lower", "left", "right", "inner", "outer", "front end", "rear end", "head", "tail", etc. is based on the orientation or positional relationship shown in the drawings. It is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as a limitation to the present invention. In addition, terms such as "first", "second", "third", etc. are only used for descriptive purposes and cannot be construed as indicating or implying relative importance.

[0089] In the description of the present invention, it should be noted that, unless otherwise clearly specified and defined, the terms "connected" and "coupled" should be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be directly connected or indirectly connected through an intermediate medium. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0090] It should be understood that the embodiments disclosed in the present invention are not limited to the specific structures, processing steps or materials disclosed herein, but should extend to equivalent alternatives of these features understood by those of ordinary skill in the relevant art. It should also be understood that the terms used herein are only for the purpose of describing specific embodiments and do not imply limitation.

[0091] The "one embodiment" or "embodiment" mentioned in the specification means that the specific features, structures or characteristics described in connection with the embodiment are included in at least one embodiment of the present invention. Therefore, the phrases "one embodiment" or "embodiment" appearing throughout the specification do not necessarily all refer to the same embodiment.

[0092] Although the embodiments disclosed in the present invention are as above, the content described is only the embodiment adopted for the convenience of understanding the present invention and is not intended to limit the present invention. Any person skilled in the art within the technical field to which the present invention pertains can make any modifications and changes in the form of implementation and details without departing from the spirit and scope disclosed by the present invention. However, the scope of patent protection of the present invention shall still be subject to the scope defined by the appended claims.

Claims

1. A method for separating the InSAR surface micro-deformation phase, characterized in that Including: Obtain radar remote sensing images of the target area, combine the precise orbit data and ground digital elevation model data of the target area, register the image data, and obtain full-time domain image information; Extract the InSAR differential interferogram according to the full-time domain image information, then extract the phase characteristics of the differential interferogram and perform optimization processing, so as to calculate the regional deformation map, wherein each pixel in the regional deformation map contains the deformation time-series phase characteristics arranged according to the image acquisition order; Perform independent component analysis and inspection on the regional deformation map to obtain several mutually independent signal source matrices, so as to separate the deformation characteristic signal sources by interpreting the several signal source matrices.

2. The method according to claim 1, characterized in that, During the process of registering the radar image data, including: Referring to the precise orbit data and the ground digital elevation model data, register the radar remote sensing images collected at different times to the preset reference image radar coordinates; According to the range direction and azimuth direction, resample the images under the same reference coordinates through the registration offset polynomial, and then perform overlapping area registration based on spectral diversity to eliminate satellite orbit errors, so as to form the full-time domain image information.

3. The method according to claim 1 or 2, wherein In the step of extracting the phase characteristics of the differential interferogram and performing optimization processing to calculate the regional deformation map, including: Preprocess the InSAR differential interferogram; Estimate the baseline information of the interferometric pair and perform baseline refinement processing on the preprocessed differential interferometric pair; Perform phase unwrapping and error rejection processing on the differentially interferometric pair after baseline refinement; Perform quality control on the processed interferometric pair; Based on the least squares criterion, obtain the displacement time series of the interferometric pair after quality control, and form a regional deformation map arranged according to the time series.

4. The method according to claim 3, wherein In the step of preprocessing the InSAR differential interferogram, including: Use the ground digital elevation model data to remove the topographic phase error in the InSAR differential interferometric pair; Adopt Goldstein filtering to filter the differential interferometric pair after topographic phase error processing.

5. The method according to claim 3 or 4, characterized in that, In the step of performing phase unwrapping and error rejection processing on the differentially interferometric pair after baseline refinement, including: Perform differential interference on the differentially interferometric pair after baseline refinement again, and estimate the average coherence of the current differential interferometric pair; According to the average coherence, determine the mask threshold, and thus make a corresponding mask file based on the mask threshold; According to the mask file based on the mask threshold, use the phase unwrapping method based on the minimum cost flow to unwrap the phase wrapped by the current interferometric pair stack to avoid unwrapping errors caused by low coherence.

6. The method according to any one of claims 3 to 5, characterized in that In the step of performing quality control on the processed interferometric pair, including: Remove the tropospheric delay error and trend error of the interferometric pair after phase unwrapping; Based on the closed phase difference of the current interferometric pair loop, locate and eliminate the interferometric pairs with poor quality, so as to complete the quality control of the interferometric pair.

7. The method according to any one of claims 1 to 6, characterized in that In the process of performing independent component analysis and inspection on the regional deformation map to obtain several mutually independent signal sources, including: Write the regional deformation map stack into a three-dimensional matrix and expand it two-dimensionally in space; Extract the eigen-components and perform independent component separation on the deformation time-series phase matrix of each pixel in the expanded deformation map; Use the projection matrices of a number of mutually independent source vectors corresponding to each pixel to separate the regional deformation map, and obtain the source signals based on each independent component.

8. The method according to claim 7, characterized in that In the process of extracting eigen-components and performing independent component separation on the deformation time-series phase matrix of each pixel, it includes: First step, centralize the deformation time-series phase matrix of each pixel; Second step, determine the number of independent components of the regional deformation map through the eigenvalues and eigenvectors of each phase matrix, and determine the number of eigenvectors corresponding to the independent components of each pixel phase matrix according to the current number of independent components; Third step, measure the non-Gaussianity of the projection matrices of the number of eigenvectors corresponding to each pixel, and update the projection matrices of the number of eigenvectors corresponding to the independent components by maximizing the measure of non-Gaussianity; Fourth step, orthogonalize the updated projection matrices of the number of eigenvectors corresponding to the independent components to ensure the mutual independence between the eigenvectors; Fifth step, according to the projection matrices of the number of orthogonalized eigenvectors corresponding to each pixel, check whether the current independent component decomposition algorithm converges. If it converges, use the projection matrices of the number of current independent components corresponding to each pixel as the mutually independent source vectors corresponding to the respective pixels, where If it does not converge, return to the second step, and re-perform signal decomposition by adjusting the number of independent components until the decomposition result converges.

9. The method according to claim 8, wherein In the second step, it includes: Step S1: Calculate the eigenvalues and eigenvectors of the covariance matrix corresponding to each phase matrix. Based on this, combined with the number of eigenvalues and the number of important eigenvalues of each pixel, determine the current component number threshold and the preset difference threshold; Step S2: Calculate the difference between the sum of the eigenvalues before the current component number threshold and the sum of all the eigenvalues after for each phase matrix respectively, so as to judge whether the difference corresponding to each pixel is greater than the preset difference threshold. If all exceed, use the current component number threshold as the current number of independent components, where If there are pixels that do not exceed, return to step S1 to re-adjust the component number threshold.

10. The method according to any one of claims 1 to 9, characterized in that In the step of separating the deformation feature source signals by interpreting the several source matrices, it includes: Perform DEM error source analysis on each separated source matrix respectively; According to the prior information about the target area, estimate the fitting accuracy of the DEM errors of each source matrix with the actual allowable error of the target area, and diagnose whether the current independent component analysis is correct, where If the decomposition result does not meet the accuracy requirements, re-perform independent component analysis.

Images