A sequential polarimetric phase optimization method and system for near real-time DS-InSAR processing

Abstract

The present invention provides a sequential polarization phase optimization method and system for near-real-time DS-InSAR processing, which includes preparing environmental files for sequential processing, dividing time series multi-polarization SAR data into spatial blocks containing multiple mini time series subsets according to the size of the segmented time series subsets and the number of segmented spatial blocks; checking the environmental files for sequential processing, reading SAR data files and preprocessing; identifying statistically homogeneous pixels (SHP) of time series subsets, including performing SHP homogeneous pixel identification on newly added time series SAR image subsets during sequential processing and using them as updated samples for phase optimization; performing SETP-EMI phase optimization processing on time series subsets, including performing TP-EMI-based phase optimization processing on the read multi-polarization SAR data and the identified homogeneous pixels, and performing data compression processing on the original multi-polarization SLC data after completing the phase optimization; sequential processing and linking of the phase optimization of time series subsets; and phase splicing of the total time series optimization of the spatial blocks to output an interference phase map.

Description

Technical Field

[0001] The present invention belongs to the field of synthetic aperture radar interferometry (InSAR) surface monitoring, and in particular relates to a distributed scatterer InSAR (DS-InSAR) near real-time update technology solution for natural surface deformation monitoring. Background Art

[0002] Time-series differential synthetic aperture radar interferometry (D-InSAR) has been widely used in fields such as early landslide detection and monitoring, mining area deformation, urban surface subsidence, and earthquake deformation. Dense vegetation cover in large mountainous areas leads to significant decorrelation in radar interferograms, necessitating the use of distributed scatterer InSAR (DS-InSAR) technology to improve phase quality. The phase optimization process requires retrieving the time-series phases (homogeneous pixels) of ground targets with similar scattering mechanisms in local space and utilizing the spatiotemporal characteristics of the time-series phases to minimize decorrelation noise while maintaining the accuracy of the phase estimate. Therefore, this technique can significantly increase the spatial density of measured scatterers, enabling more effective time-series D-InSAR deformation analysis.

[0003] The development of DS-InSAR has undergone significant evolution over the past few decades. Ferretti et al. introduced a maximum likelihood phase linking algorithm, known as the phase triangulation algorithm (PTA), under the assumption of a complex circular Gaussian (CCG) distribution. Building on this work, Ansari et al. proposed a computationally efficient method, known as the eigendecomposition-based interferometric phase maximum likelihood estimator (EMI). This method aims to achieve computational efficiency while optimizing the estimation of the interferometric phase. These algorithms have helped advance phase linking techniques, improve phase quality, and enable more accurate deformation analysis in various InSAR applications.

[0004] However, the above-mentioned techniques are mainly based on single-polarization data and cannot obtain satisfactory results in some natural scenes covered by vegetation. Its limitation stems from the fact that the information provided by single-polarization data is insufficient to accurately characterize the target scattering behavior. Therefore, these single-polarization methods are challenging in capturing the full complexity of the scattering process and accurately estimating the phase information in different scenes. Due to the different sensitivities of polarization channels to different scattering characteristics and different shapes of objects, there are large differences in the interferometric phase coherence of different polarization channels. The interferometric phase quality obtained using only a single polarization channel is not optimal. To address this challenge, the additional polarization information of synthetic aperture radar (SAR) satellites can be used to alleviate the decorrelation problem and obtain more accurate and reliable phase series estimates. At present, related studies have adopted two different strategies to incorporate multi-polarization information into time series D-InSAR. The first method is based on polarization optimization, while the second method is based on stacking interferometric coherence matrices.

[0005] Polarization-based phase optimization is achieved by selecting the channel with the highest temporal coherence. For example, the exhaustive search polarimetric optimization (ESPO) method uses two evaluation metrics to optimize the selection of polarization channels. By using the metric with the best average coherence for distributed scatterers (DS) to calculate the polarimetric projection angle, time-series D-InSAR results are further improved. However, the ESPO method requires a large number of brute-force search calculations, resulting in extremely low computational efficiency and unsuitable for large-scale data processing. To address the limitations of ESPO, a suboptimal solution called PolPSI based on coherence matrix decomposition and a polarization phase optimization method based on scattering mechanism filtering were subsequently proposed. This method aims to effectively identify and process persistent scatterers (PS) and distributed scatterers (DS) in an adaptive manner. Compared to the polarization optimization method, another method that utilizes the total power (TP) coherence matrix method simply superimposes the interferometric coherence matrices from different polarization channels to obtain a more accurate interferometric coherence matrix. This method allows the phase sequence to be obtained directly from the off-diagonal elements of the TP coherence matrix or by applying a maximum likelihood estimation algorithm. The TP coherence matrix-based maximum likelihood estimation algorithm, also known as TP-EMI, has been proven to be an effective method for multi-polarization phase optimization, surpassing the ESPO method in terms of efficiency and accuracy. Unfortunately, the TP-EMI method lacks the ability to dynamically process newly added SAR data sequentially.

[0006] In the era of SAR big data, maintaining a balance between phase optimization accuracy and efficient dynamic processing is crucial. This is especially true for nonlinearly deformable surfaces, which require continuous dynamic monitoring. To address this issue, the Sequential Estimator (SE) recursively processes newly added datasets, reducing data storage requirements and ensuring phase optimization accuracy for new data. However, the SE is designed for single-polarization data and does not account for temporal variations in the scattering behavior of ground targets.

[0007] Therefore, in order to address the shortcomings of the existing technology, the present invention proposes to combine polarimetric SAR information and perform phase optimization based on a sequential estimation framework, thereby achieving near real-time and efficient phase optimization processing while ensuring the accuracy of phase optimization. Summary of the Invention

[0008] The purpose of the present invention is to realize sequential dynamic processing of phase optimization by using multi-polarization SAR data and provide high-quality interferometric phase images in near real time.

[0009] To achieve the above objectives, the present invention proposes a sequential polarization phase optimization method for near real-time DS-InSAR processing, which includes the following processing steps:

