A multi-scale high-coherent scatterer combined InSAR precise measurement method and system

Abstract

This invention discloses a multi-scale high-coherence scatterer joint InSAR precision measurement method, comprising: 1. Initially constructing an initial layer PS node network; 2. Updating the PS candidate point set with tight to wide constraints, constructing a new high-coherence PS network of the same scale fused with the previous layer PS network; 3. Identifying adjacent homogeneous points of image pixels of non-PS points, and performing phase optimization on image pixels of non-PS points; 4. Updating the DS candidate point set with tight to wide constraints, constructing new high-coherence networks of different scales for PS and DS; 5. Joint inversion of elevation or deformation rate of multi-scale high-coherence scatterers of PS and DS. This invention also discloses a multi-scale high-coherence scatterer joint InSAR precision measurement system. This invention is used to improve the monitoring point density, coverage, and accuracy of traditional PSInSAR and DSInSAR technologies, and can be widely applied in the field of synthetic aperture radar data processing technology.

Description

Technical Field

[0001] This invention relates to the field of synthetic aperture radar data processing technology, and in particular to an InSAR precision measurement method and system that combines multi-scale highly coherent scatterers. Background Technology

[0002] Synthetic Aperture Radar Interferometry (InSAR) technology has become one of the important remote sensing methods for geodesy using space imagery after decades of development. The measurement accuracy of traditional differential InSAR (D-InSAR) is limited by factors such as spatiotemporal decoherence, atmospheric delay, and DEM errors. To overcome these limitations, temporal InSAR technology introduces multi-temporal SAR images in the time dimension and selects pixels with less decoherence impact for phase modeling, thus giving rise to two major technical paths: Permanent Scatterer Interferometry (PS-InSAR) and Distributed Scatterer Interferometry (DS-InSAR).

[0003] PS-InSAR extracts deformation information from point targets (such as buildings, rocks, and artificial corner reflectors) with stable scattering characteristics within a resolution cell, offering advantages such as high signal-to-noise ratio and measurement accuracy down to the millimeter level. However, in non-urban environments such as vegetated areas and bare soil areas, PS points are often sparse or even missing, leading to insufficient spatial sampling. DS-InSAR, on the other hand, identifies statistically homogeneous pixel groups and uses spatial-temporal filtering to improve the phase quality of low signal-to-noise ratio targets, thereby significantly increasing the density of measurement points in sparse PS regions. Multi-scale scatterer InSAR technology, which combines these two types of scatterers with different spatial scales, has become one of the cutting-edge technologies in this field. The main PS-DS joint processing methods currently include SqueeSAR (Ferretti et al. (2011)), JSInSAR (Lv et al. (2014)), and CSI (Wang et al. (2018)). The core of PS-DS joint processing lies in: using the spatial statistical characteristics of DS to preprocess it into PS-like targets, and then incorporating PS and DS into the parameter estimation framework of PS-InSAR. However, the selection of PS and DS targets is extremely sensitive to the screening threshold parameters. The ICOPS technique (Fadhillah et al. (2022)) points out that direct joint temporal InSAR measurements face the problem of insufficient robustness. Therefore, PS and DS targets need to be further optimized and identified before temporal InSAR analysis. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the aforementioned background technology and provide a multi-scale, highly coherent scatterer combined InSAR precision measurement method and system, which can be used to improve the monitoring point density, coverage and accuracy of traditional PSInSAR and DSInSAR technologies.

[0005] This invention provides a multi-scale highly coherent scatterer joint InSAR precision measurement method, comprising the following steps: S1, inputting multi-temporal interferometric SAR image data and initially constructing an initial layer PS node network; S2, setting tight to wide constraints to update the PS candidate point set and constructing a new high-coherence PS network of the same scale fused with the previous layer PS network; S3, identifying adjacent homogeneous points of image pixels that are not PS points and performing phase optimization on the image pixels that are not PS points; S4, setting tight to wide constraints to update the DS candidate point set and constructing a new high-coherence network of different scales of PS and DS that fused the latest layer PS network and the previous layer DS network; S5, jointly inverting the elevation or deformation rate of the multi-scale highly coherent scatterers of PS and DS.