[0010] Preparation of environmental files for sequential processing, including dividing the time series multi-polarimetric SAR data into spatial blocks containing multiple mini-time series subsets according to the size of the time series subsets and the number of spatial blocks;

[0011] Sequential processing of environmental file checks, including the current processing of historical data, so that when new SAR data is added, there is no need to repeat the processing of historical SAR data;

[0012] Read SAR data files and preprocess them, including reading SLC files and differential interferogram files according to the calculated spatial block size and temporal subset size, and removing the terrain phase contribution in the SLC phase to reduce the fringe slope in the local window and improve the spatial stability of the local phase;

[0013] Statistically homogeneous pixels (SHP) identification of time series subsets, including SHP homogeneous pixel identification of newly added time series SAR image subsets during sequential processing, and using them as update samples for phase optimization to reduce the impact of dynamic changes in ground target scattering characteristics on the accuracy of coherence matrix estimation;

[0014] SETP-EMI phase optimization processing of the time series subset includes performing TP-EMI-based phase optimization processing based on the read multi-polarization SAR data and the identified homogeneous pixels, and performing data compression processing on the original multi-polarization SLC data after the phase optimization is completed; the data compression processing is implemented by normalizing the optimized SLC phase and using it as the basis for linear transformation to compress the original multi-polarization SLC data;

[0015] Phase-optimized sequential processing and chaining of temporal subsets;

[0016] The total time series of spatial blocks is optimized for phase stitching and the interferometric phase map is output.

[0017] Furthermore, the implementation of the sequential processing environment file preparation is as follows:

[0018] Read the file name of the time series radar image in the SAR file path to determine the total number of images in the file path; use the size of the split time series subset as the number of scenes to calculate the number of mini-subsets;

[0019] Read the first image in the SAR file path to obtain the image row and column number range; calculate the row and column number range of each spatial block based on the number of spatial blocks defined by the user; while calculating the row and column number range of the spatial block, based on the window size of the statistical homogeneous pixel (SHP) recognition set by the user, add a buffer of half the window size to the spatial block to ensure that each spatial block has overlap at the edge;

[0020] According to the number of spatial blocks, corresponding subfolders are created for storing and sequentially processing the block data results; the spatial block list file is saved in the processing file path.

[0021] Furthermore, the statistically homogeneous pixels SHP recognition of the time series subset is implemented as follows:

[0022] After reading and preprocessing the SAR data, read the intensity data of the SAR dataset;

[0023] Homogeneous pixels are identified on the intensity data according to the user-defined SHP homogeneous pixel identification window.

[0024] Moreover, when the original multi-polarization SLC data is compressed after the phase optimization is completed, the normalized maximum likelihood estimated phase is used as the basis of the linear transformation to convert the multi-polarization SLC complex signal from the original high-dimensional data to low-dimensional data. The implementation method is as follows:

[0025] Assume that there are N SLC images with h polarization channels; the number of time series subsets is τ, and each time series subset contains ξ images; at the same time, there are l statistically homogeneous pixels SHP identified in the local window; linear transformation is used to transform the high-dimensional space data Z into ξ·l·h Projection to low-dimensional space Z η·l·h , where Z is the temporal subset matrix of a pixel, ξ·l·h is the dimension of the temporal subset matrix of the pixel, ξ is the time dimension, which represents the number of temporal images, l is the spatial dimension, which represents the number of homogeneous pixels in the SHP area, h is the polarization dimension, which represents the number of polarization channels, and η represents the time dimension of the matrix after compression; the basis of its linear transformation is B = {v1; ...; v ξ}, where v1;…;v ξ Refers to the elements of the transformation basis, there are ξ;

[0026] The complex Pauli niche time series scattering vectors formed by each polarization channel need to be compressed separately. For full polarization data, the three Pauli scattering vectors exist in three compressed low-rank subspaces, h = 3; the ξ time series maximum likelihood estimation signals u estimated by the phase optimization processing based on TP-EMI are EMI Normalize and get the basis B of linear transformation, and transform the original multi-polarization SAR signal of the first time series subset into the following:

[0027]

[0028] in, The SLC after compression of the first timing subset is It is the original timing multipolarization SLC.

[0029] On the other hand, the present invention also provides a sequential polarization phase optimization system for near real-time DS-InSAR processing, which is used to implement the above-mentioned sequential polarization phase optimization method for near real-time DS-InSAR processing.

[0030] Furthermore, the following modules are included,

[0031] The first module is the preparation of environmental files for sequential processing, which includes dividing the time series multi-polarimetric SAR data into spatial blocks containing multiple mini time series subsets according to the size of the time series subsets and the number of spatial blocks;

[0032] The second module is used for environmental file checking for sequential processing, including the processing data check of historical data, so that there is no need to repeat the processing of historical SAR data when new SAR data is added;

[0033] The third module is used to read SAR data files and preprocess them, including reading SLC files and differential interferogram files according to the calculated spatial block size and temporal subset size, and removing the terrain phase contribution in the SLC phase to reduce the fringe slope in the local window and improve the spatial stability of the local phase;

[0034] The fourth module is used for statistically homogeneous pixel (SHP) identification of time-series subsets. This includes SHP homogeneous pixel identification of newly added time-series SAR image subsets during sequential processing and using them as updated samples for phase optimization to reduce the impact of dynamic changes in ground target scattering characteristics on the accuracy of coherence matrix estimation.

[0035] The fifth module is used for SETP-EMI phase optimization processing of the time series subset, including performing TP-EMI-based phase optimization processing based on the read multi-polarization SAR data and the identified homogeneous pixels, and compressing the original multi-polarization SLC data after the phase optimization is completed. The data compression processing is implemented by normalizing the optimized SLC phase and using it as the basis for linear transformation to compress the original multi-polarization SLC data.

[0036] The sixth module is used for sequential processing and linking of phase-optimized timing subsets;

[0037] The seventh module is used to optimize the phase stitching of the total time series of spatial blocks and output the interference phase map.

[0038] Alternatively, the system comprises a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the sequential polarization phase optimization method for near real-time DS-InSAR processing as described above.

[0039] Alternatively, the method comprises a readable storage medium having a computer program stored thereon, and when the computer program is executed, the method implements the sequential polarization phase optimization method for near real-time DS-InSAR processing as described above.

[0040] To address the shortcomings of existing technologies, this paper proposes a SETP-EMI method, combining a sequential estimator (SE) with a total power (TP) polarization coherence matrix approach and using an eigendecomposition-based maximum likelihood estimator (EMI) for the solution. Furthermore, the phase optimization process considers the dynamic changes in the scattering characteristics of ground targets, dynamically identifying statistically homogeneous pixels (SHPs) during sequential processing.

[0041] The solution of the present invention is simple and convenient to implement and has strong practicality. It solves the problems of low practicality and inconvenience in actual application existing in related technologies, can improve user experience, and has important market value. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 4 is a flowchart of the DS phase optimization process according to an embodiment of the present invention.

[0043] Figure 2 1 is a phase optimization processing result diagram of an embodiment of the present invention, wherein the first row is the original interferogram (a1) and the interferogram after EMI (b1), TPEMI (c1) and SETP-EMI (d1) optimization processing, the second row (a2) to (d2) are the average spatial coherence SPC estimated from all interferograms, which are the original interferogram (a2) and the average spatial coherence SPC of the interferogram after EMI (b2), TPEMI (c2) and SETP-EMI (d2) optimization processing, respectively, the third row (a3) ​​to (d3) are the temporal coherence TPC estimated from all interferograms, which are the original interferogram (a3) ​​and the temporal coherence TPC of the interferogram after EMI (b3), TPEMI (c3) and SETP-EMI (d3) optimization processing, respectively. DETAILED DESCRIPTION

[0044] The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0045] In this embodiment of the present invention, a sequential polarization phase optimization method for near-real-time DS-InSAR processing is proposed. The file to be processed is a registered and generated multi-polarimetric SAR dataset containing noise. The SAR data here refers to the first-level product of radar imagery: single-look complex data (SLC). Multi-polarization time-series single-primary image differential interferograms and multi-polarization time-series SAR intensity maps can be obtained through registration, interferometry, and terrain phase removal using various SAR interferometry processing software such as GAMMA and SNAP. See [1] for more information. Figure 1The polarization sequential DS phase optimization process of the embodiment includes the following steps:

[0046] Step 1: Prepare the environment file for sequential processing. Divide the time series multi-polarization single-look complex (SLC) data into spatial blocks containing multiple mini-time series subsets according to the user-defined time series subset size (number of scenes) and the number of spatial blocks (number of blocks in range and azimuth).

[0047] The purpose of spatially partitioning the dataset is to reduce memory consumption, thereby enabling the processing of full-scene SAR images. The purpose of temporal partitioning of time-series images is, on the one hand, to further reduce memory consumption, and on the other hand, to enable the update and identification of homogeneous pixels and improve the accuracy and efficiency of phase optimization. At the same time, it facilitates the update processing of new data without having to reprocess historical data.

[0048] The preferred implementation method of the embodiment includes the following sub-steps:

[0049] (1) Read the file name of the time series radar image in the SAR file path to determine the total number of images in the file path. According to the user-defined time series subset size (number of scenes), calculate the number of mini-subsets: the total number of scenes divided by the number of scenes in the time series subset, and round up.

[0050] (2) Read the first image in the SAR file path to obtain the image row and column number range. Calculate the row and column number range of each spatial block based on the number of spatial blocks defined by the user (the number of blocks in the range and azimuth directions). While calculating the row and column number range of the spatial blocks, add a buffer of half the window size to the spatial blocks based on the window size of the statistically homogeneous pixel (SHP) recognition set by the user, so that each spatial block has overlap at the edge. This processing can further ensure the consistency of spatial stitching after phase optimization.

[0051] (3) Create corresponding subfolders according to the number of spatial blocks for storing and sequentially processing the block data results.

[0052] Finally, save the spatial block list file in the processed file path.

[0053] Step 2: Check the environment files for sequential processing.

[0054] This step checks the historical data for processed data, eliminating the need to reprocess the historical SAR data when adding new SAR data. (These files are the processing results of the SETP-EMI method of the present invention. This allows only compressed SLCs to be read when adding new SLC data, avoiding reprocessing the uncompressed time-series historical SLC data and ensuring smooth dynamic update processing.)

[0055] The specific implementation method preferably adopted in the embodiment includes the following sub-steps:

[0056] (1) Traverse each subfolder according to the space block list file and enter the subfolder path.

[0057] (2) Check the statistical homogeneous pixel (SHP) identification result file and read it if it exists.

[0058] (3) Check the compressed single-view complex (SLC) data result file and read it if it exists.

[0059] (4) Check the file recording the number of processed files and read the file if it exists.

[0060] Step 3: Read SAR data files and preprocess.

[0061] This step reads the SLC file and the differential interferogram file based on the calculated spatial block size and temporal subset size and performs preprocessing. Preprocessing involves removing the topographic phase contribution from the SLC phase, reducing the interference fringe density within the local window, thereby improving the spatial stationarity of the local phase and ultimately increasing the accuracy of phase optimization.

[0062] The specific implementation method preferably adopted in the embodiment includes the following sub-steps:

[0063] (1) After checking and reading the environment file in the subfolder path in step 2, determine the serial number corresponding to the time series SAR file to be read based on the number of processed files, the time series subset size setting, and the number of mini-subsets.

[0064] (2) Read the SLC file and differential interference pattern file according to the serial number and the spatial block size corresponding to the subfolder.

[0065] (3) Perform interference processing on the SLC file to obtain the interference pattern file.

[0066] (4) Differentiate the differential interferogram and the interferogram file to obtain the terrain phase file.

[0067] (5) Differentiate the SLC file and the terrain phase file to obtain the SLC file without terrain phase, completing the preprocessing.

[0068] Step 4: Identification of statistically homogeneous pixels (SHP) of the time series subset.