[0006] According to the above scheme, the specific process in step S1 is as follows: For the phase-preprocessed multi-temporal interferometric SAR image data, set the initial amplitude deviation threshold to a small value, and retain the points below the threshold as the initial PS candidate point set; construct a Delaunay triangulation network with neighborhood distance constraints, calculate the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of the network nodes; set a coherence threshold, remove points with low average set coherence coefficient and edges with low set coherence coefficient, retain high coherence points connected by high coherence edges, form a high coherence network of initial layer PS points, and at the same time retain the set coherence coefficient, residual elevation gradient and deformation rate gradient of the high coherence edges connected by the network.

[0007] Furthermore, in step S1, the specific steps are as follows: S11, input the phase-preprocessed multi-temporal SAR image data and calculate the amplitude deviation of the multi-temporal interferometric SAR image; S12, set the initial amplitude deviation threshold to a small value and retain the pixels in the image below the threshold as the initial PS candidate point set; S13, set the distance constraint threshold of the edges connecting the initial PS candidate points and construct a neighborhood distance-constrained Delaunay triangulation network; S14, quickly search for the parameters of the maximizing set coherence coefficient model using CPU parallelism and GPU computing strategies. (1) Obtain the set coherence coefficient of any connected edge in the candidate point set of PS. Residual elevation gradient Average deformation rate gradient And the average set coherence coefficient of the edges connected to any PS candidate point, in the optimization formula (1) For any two adjacent PS candidate points and The phase difference operator, Indicates the first One PS candidate point in Differential interferometric phase in spectral interferometric SAR images and They represent the first The PS candidate point at the th The gradient form of the simulated residual terrain and linear deformation interferometric phase in the interferogram can be expressed as follows: (2), where, and They represent the first Vertical baseline and observation time of each image. and S15: Set a coherence threshold, remove points with low average set coherence coefficients and edges with low set coherence coefficients, retain high coherence points connected by high coherence edges, form a high coherence network of PS nodes in the initial layer, and retain the set coherence coefficient, residual elevation gradient and deformation rate gradient of the high coherence edges connected to the network.

[0008] According to the above scheme, the specific process in step S2 is as follows: gradually relax the amplitude deviation threshold constraint from small to large, and select a new PS candidate point set; combine with the previous layer PS node network, and gradually execute the triangulation network construction in step S1, as well as the calculation steps of the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; set a coherence threshold, and select candidate points that are highly coherently fused with the previous layer PS network from the PS candidate point set of this layer network, as the PS point set of this layer network, while retaining the highly coherent edges connected between the PS point sets of this layer network to update the previous layer PS high coherence network, as well as the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the network; gradually execute this process until the most relaxed amplitude deviation threshold is reached, while retaining the network fused in the last step as a new highly coherent PS network of the same scale, and recording the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the network.

[0009] Furthermore, in step S2, the specific steps are as follows: S21, reset the PS candidate point screening threshold to a smaller value, and select a new set of PS candidate points; S22, combine the previous layer's PS point high coherence network, execute steps S13-S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; S23, set a coherence threshold, and select candidate points from the PS candidate point set of this layer that are highly coherently fused with the previous layer's PS network, as the candidates for this layer's network. The PS point set is used to update the previous layer's PS high coherence network, while retaining the highly coherent edges connecting the PS point sets in this layer. The set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected in the network are also recorded. In step S24, the amplitude deviation threshold constraint is gradually relaxed from small to large, and the PS candidate point set is updated. Steps S22 to S23 are executed repeatedly until the most relaxed amplitude deviation threshold is reached. At the same time, the network fused in the last step is retained as a new highly coherent PS network of the same scale, and the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected in the network are recorded.

[0010] According to the above scheme, the specific process in step S3 is as follows: Based on the temporal similarity metric, the adjacent homogeneous points of the image pixels of the non-PS points are identified, and the number of adjacent homogeneous points of the non-PS points is recorded; based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points, the coherence matrix of the non-PS points is calculated; singular value decomposition is performed on the coherence matrix of the non-PS points, and the phase of the strongest scattering mechanism within the image pixels of the non-PS points is separated using the maximum scattering mechanism separation principle, and this phase is used as the optimized phase, while the phase fitting goodness of the non-PS points is recorded.