[0069] Statistically homogeneous pixels (SHP) are identified for time series subsets. The specific SHP identification method is not within the scope of this invention. The difference is that this invention performs SHP identification on time series SAR image subsets. Its characteristic is that during sequential processing, SHP homogeneous pixels are identified for newly added time series SAR image subsets and used as updated samples for phase optimization. This can reduce the impact of dynamic changes in the scattering characteristics of ground targets on the accuracy of coherence matrix estimation. Continuously performing SHP updates on time series subsets is one of the characteristics of the SETP-EMI method.

[0070] The preferred implementation method of identifying homogeneous pixels in a time series subset includes the following two sub-steps:

[0071] (1) After reading and preprocessing the SAR data, read the intensity data of the SAR data set.

[0072] (2) Perform SHP recognition on the intensity data according to the user-defined SHP recognition window.

[0073] Step 5: SETP-EMI phase optimization processing of the timing subset.

[0074] Based on the read multi-polarimetric SAR data and identified homogeneous pixels, a phase optimization process using the eigendecomposition maximum likelihood estimator (EMI) based on the total power coherence matrix (TP) is performed, namely the TP-EMI method. However, after the phase optimization, the original multi-polarimetric SLC data is compressed, which is one of the characteristics of the SETP-EMI method.

[0075] The specific implementation method preferably adopted in the embodiment includes the following sub-steps: Sub-steps 1 to 5 are the same as the TP-EMI method and do not belong to the innovative content of the present invention.

[0076] 1. According to the polarization channel set by the user, the multi-polarization data read is combined and converted into a complex bubble niche scattering vector

[0077] k Pol . Reciprocity of backscattering is assumed here.

[0078] When the polarization channel is fully polarized, the complex bubble niche scattering vector is expressed as:

[0079]

[0080] When the polarization channel xx is dual-polarization or cross-polarization, it is expressed as: Here, xx can be replaced by VV or HH.

[0081] k Pol =[S xx 2S HV ] T

[0082] Here the intermediate variable S HH+VV =S HH +S VV ,S HH-VV =S HH -S VV .S HH , S VV and S VH are the complex scattering vectors of the HH polarization channel, VV polarization channel, and VH polarization channel, that is, the complex values ​​of the SLC pixels under different polarization channels, and T represents the matrix transpose.

[0083] 2. Extraction of time-series polarization complex scattering vectors of statistically homogeneous pixels SHP. The embodiment traverses each pixel and extracts the complex bubble niche scattering vector of each pixel corresponding to the SHP based on the SHP identification results. The result is the time-series complex scattering vector k in the statistically homogeneous pixel (SHP) area. TSUn , k TSIn =[s 1 ,…,s N ] T The complex scattering vector s including the N time series polarization SAR data collected 1 ,…,s N If there are l homogeneous pixels in the SHP region, there are l time-series complex scattering vectors k TSIn .

[0084] 3. Temporal Interference Coherence Matrix T of Complex Bubble Niche TSIn The embodiment generates the l time series scattering vectors k under statistical homogeneous pixels (SHP) and each complex bubble niche. TSIn Generate the interference coherence matrix of each complex bubble niche separately.

[0085] Using the N collected time series polarimetric SAR data, the SAR data of each complex bubble niche can generate two-to-two interferences to form a maximum of N(N-1) / 2 interference pairs. The time series interferometric coherence matrix of each pixel in each polarimetric complex bubble niche SAR data contains the measurement information of all possible interference pair combinations. This matrix can be obtained by identifying the l time series complex scattering vectors k in the statistically homogeneous pixel SHP region. TSIn get:

[0086]

[0087] Where T TSIn is the N×N dimensional interference coherence matrix. TSIn N time series signals representing the scattering vector of each polarization complex bubble niche. represents the complex conjugate transpose, and <·> denotes multilook processing, which is obtained by summing the complex signals after conjugate multiplication within the SHP spatial region. The scattering vector of the complex bubble niche consists of three or two polarization channels, so each channel has its own corresponding coherence matrix. Multilook processing can mitigate speckle noise in the scattering vector and improve the maximum accuracy of single master image (SM) phase optimization.

[0088] 4. Time series total power (TP) coherence matrix T TSTP For multi-polarization data, using all time series interference coherence matrices from multiple polarization complex bubble niches as statistical samples can obtain more accurate coherence estimation, thereby further optimizing the SM phase. The time series total power (TSTP) coherence matrix T can be formed by simply adding the coherence matrices of the complex bubble niches of each polarization channel. TSTP ,as follows:

[0089]

[0090] Where Pol refers to the complex bubble niche of polarization SLC data. For SAR data with h polarization channels, h complex bubble niche scattering vectors can be generated according to sub-step 1. Here, the temporal interferometric coherence matrix of each complex bubble niche is T TSIn The method for constructing the TP coherence matrix has been proposed and is not part of the present invention.

[0091] 5. Use the EMI estimator to process the TP coherence matrix and obtain the optimized SLC phase.

[0092] Based on the phase triangulation (PTA) triangular phase closure theory, the total power coherence matrix T of the time series is used. TSTP Instead of the single polarization coherence matrix, the maximum likelihood estimate (MLE) of the single main image (SM) can be obtained as follows under the complex Wishart distribution statistics:

[0093]

[0094] here represents the Hadamard product, θ=[θ1,θ2,…,θ N ] T is the phase of the single master image (SM) of the N-scene time series. γ represents the real complex coherence matrix, |·| represents the modulo operation of the complex matrix, |γ| -1 Represents the inverse of the real complex coherence matrix. Since the real coherence matrix γ is unknown, the estimated coherence matrix T is usually used instead of the real coherence matrix γ. In the present invention, the time series total power (TSTP) coherence matrix T is used. TSTPIn order to efficiently solve the function, the maximum likelihood estimation method based on eigendecomposition, i.e. EMI, is used to obtain the minimum eigenvector u EMI As the estimation result, the above formula is rewritten as:

[0095]

[0096] Here u is the matrix after Hadamard product The smallest eigenvector of .

[0097] 6. SLC data compression: Specifically, the optimized SLC phase is normalized and used as the basis for linear transformation to compress the original multi-polarization SLC data.

[0098] Assume that there are N SLC images with h polarization channels. The number of time series subsets is τ, and each time series subset contains ξ images. At the same time, there are l statistically homogeneous pixels (SHP) identified in the local window. The high-dimensional space data Z can be transformed into ξ·l·h Projection to low-dimensional space Z η·l·h Where Z is the temporal subset matrix of a pixel, ξ·l·h is the dimension of the temporal subset matrix of the pixel. ξ is the time dimension (the number of temporal images), l is the spatial dimension (the number of homogeneous pixels in the SHP area), h is the polarization dimension (the number of polarization channels), and η represents the time dimension of the matrix after compression. The basis of its linear transformation is B = {v1; ...; v ξ}, where v1;…;v ξ Refers to the elements of the transformation basis, and there are ξ. Assuming that there is only a single dominant scattering mechanism in the SHP region, the compression dimension η can be set to 1. Unlike the compression of single-polarization data, since the present invention uses multiple polarization channels, the complex Pauli time-series scattering vectors formed by each polarization channel need to be compressed separately. For full polarization data, there are three compressed low-rank subspaces (h = 3) for the three Pauli scattering vectors. The ξ time-series maximum likelihood estimation signals u estimated in sub-step 5 are EMI Normalize it and get As the basis of linear transformation. EMI Optimize the phase of the ξ time series for maximum likelihood estimation. exp(j*u EMI ) means u EMI The phase is converted to a complex number, j is the imaginary unit. ‖exp(j*u EMI )‖ represents the difference between exp(j*u EMI ) Take the norm. Using the basis B of linear transformation, the original multi-polarization SAR signal of the first time series subset can be transformed into the following:

[0099]

[0100] in The SLC after compression of the first timing subset is It is the original timing multipolarization SLC.

[0101] 7. Save the optimized time-series single master image (SM) phase, compressed SLC data, identified homogeneous pixel set SHP file, and the number of processed SAR files to a folder.

[0102] Step 6: Phase-optimized sequential processing and linking of timing subsets.

[0103] After completing phase optimization and compression for one time-series SAR subset, sequential recursive processing is required for subsequent time-series SAR subsets without processing historical data. Sequential recursive processing combines the compressed SLC with the newly added time-series subsets for phase optimization. During phase linking, all compressed SLC data is phase-optimized to obtain a baseline phase, completing the phase linking of all time-series subsets.

[0104] The specific implementation method preferably adopted in the embodiment includes the following sub-steps:

[0105] (1) Phase optimization of subsequent time series subsets. Repeat steps 2 to 5, except that the compressed SLC generated by the previous time series subset is read and used as the first scene image of the subsequent time series subset. τ ,

[0106] (τ>1), the compressed timing subset SLC Equal combinations form new time series subsets as follows:

[0107]

[0108] in For the compressed SLC of the first timing subset, For the compressed SLC of the second timing subset, For the

[0109] The compressed SLC of τ-1 time series subsets. The difference is that Z τ is the uncompressed SLC of the newly added τth timing subset.

[0110] In this way, the compressed SLC and the newly added SLC can be combined into a new timing subset New time series subset It is used as the SAR dataset to be input into the processing of steps 2 to 5.

[0111] (2) Phase optimization of the last time series subset. If the size of the time series subset is less than the number of time series subset images defined by the user, the SHP result of the previous time series subset is used to ensure the robustness of SHP recognition. Finally, the phase optimization and SLC compression of the last time series subset are completed. After the phase optimization is finally completed, only the time series subset Z is used. τ Compress and generate compressed SLC data of the new timing subset

[0112] (3) Compressed SLC phase benchmark link processing. After completing the phase optimization and SLC compression of all time series subsets, EMI estimation is performed again on all compressed SLCs to obtain the benchmark phase. The SHP homogeneous pixel set here is the homogeneous pixel identification result of the last time series subset. The phases of the optimized results of multiple time series subsets are then benchmarked to obtain the optimized phase of the total time series.

[0113] For τ compressed timing subsets SLC, the optimized phase is calculated by the EMI method, that is, the reference offset phase between each timing subset The calculation formula is as follows:

[0114]

[0115] Where T TSTP_sub Represents a compressed time series subset The TP coherence matrix constructed as the SLC data set. The reference offset phase is obtained Then, the corrected phase can be obtained by adding the reference offset phase to the original optimized phase. The correction phase for the i-th time series subset is written as:

[0116]

[0117] here is the optimized phase of the maximum likelihood estimate of the i-th time series subset, is the reference offset phase estimated by the maximum likelihood of the compressed time series subset, is the correction phase of the i-th time series subset. In this way, the optimized phase of the complete time series can be obtained.

[0118] Step 7: Phase stitching of the total time series of spatial blocks. After completing the phase optimization, SLC compression, and reference linking of all spatial blocks, the results are spatially stitched, and the interferometric phase map is output and saved as a file. The specific implementation method includes the following sub-steps:

[0119] (1) Read the spatial block list file and determine the block cropping row and column number range and the mosaic row and column number range.

[0120] (2) Check the phase optimization results in the result folder and determine the file name of the newly added single principal interferogram.

[0121] (3) Mosaic the phase optimization results of the spatial blocks, save the file with the file name of the newly added single main interference pattern, and complete all processing.

[0122] The present invention primarily utilizes polarization information to increase the number of phase optimization samples, increases the homogeneity of phase optimization samples by dynamically updating homogeneous sets, and increases phase optimization efficiency and provides dynamic update capabilities through a sequential subset recursive method, ultimately achieving efficient and high-precision sequential phase optimization processing. In specific implementations, the method proposed by the technical solution of the present invention can be automatically run by those skilled in the art using computer software technology, such as using Matlab to program a SETP-EMI processing system. System devices implementing the method, such as a computer-readable storage medium storing a computer program corresponding to the technical solution of the present invention and a computer device that runs the corresponding computer program, should also fall within the scope of protection of the present invention.