[0011] Furthermore, in step S3, the specific steps are as follows: S31, based on the temporal similarity metric, identify the adjacent homogeneous points of the image pixels of the non-PS points and record the number of adjacent homogeneous points of the non-PS points; S32, based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points, calculate the coherence matrix of the non-PS points; S33, perform singular value decomposition on the coherence matrix of the non-PS points, and use the maximum scattering mechanism separation principle to separate the phase of the strongest scattering mechanism within the image pixels of the non-PS points, use it as the optimized phase, and record the phase fitting goodness of the non-PS points.

[0012] According to the above scheme, the specific process in step S4 is as follows: Set a constraint on the number of homogeneous points, gradually relax the goodness-of-fit threshold constraint from large to small, and select a new set of candidate DS points excluding the latest PS point set; combine the latest layer PS network with the previous layer DS network, and progressively execute the triangulation construction in step S1, the calculation steps for the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; set a coherence threshold, and select points from the DS candidate point set that are consistent with the latest layer PS network and the previous layer DS network. Candidate points for coherent fusion are used as the DS point set of the network layer. At the same time, the highly coherent edges connecting the PS and DS point sets of the network layer are retained to update the latest PS network and the previous DS network, as well as the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the new network. This process is performed stepwise until the minimum allowable goodness-of-fit threshold is reached. Meanwhile, the network fused in the last step is retained as a highly coherent new network of PS and DS at different scales, and the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network are recorded.

[0013] Furthermore, in step S4, the specific steps are as follows: S41, set a constraint on the number of homogeneous points, gradually relax the goodness-of-fit threshold constraint from large to small, and select a new set of DS candidate points after excluding the latest PS point set; S42, combine the latest layer PS network with the previous layer DS network, execute steps S13-S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of the network nodes; S43, set a coherence threshold, and select candidate points from the DS candidate point set in this layer that have high coherence with the latest layer PS network and the previous layer DS network. The DS point set is used as the network of this layer. At the same time, the highly coherent edges connecting the PS and DS point sets of this layer are retained to update the latest PS network and the previous DS network, as well as the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the new network. S44, gradually relax the goodness-of-fit threshold constraint from large to small, update the DS candidate point set, and repeat steps S42 to S43 until the smallest allowable goodness-of-fit threshold is reached. At the same time, the network fused in the last step is retained as a highly coherent new network of PS and DS at different scales, and the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the network are recorded.

[0014] According to the above scheme, the specific process in step S5 is as follows: Based on the ensemble coherence weighted least squares method, the residual topography and deformation rate information of all multi-scale highly coherent scatterers are calculated. Further, the specific steps in step S5 are: S51, using the residual elevation gradients and deformation rate gradients of all connected PS and DS points obtained in step 4, arrange them in rows to form a differential elevation vector. Differential deformation rate vector The set coherence coefficients of highly coherent edges in the network are used as the weight matrix. The residual topography and deformation rate of all PS and DS multiscale highly coherent scatterers are estimated using the weighted least squares method of equation (3): (3), among which, This represents the mapping matrix of highly coherent edges connecting PS and DS points to the new multi-scale scatterer network. Finally, the residual terrain and deformation rate information of all multi-scale highly coherent scatterers can be calculated using equation (3).

[0015] The present invention also provides an InSAR precision measurement system for multi-scale highly coherent scatterers, which has a computer program that can execute the InSAR precision measurement method for multi-scale highly coherent scatterers.

[0016] The InSAR precision measurement method and system based on multi-scale highly coherent scatterers of the present invention has the following beneficial effects: 1. The present invention provides a multi-scale highly coherent scatterer joint InSAR precision measurement method, which offers a new framework for joint scatterer processing for accurate estimation of time-series InSAR parameters.

[0017] 2. The InSAR precision measurement method of multi-scale highly coherent scatterers of the present invention can be used as a further optimization step of the current PS and DS identification algorithm and as a subsequent time-series InSAR analysis process.

[0018] 3. The InSAR precision measurement method of the present invention, which combines multi-scale highly coherent scatterers, not only maintains the accuracy of measurement points but also greatly improves the density and coverage of precision measurement points.