[0123] In some possible embodiments, a sequential polarization phase optimization system for near real-time DS-InSAR processing is provided, comprising the following modules:

[0124] The first module is the preparation of environmental files for sequential processing, which includes dividing the time series multi-polarimetric SAR data into spatial blocks containing multiple mini time series subsets according to the size of the time series subsets and the number of spatial blocks;

[0125] The second module is used for environmental file checking for sequential processing, including the processing data check of historical data, so that there is no need to repeat the processing of historical SAR data when new SAR data is added;

[0126] The third module is used to read SAR data files and preprocess them, including reading SLC files and differential interferogram files according to the calculated spatial block size and temporal subset size, and removing the terrain phase contribution in the SLC phase to reduce the fringe slope in the local window and improve the spatial stability of the local phase;

[0127] The fourth module is used for statistically homogeneous pixel (SHP) identification of time-series subsets. This includes SHP homogeneous pixel identification of newly added time-series SAR image subsets during sequential processing and using them as updated samples for phase optimization to reduce the impact of dynamic changes in ground target scattering characteristics on the accuracy of coherence matrix estimation.

[0128] The fifth module is used for SETP-EMI phase optimization processing of the time series subset, including performing TP-EMI-based phase optimization processing based on the read multi-polarization SAR data and the identified homogeneous pixels, and compressing the original multi-polarization SLC data after the phase optimization is completed. The data compression processing is implemented by normalizing the optimized SLC phase and using it as the basis for linear transformation to compress the original multi-polarization SLC data.

[0129] The sixth module is used for sequential processing and linking of phase-optimized timing subsets;

[0130] The seventh module is used to optimize the phase stitching of the total time series of spatial blocks and output the interference phase map.

[0131] During specific implementation, the system can also be divided into four modules for implementation, namely user processing parameter setting module, environment file preparation module, phase optimization processing module, and result splicing and storage module.

[0132] This example uses dual-polarization data from the Sentine-1 SAR satellite as an example to illustrate this. First, prepare the files to be processed: the single principal differential interferogram (diff) generated using GAMMA software, the registered SLC data (slc), and the time-series intensity map (mli). These files can also be processed using other InSAR processing software such as ISCE and SNAP.

[0133] Module 1: User Processing Parameter Setting The module contains various parameters that need to be defined by the user, as follows:

[0134] Set the working path for storing files to be processed; set the polarization channel to select the SAR polarization channel data to be used. In this embodiment, set channels = ["vv"; "vh"], which means dual polarization channels vv and vh are used. Set the row number to determine the file size to be read. The reference master image number is the reference master image number of the single main differential interferogram. The image cropping row and column numbers are used to crop the original image. The number of spatial blocks (the number in azimuth and the number in range) is set in this embodiment as paz = 2 in azimuth and prg = 4 in range, resulting in a final partitioning of paz * prg = 8 spatial sub-blocks. The time-series SAR subset size (number of frames) interval = 10, meaning every 10 SAR images are considered a time-series subset. The homogeneous point identification parameter (local window size) is set in this embodiment as hW_w = 5 for half the window width and hW_l = 3 for half the window height. Set the number of CPU parallel cores to 12, meaning that 12 threads are processed in parallel when traversing each pixel for phase optimization, increasing processing efficiency.

[0135] Module 2: Environmental file preparation module, which implements the function of step 1 in the invention.

[0136] Create a blank result folder to store the final phase-optimized differential interferogram.

[0137] Read the file name in the SAR dataset path and determine the total number of images in the file path. Calculate the number of mini-subsets based on the user-defined time series subset size (number of scenes).

[0138] The time series multi-polarimetric SAR data are divided into spatial blocks containing multiple mini-subsets according to the user-defined time series subset size (number of scenes) and the number of spatial blocks (number of blocks in the range and azimuth directions), and a block row and column number file is generated.

[0139] Create corresponding subfolders based on the number of spatial blocks for storing and sequentially processing the block data results. Finally, save the spatial block list file in the processing file path.

[0140] Module 3: SETP-EMI phase optimization processing module. Implements the functions of steps 2 to 6 in the invention.

[0141] Sequential processing environment file check. This checks the processed historical data so that when new SAR data is added, there is no need to reprocess the historical SAR data. It includes the following sub-steps:

[0142] Traverse each subfolder according to the space block list file and enter the subfolder path.

[0143] Check the homogeneous pixel identification result file and read it if it exists.

[0144] Checks for a compressed SLC results file and reads it if it exists.

[0145] Check the file that records the number of processed files and read it if it exists.

[0146] Read SAR data files and preprocess them. Read SLC files and differential interferogram files based on the calculated spatial block size and temporal subset size, and perform preprocessing. This includes the following sub-steps:

[0147] The serial number corresponding to the time series SAR file to be read is determined based on the number of processed files, the time series subset size setting, and the number of mini-subsets.

[0148] Read the SLC file and differential interferogram file according to the sequence number and the spatial block size corresponding to the subfolder.

[0149] Perform interference processing on the SLC file to obtain the interference pattern file.

[0150] The differential interferogram and the interferogram file are differentiated to obtain the terrain phase file.

[0151] The SLC file and the terrain phase file are differentiated to obtain the SLC file without the terrain phase.

[0152] SHP homogeneous pixel identification of time series subsets includes the following two sub-steps:

[0153] After reading and preprocessing the SAR data, read the intensity data of the SAR dataset.

[0154] Homogeneous pixels are identified on the intensity data according to the user-defined SHP homogeneous pixel identification window. The open-source FaSHP homogeneous pixel identification code is used to implement SHP homogeneous pixel identification. The specific SHP homogeneous pixel identification method is beyond the scope of this invention.

[0155] SETP-EMI phase optimization of the timing subset includes the following sub-steps:

[0156] The multi-polarization data is combined according to the user-defined polarization channel and converted into a complex bubble niche scattering vector k_Pol. The time-series polarization complex scattering vector of homogeneous pixels is extracted. CPU parallel processing is used to traverse each pixel and extract the corresponding complex bubble niche scattering vector of each homogeneous pixel based on the SHP homogeneous pixel identification results.