[0019] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0020] Figure 1 This is the overall flowchart of Example 1 of the InSAR precision measurement method based on multi-scale highly coherent scatterers of the present invention; Figure 2 This is an average ensemble coherence coefficient diagram of the PS and DS points in Example 3 of the InSAR precision measurement method for multi-scale high-coherence scatterers of the present invention. Figure 3 This is a constrained high-coherence network diagram showing the connection between PS and DS points in Example 3 of the InSAR precision measurement method for multi-scale high-coherence scatterers of the present invention. Figure 4Figures (a) and (b) are respectively the elevation map of the three-dimensional reconstruction of the PS and DS multi-scale scatterers and the deformation map of the interferometric measurement in Example 3 of the InSAR precision measurement method of multi-scale highly coherent scatterers of the present invention. Figure 5 This is a schematic diagram of the architecture of the InSAR precision measurement system based on multi-scale highly coherent scatterers of the present invention. Detailed Implementation

[0021] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments, but these embodiments should not be construed as limiting the present invention.

[0022] This invention proposes a multi-scale, highly coherent scatterer-based InSAR precision measurement technique. This technique is unaffected by the initial identification accuracy of the input PS and DS targets. Unlike the ICOPS technique, this invention is a novel joint scatterer processing framework oriented towards accurate estimation of time-series InSAR parameters. It organically combines robust identification of PS and DS targets with parameter estimation, identifying PS and DS targets whose parameters can be accurately inverted. This processing framework can also serve as a further optimization step for current PS and DS identification algorithms and a time-series InSAR processing framework.

[0023] Example 1 See Figure 1 The specific steps of the InSAR precision measurement method based on multi-scale highly coherent scatterers of the present invention are as follows: S1. Input multi-temporal interferometric SAR image data and initially construct the initial layer PS node network; S2. Set constraints from tight to wide to update the PS candidate point set and construct a new PS network of the same scale with high coherence that merges the previous layer PS network; S3. Identify adjacent homogeneous points of image pixels that are not PS points, and perform phase optimization on image pixels that are not PS points; S4. Set constraints from tight to wide to update the DS candidate point set, and construct a new highly coherent network of PS and DS at different scales that integrates the latest layer PS network and the previous layer DS network. Joint inversion of elevation or deformation rate of multi-scale highly coherent scatterers using S5, PS, and DS methods.

[0024] Example 2 The steps in this embodiment are the same as in Embodiment 1, except that each step is applied to a specific instance. Specifically, it includes the following steps: S1. For the phase-preprocessed multi-temporal interferometric SAR image data, set an initial amplitude deviation threshold to a small value, and retain points below the threshold as the initial PS candidate point set; construct a Delaunay triangulation network with neighborhood distance constraints, calculate the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; set a coherence threshold, remove points with low average set coherence coefficients and edges with low set coherence coefficients, and retain high coherence points connected by high coherence edges to form a high coherence network of initial layer PS points, while retaining the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the high coherence edges connected by the network. The specific steps are as follows: S11. Input the phase-preprocessed multi-temporal SAR image data and calculate the amplitude deviation of the multi-temporal interferometric SAR image; S12. Set the initial amplitude deviation threshold to a small value, and retain the pixels in the image below the threshold as the initial PS candidate point set; S13. Set the distance constraint threshold for the edges connecting the initial PS candidate points and construct a Delaunay triangulation network with neighborhood distance constraints. S14. Utilize CPU parallelism and GPU computing strategies to rapidly search for parameters of the model maximizing the ensemble coherence coefficient. (1) Obtain the set coherence coefficient of any connected edge in the candidate PS point set. Residual elevation gradient Average deformation rate gradient And the average set coherence coefficient of the edges connected to any PS candidate point.

[0025] In the optimization formula (1) For any two adjacent PS candidate points and The phase difference operator, Indicates the first One PS candidate point in Differential interferometric phase in spectral interferometric SAR images and They represent the first The PS candidate point at the th The gradient form of the simulated residual terrain and linear deformation interferometric phase in the interferogram can be expressed as follows: (2) in, and They represent the first Vertical baseline and observation time of each image. and These represent the radar incident angle and wavelength, respectively.

[0026] S15. Set a coherence threshold, remove points with low average set coherence coefficients and edges with low set coherence coefficients, retain high coherence points connected by high coherence edges, and form a high coherence network of PS nodes in the initial layer. At the same time, retain the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the high coherence edges connected to the network.

[0027] S2. Gradually relax the amplitude deviation threshold constraint from small to large, and select a new PS candidate point set; combine with the previous layer PS node network, and progressively execute the triangulation construction in step S1, as well as the calculation steps of the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; set a coherence threshold, and select candidate points from the PS candidate point set of this layer network that are highly coherently fused with the previous layer PS network as the PS point set of this layer network. At the same time, retain the highly coherent edges connected between the PS point sets of this layer network to update the previous layer PS high coherence network, as well as the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network; progressively execute this process until the most relaxed amplitude deviation threshold is reached, and retain the network fused in the last step as a new highly coherent PS network of the same scale, and record the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network. The specific steps are as follows: S21. Reset the PS candidate point screening threshold to a smaller value and select a new PS candidate point set; S22. Combine the high coherence network of the PS points in the previous layer, execute steps S13 to S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between the network edges, and the average set coherence coefficient of the network nodes. S23. Set a coherence threshold, select candidate points that are highly coherently fused with the previous layer PS network from the PS candidate point set of this layer network, and use them as the PS point set of this layer network. At the same time, retain the highly coherent edges that connect the PS point sets of this layer network to update the previous layer PS highly coherent network, as well as the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges that connect the network. S24. Gradually relax the amplitude deviation threshold constraint from small to large, update the PS candidate point set, and repeat steps S22 to S23 until the most relaxed amplitude deviation threshold is reached. At the same time, retain the network fused in the last step as a new PS network with high coherence at the same scale, and record the set coherence coefficient, residual elevation gradient and deformation rate gradient of the high coherence edges connected to the network.

[0028] S3. Based on the temporal similarity metric, identify the adjacent homogeneous points of image pixels of non-PS points and record the number of adjacent homogeneous points of non-PS points; calculate the coherence matrix of non-PS points based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points; perform singular value decomposition on the coherence matrix of non-PS points, and use the maximum scattering mechanism separation principle to separate the phase of the strongest scattering mechanism within the image pixels of non-PS points, using it as the optimized phase, and simultaneously record the phase fitting goodness of non-PS points. The specific steps are as follows: S31. Based on the temporal similarity criterion, identify the adjacent homogeneous points of image pixels that are not PS points, and record the number of adjacent homogeneous points of non-PS points. S32. Calculate the coherence matrix of non-PS points based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points; S33. Perform singular value decomposition on the coherence matrix of non-PS points, and use the maximum scattering mechanism separation principle to separate the phase of the strongest scattering mechanism in the image pixels of non-PS points, and use it as the optimized phase. At the same time, record the phase fitting goodness of non-PS points.

[0029] S4. Set a constraint on the number of homogeneous points, and gradually relax the goodness-of-fit threshold constraint from large to small, selecting a new set of DS candidate points after excluding the latest PS point set; combine the latest PS network and the previous DS network, and progressively execute the triangulation construction in step S1, the calculation steps of the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes; set a coherence threshold, and select candidate points from the DS candidate point set of this layer that are highly coherently integrated with the latest PS network and the previous DS network. The DS point set is used as the network layer, and the highly coherent edges connecting the PS and DS point sets of this layer are preserved to update the latest PS network and the previous DS network, as well as the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the new network. This process is performed stepwise until the minimum allowable goodness-of-fit threshold is reached. The network fused in the last step is retained as a new highly coherent network of PS and DS at different scales, and the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network are recorded. The specific steps are as follows: S41. Set a constraint on the number of homogeneous points, gradually relax the goodness-of-fit threshold constraint from large to small, and select a new DS candidate point set after excluding the latest PS point set. S42. Combine the latest layer PS network with the previous layer DS network, execute steps S13 to S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of the network nodes. S43. Set a coherence threshold, select candidate points from the DS candidate point set in the network layer that are highly coherently fused with the latest PS network and the previous DS network, and use them as the DS point set of the network layer. At the same time, retain the highly coherent edges connecting the PS and DS point sets of the network layer to update the coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connecting the latest PS network, the previous DS network, and the new network. S44. Gradually relax the goodness-of-fit threshold constraint from large to small, update the DS candidate point set, and repeat steps S42 to S43 until the smallest allowable goodness-of-fit threshold is reached. At the same time, retain the network fused in the last step as a new highly coherent network of PS and DS at different scales, and record the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the network.