[0157] Generation of temporal interferometric coherence matrix. The interferometric coherence matrix is ​​generated from a homogeneous pixel set and the temporal scattering vectors under each complex bubble niche.

[0158] Generation of the Total Power Coherence Matrix (TP): The Time Series Total Power (TSTP) coherence matrix can be formed by simply adding the coherence matrices of each polarization channel.

[0159] The EMI estimator is used to process the TP coherence matrix and obtain the optimized SLC phase.

[0160] Based on the Phase Linking triangular phase closure theory, the time series total power (TSTP) coherence matrix is ​​used instead of the single-polarization coherence matrix. Eigendecomposition-based maximum likelihood estimation (EMI) is employed to obtain the minimum eigenvector as the estimation result. EMI phase optimization methods are not within the scope of this invention.

[0161] Compression of the optimized SLC: The optimized SLC phase is normalized and used as the basis for linear transformation to compress the original multi-polarization SLC data.

[0162] Save the optimized time-series single master image (SM) phase, compressed SLC data, identified homogeneous pixel set SHP file, and the number of processed SAR files to a folder.

[0163] Sequential processing and chaining of phase optimization of timing subsets. This includes the following sub-steps:

[0164] The phase optimization process of the subsequent time series subset is repeated. The process (2) to (5) is repeated for the subsequent time series subset, except that the compressed SLC generated by the previous time series subset is read and used as the first scene image of the subsequent time series subset.

[0165] After the phase optimization of the subsequent timing subset is finally completed, only the timing subset is compressed to generate compressed SLC data of the new timing subset.

[0166] Phase optimization processing of the last time series subset. If the time series subset size is less than the number of time series subset images defined by the user, the SHP result of the previous time series subset is used to ensure the robustness of SHP recognition. Finally, the phase optimization processing and SLC compression of the last time series subset are completed.

[0167] Compressed SLC phase benchmark linking. After completing phase optimization and SLC compression for all time series subsets, EMI estimation is performed again on all compressed SLCs to obtain a benchmark phase. The SHP homogeneous pixel set here is the homogeneous pixel identification result of the last time series subset. The phases of the optimized results of multiple time series subsets are then benchmarked to obtain the optimized phase of the complete time series.

[0168] Module 4: Result splicing and storage module. After completing phase optimization, SLC compression, and benchmark linking for all spatial blocks, the results are spatially spliced ​​and saved as a file. The specific implementation method includes the following sub-steps:

[0169] (1) Read the spatial block list file and determine the block cropping row and column number range and the mosaic row and column number range.

[0170] (2) Check the phase optimization results in the result folder and determine the file name of the newly added single principal interferogram.

[0171] (3) Mosaic the phase optimization results of the spatial blocks, save the file with the file name of the newly added single main interference pattern, and complete all processing.

[0172] The present invention utilizes multi-polarization channel information to increase maximum likelihood estimation samples and combines a sequential estimator to dynamically update the phase optimization result of a data set, which not only reduces memory consumption but also improves the accuracy and efficiency of phase optimization.

[0173] The embodiment area is a landslide deformation area on the bank of a reservoir. The comparison of the results of the original interference pattern and other traditional methods, such as the maximum likelihood estimator (EMI) based on eigendecomposition and the maximum likelihood estimator (TP-EMI) based on the total power coherence matrix, is shown in the attached figure. Figure 2 . Figure 2 The first row is the original interferogram (a1) and the interferogram after optimization by EMI (b1), TPEMI (c1) and SETP-EMI (d1). The second row (a2) to (d2) are the average spatial coherence SPC estimated from all interferograms, which are the original interferogram (a2) and the average spatial coherence (SPC) of the interferogram after optimization by EMI (b2), TPEMI (c2) and SETP-EMI (d2). The third row (a3) ​​to (d3) are the temporal coherence (TPC) estimated from all interferograms, which are the original interferogram (a3) ​​and the temporal coherence TPC of the interferogram after optimization by EMI (b3), TPEMI (c3) and SETP-EMI (d3). The spatial coherence SPC and temporal coherence TPC here are widely used evaluation indicators for evaluating the quality of interferograms. It can be seen that the interferogram obtained by the invented SETP-EMI method has less noise. Appendix Figure 2 The interference fringes in the landslide deformation area marked by the red circle are clearer. The spatial coherence SPC and temporal coherence TPC indicators are used to evaluate the interference pattern. It is found that the coherence after optimization of the SETP-EMI method is higher.

[0174] In some possible embodiments, a sequential polarization phase optimization system for near-real-time DS-InSAR processing is provided, including a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute a sequential polarization phase optimization method for near-real-time DS-InSAR processing as described above.

[0175] In some possible embodiments, a sequential polarization phase optimization system for near-real-time DS-InSAR processing is provided, including a readable storage medium having a computer program stored thereon. When the computer program is executed, the sequential polarization phase optimization method for near-real-time DS-InSAR processing as described above is implemented.

[0176] The specific embodiments described herein are merely illustrative of the spirit of the present invention. Persons skilled in the art may make various modifications, additions, or substitutions to the described specific embodiments without departing from the spirit of the present invention or exceeding the scope of the appended claims.

Claims