[0030] S5. Based on the ensemble coherence weighted least squares method, the residual topography and deformation rate information of all multi-scale highly coherent scatterers are calculated. The specific steps are as follows: S51. Using the residual elevation gradients and deformation rate gradients of all connected PS and DS points obtained in step 4, arrange them in rows to form a differential elevation vector. Differential deformation rate vector The set coherence coefficients of highly coherent edges in the network are used as the weight matrix. The residual topography and deformation rate of all PS and DS multiscale highly coherent scatterers are estimated using the weighted least squares method of equation (3): (3) in, This represents the mapping matrix of highly coherent edges connecting PS and DS points to the new multi-scale scatterer network. Finally, the residual terrain and deformation rate information of all multi-scale highly coherent scatterers can be calculated using equation (3).

[0031] It should be understood that the sequence number of each step in the above embodiments does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0032] Example 3 This embodiment is basically the same as embodiment 1, except that this embodiment uses the steps of embodiment 1 to perform three-dimensional reconstruction of urban buildings and deformation rate measurement on 37 multi-temporal SAR images to obtain the measurement effect image. Figure 2 This is the average set coherence coefficient of the PS and DS points calculated in this embodiment of the invention; Figure 3 This is a constrained highly coherent network connecting the PS and DS points in an embodiment of the present invention; Figure 4In the figures (a) and (b), respectively, the elevation map of the three-dimensional reconstruction of the PS and DS multi-scale scatterers and the deformation map of the interferometric measurement are shown. It can be seen that after the construction of the multi-scale scatterer constraint network of the present invention, the elevation and deformation rate of the complex urban scene are accurately measured, which shows the feasibility of the technical solution of the present invention.

[0033] Example 4 See Figure 5 The InSAR precision measurement system based on multi-scale highly coherent scatterers of the present invention includes the following: Initial layer construction module: Input multi-temporal interferometric SAR image data, and initially construct the initial layer PS node network; Candidate point set update module: Sets constraints from tight to wide to update the PS candidate point set, and constructs a new PS network of the same scale with high coherence that merges the previous layer PS network; Phase optimization module: Identifies adjacent homogeneous points of image pixels that are not PS points, and performs phase optimization on the image pixels that are not PS points; Constraint Update Module: Sets constraints from tight to wide to update the DS candidate point set, and constructs a new highly coherent PS and DS network with different scales that integrates the latest PS network and the previous DS network; Joint Inversion Module: Joint inversion of elevation or deformation rate of multi-scale highly coherent scatterers using PS and DS.

[0034] Example 5 Those skilled in the art will understand that embodiments of this application can be provided as methods or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects.

[0035] Furthermore, this application may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0036] This embodiment is described with reference to the flowchart of the method and computer program product according to Embodiment 1 of this application and the method block diagram of Embodiment 3. It should be understood that each step or block in the flowchart or block diagram, as well as combinations of steps or blocks in the flowchart or block diagram, can be implemented by computer program instructions.

[0037] These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing device to produce a machine, such that the instructions, which are executed by the processor of the computer or other programmable data processing device, produce a multi-scale, highly coherent scatterer joint InSAR precision measurement system for implementing the specified functions.

[0038] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means that perform the specified function.

[0039] These computer program instructions may also be loaded onto a computer or other programmable data processing device to cause a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process, thereby providing steps for implementing a multi-scale, highly coherent scatterer-integrated InSAR precision measurement method.

[0040] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

[0041] The contents not described in detail in this specification are existing technologies known to those skilled in the art.

Claims

1. A precise InSAR measurement method using multi-scale highly coherent scatterers, characterized in that: Includes the following steps: S1. Input multi-temporal interferometric SAR image data and initially construct the initial layer PS node network; S2. Set constraints from tight to wide to update the PS candidate point set and construct a new PS network of the same scale with high coherence that merges the previous layer PS network; S3. Identify adjacent homogeneous points of image pixels that are not PS points, and perform phase optimization on image pixels that are not PS points; S4. Set constraints from tight to wide to update the DS candidate point set, and construct a new highly coherent network of PS and DS at different scales that integrates the latest layer PS network and the previous layer DS network. Joint inversion of elevation or deformation rate of multi-scale highly coherent scatterers using S5, PS, and DS methods.

2. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 1, characterized in that: The specific process of step S1 is as follows: For multi-temporal interferometric SAR image data with phase preprocessing, an initial amplitude deviation threshold is set to a small value, and points below the threshold are retained as the initial PS candidate point set. A Delaunay triangulation network with neighborhood distance constraints is constructed, and the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of the network nodes are calculated. A coherence threshold is set, and points with low average set coherence coefficients and edges with low set coherence coefficients are removed, while high coherence points connected by high coherence edges are retained to form a high coherence network of initial layer PS points. At the same time, the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the high coherence edges connected by the network are retained.

3. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 2, characterized in that: The specific steps of step S1 are as follows: S11. Input the phase-preprocessed multi-temporal SAR image data and calculate the amplitude deviation of the multi-temporal interferometric SAR image; S12. Set the initial amplitude deviation threshold to a small value, and retain the pixels in the image below the threshold as the initial PS candidate point set; S13. Set the distance constraint threshold for the edges connecting the initial PS candidate points and construct a Delaunay triangulation network with neighborhood distance constraints. S14. Utilize CPU parallelism and GPU computing strategies to rapidly search for parameters of the model maximizing the ensemble coherence coefficient. （1） Obtain the set coherence coefficient of any connected edge in the candidate PS point set. Residual elevation gradient Average deformation rate gradient And the average set coherence coefficient of the edges connected to any PS candidate point. In the optimization formula (1), For any two adjacent PS candidate points and The phase difference operator, Indicates the first One PS candidate point in Differential interferometric phase in spectral interferometric SAR images and They represent the first The PS candidate point at the th The gradient form of the simulated residual terrain and linear deformation interferometric phase in the amplitude interferogram can be expressed as follows: （2） in, and They represent the first Vertical baseline and observation time of each image. and These represent the radar incident angle and wavelength, respectively. S15. Set a coherence threshold, remove points with low average set coherence coefficients and edges with low set coherence coefficients, retain high coherence points connected by high coherence edges, and form a high coherence network of PS nodes in the initial layer. At the same time, retain the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the high coherence edges connected to the network.

4. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 3, characterized in that: The specific process of step S2 is as follows: Gradually relax the amplitude deviation threshold constraint from small to large to select a new set of PS candidate points; combine the PS node network of the previous layer, and gradually execute the Delaunay triangulation construction, as well as the calculation steps of the set coherence coefficient of network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of network nodes. Set a coherence threshold, select candidate points from the PS candidate point set of this layer network that are highly coherently fused with the PS network of the previous layer, and use them as the PS point set of this layer network. At the same time, retain the highly coherent edges connected between the PS point sets of this layer network to update the highly coherent PS network of the previous layer, as well as the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network. This process is executed step by step until the most lenient magnitude deviation threshold is reached. At the same time, retain the network fused in the last step as a new highly coherent PS network of the same scale, and record the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network.

5. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 4, characterized in that: The specific steps of step S2 are as follows: S21. Reset the PS candidate point screening threshold to a smaller value and select a new PS candidate point set; S22. Combine the high coherence network of the PS points in the previous layer, execute steps S13 to S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between the network edges, and the average set coherence coefficient of the network nodes. S23. Set a coherence threshold, select candidate points that are highly coherently fused with the previous layer PS network from the PS candidate point set of this layer network, and use them as the PS point set of this layer network. At the same time, retain the highly coherent edges that connect the PS point sets of this layer network to update the previous layer PS highly coherent network, as well as the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges that connect the network. S24. Gradually relax the amplitude deviation threshold constraint from small to large, update the PS candidate point set, and repeat steps S22 to S23 until the most relaxed amplitude deviation threshold is reached. At the same time, retain the network fused in the last step as a new PS network with high coherence at the same scale, and record the set coherence coefficient, residual elevation gradient and deformation rate gradient of the high coherence edges connected to the network.

6. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 5, characterized in that: The specific process of step S3 is as follows: Based on the temporal similarity metric, neighboring homogeneous points of image pixels of non-PS points are identified, and the number of neighboring homogeneous points of non-PS points is recorded. Based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points, the coherence matrix of non-PS points is calculated. Singular value decomposition is performed on the coherence matrix of non-PS points, and the phase of the strongest scattering mechanism within the image pixels of non-PS points is separated using the maximum scattering mechanism separation principle. This phase is used as the optimized phase, and the phase fitting goodness of non-PS points is recorded. The specific steps of step S3 are as follows: S31. Based on the temporal similarity criterion, identify the adjacent homogeneous points of image pixels that are not PS points, and record the number of adjacent homogeneous points of non-PS points. S32. Calculate the coherence matrix of non-PS points based on the pixel values ​​of the multi-temporal SAR images of the identified homogeneous points; S33. Perform singular value decomposition on the coherence matrix of non-PS points, and use the maximum scattering mechanism separation principle to separate the phase of the strongest scattering mechanism in the image pixels of non-PS points, and use it as the optimized phase. At the same time, record the phase fitting goodness of non-PS points.

7. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 6, characterized in that: The specific process of step S4 is as follows: Set a constraint on the number of homogeneous points, gradually relax the goodness-of-fit threshold constraint from large to small, and select a new DS candidate point set after excluding the latest PS point set; combine the latest layer PS network with the previous layer DS network, and gradually execute the steps of triangulation construction in step S1, calculation of the set coherence coefficient of network edges, residual elevation gradient and deformation rate gradient between network edges, and average set coherence coefficient of network nodes. Set a coherence threshold, and select candidate points from the DS candidate point set in the network layer that are highly coherently fused with the latest PS network and the previous DS network. These points are then used as the DS point set of the network layer. At the same time, retain the highly coherent edges connecting the PS and DS point sets of the network layer to update the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the new network. This process is executed step by step until the minimum allowable goodness-of-fit threshold is reached. Meanwhile, the network fused in the last step is retained as a new highly coherent network of PS and DS at different scales, and the set coherence coefficient, residual elevation gradient, and deformation rate gradient of the highly coherent edges connected to the network are recorded.

8. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 7, characterized in that: The specific steps of step S4 are as follows: S41. Set a constraint on the number of homogeneous points, gradually relax the goodness-of-fit threshold constraint from large to small, and select a new DS candidate point set after excluding the latest PS point set. S42. Combine the latest layer PS network with the previous layer DS network, execute steps S13 to S14, construct a new constrained triangular network, calculate and record the set coherence coefficient of the network edges, the residual elevation gradient and deformation rate gradient between network edges, and the average set coherence coefficient of the network nodes. S43. Set a coherence threshold, select candidate points from the DS candidate point set in the network layer that are highly coherently fused with the latest PS network and the previous DS network, and use them as the DS point set of the network layer. At the same time, retain the highly coherent edges connecting the PS and DS point sets of the network layer to update the coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connecting the latest PS network, the previous DS network, and the new network. S44. Gradually relax the goodness-of-fit threshold constraint from large to small, update the DS candidate point set, and repeat steps S42 to S43 until the smallest allowable goodness-of-fit threshold is reached. At the same time, retain the network fused in the last step as a new highly coherent network of PS and DS at different scales, and record the set coherence coefficient, residual elevation gradient and deformation rate gradient of the highly coherent edges connected to the network.

9. The InSAR precision measurement method based on multi-scale highly coherent scatterers according to claim 8, characterized in that: The specific process of step S5 is as follows: Based on the least squares method with coherence weighting, the residual topography and deformation rate information of all highly coherent scatterers are calculated. The S5 step is further detailed as follows: Using the residual elevation gradients and deformation rate gradients of all connected PS and DS points obtained, a differential elevation vector is formed by arranging them in rows. Differential deformation rate vector The set coherence coefficients of highly coherent edges in the network are used as the weight matrix. The residual topography and deformation rate of all PS and DS multiscale highly coherent scatterers are estimated using the weighted least squares method of formula (3): （3） in, The mapping matrix represents the highly coherent edges connecting PS and DS points to the new multi-scale scatterer network. The residual terrain and deformation rate information of all multi-scale highly coherent scatterers can be calculated by formula (3).

10. A multi-scale, highly coherent scatterer combined InSAR precision measurement system, comprising a computer program, characterized in that: The computer program is capable of performing the InSAR precision measurement method using multi-scale highly coherent scatterers as described in any one of claims 1 to 9.

Images