1. A sequential polarimetric phase optimization method for near real-time DS-InSAR processing, characterized in that: The following processing steps are included: Preparation of environmental files for sequential processing, including dividing the time series multi-polarimetric SAR data into spatial blocks containing multiple mini-time series subsets according to the size of the time series subsets and the number of spatial blocks; Sequential processing of environmental file checks, including checking of processed historical data, so that historical SAR data does not need to be reprocessed when new SAR data is added; Read SAR data files and preprocess them, including reading SLC files and differential interferogram files according to the calculated spatial block size and temporal subset size, and removing the terrain phase contribution in the SLC phase to reduce the fringe slope in the local window and improve the spatial stability of the local phase; Statistically homogeneous pixels (SHP) identification of time series subsets, including SHP homogeneous pixel identification of newly added time series SAR image subsets during sequential processing, and using them as update samples for phase optimization to reduce the impact of dynamic changes in ground target scattering characteristics on the accuracy of coherence matrix estimation; SETP-EMI phase optimization processing of the time series subset, including TP-EMI-based phase optimization processing based on the read multi-polarization SAR data and the identified homogeneous pixels, and data compression processing of the original multi-polarization SLC data after the phase optimization is completed; The data compression process is implemented by normalizing the optimized SLC phase and using it as a basis for linear transformation to perform data compression on the original multi-polarization SLC; Phase-optimized sequential processing and chaining of temporal subsets; The total time series of spatial blocks is optimized for phase stitching and the interferometric phase map is output.

2. The sequential polarization phase optimization method for near real-time DS-InSAR processing according to claim 1, characterized in that: The implementation of the sequential processing environment file preparation is as follows: Read the file name of the time series radar image in the SAR file path to determine the total number of images in the file path; use the size of the split time series subset as the number of scenes to calculate the number of mini-subsets; Read the first image in the SAR file path to obtain the image row and column number range; calculate the row and column number range of each spatial block based on the number of spatial blocks defined by the user; while calculating the row and column number range of the spatial block, based on the window size of the statistical homogeneous pixel (SHP) recognition set by the user, add a buffer of half the window size to the spatial block to ensure that each spatial block has overlap at the edge; Create corresponding subfolders according to the number of spatial blocks for storing and sequentially processing the block data results; Save the spatial chunk list file in the processed file path.

3. The sequential polarization phase optimization method for near real-time DS-InSAR processing according to claim 1, characterized in that: The statistical homogeneous pixel SHP recognition of the time series subset is implemented as follows: After reading and preprocessing the SAR data, read the intensity data of the SAR dataset; Homogeneous pixels are identified on the intensity data according to the user-defined SHP homogeneous pixel identification window.

4. The sequential polarization phase optimization method for near real-time DS-InSAR processing according to claim 1, characterized in that: When the original multi-polarization SLC data is compressed after phase optimization, the normalized maximum likelihood estimated phase is used as the basis of linear transformation to convert the multi-polarization SLC complex signal from the original high-dimensional data to low-dimensional data. The implementation method is as follows: Assume that there are N SLC images with h polarization channels; the number of time series subsets is τ, and each time series subset contains ξ images; at the same time, there are l statistically homogeneous pixels SHP identified in the local window; linear transformation is used to transform the high-dimensional space data Z into ξ·l·h Projection to low-dimensional space Z η·l·h , where Z is the temporal subset matrix of a pixel, ξ·l·h is the dimension of the temporal subset matrix of the pixel, ξ is the time dimension, which represents the number of temporal images, l is the spatial dimension, which represents the number of homogeneous pixels in the SHP area, h is the polarization dimension, which represents the number of polarization channels, and η represents the time dimension of the matrix after compression; the basis of its linear transformation is B = {v1; ...; v ξ }, where v1;…;v ξ Refers to the elements of the transformation basis, there are ξ; The complex Pauli niche time series scattering vectors formed by each polarization channel need to be compressed separately. For full polarization data, there are three compressed low-rank subspaces for the three Pauli scattering vectors, h = 3; The maximum likelihood estimation signal u of the time series estimated by the phase optimization processing based on TP-EMI is EMI Normalize and get the basis B of linear transformation, and transform the original multi-polarization SAR signal of the first time series subset into the following: in, The SLC after compression of the first timing subset is It is the original timing multipolarization SLC.

5. A sequential polarimetric phase optimization system for near real-time DS-InSAR processing, characterized by: Used to implement a sequential polarization phase optimization method for near real-time DS-InSAR processing as described in any one of claims 1-4.

6. The sequential polarization phase optimization system for near real-time DS-InSAR processing according to claim 5, characterized in that: Includes the following modules, The first module is the preparation of environmental files for sequential processing, which includes dividing the time series multi-polarimetric SAR data into spatial blocks containing multiple mini time series subsets according to the size of the time series subsets and the number of spatial blocks; The second module is used for environmental file checking for sequential processing, including checking the processed historical data, so that the historical SAR data does not need to be processed repeatedly when new SAR data is added; The third module is used to read SAR data files and preprocess them, including reading SLC files and differential interferogram files according to the calculated spatial block size and temporal subset size, and removing the terrain phase contribution in the SLC phase to reduce the fringe slope in the local window and improve the spatial stability of the local phase; The fourth module is used for statistically homogeneous pixel (SHP) identification of time-series subsets. This includes SHP homogeneous pixel identification of newly added time-series SAR image subsets during sequential processing and using them as updated samples for phase optimization to reduce the impact of dynamic changes in ground target scattering characteristics on the accuracy of coherence matrix estimation. The fifth module is used for SETP-EMI phase optimization processing of the time series subset, including TP-EMI-based phase optimization processing based on the read multi-polarization SAR data and the identified homogeneous pixels, and data compression processing of the original multi-polarization SLC data after the phase optimization is completed; The data compression process is implemented by normalizing the optimized SLC phase and using it as a basis for linear transformation to perform data compression on the original multi-polarization SLC; The sixth module is used for sequential processing and linking of phase-optimized timing subsets; The seventh module is used to optimize the phase stitching of the total time series of spatial blocks and output the interference phase map.

7. The sequential polarization phase optimization system for near real-time DS-InSAR processing according to claim 5, characterized in that: The method comprises a processor and a memory, wherein the memory is used to store program instructions, and the processor is used to call the stored instructions in the memory to execute the sequential polarization phase optimization method for near real-time DS-InSAR processing according to any one of claims 1 to 4.

8. The sequential polarimetric phase optimization system for near real-time DS-InSAR processing according to claim 5, characterized in that: The method comprises a readable storage medium having a computer program stored thereon, wherein when the computer program is executed, the method realizes a sequential polarization phase optimization method for near real-time DS-InSAR processing according to any one of claims 1 to 4.

Images