Multi-channel InSAR phase unwrapping method, system, terminal and storage medium

Abstract

The present invention discloses a multi-channel InSAR phase unwrapping method, system, terminal, and storage medium. The method comprises: constructing a mathematical model of a congruence equation system based on the basic principles of InSAR; selecting multiple integer congruence equations corresponding to K different moduli to equally divide the multiple integer congruence equation systems; based on the acquired multi-channel InSAR data, calculating the deviation value sets that may appear in the search process of the multiple integer congruence equation systems and finding the intersection, determining whether the intersection is {0, 1, 2}; if so, directly solving the terrain data estimation value set, and obtaining the final terrain data by finding the intersection and averaging; otherwise, selecting a new congruence equation with a different modulus to replace any congruence equation in the multiple integer congruence equation systems, resolving the terrain data estimation value set, and obtaining the final terrain data by finding the intersection and averaging. The present invention significantly improves the unwrapping success rate and the accuracy of the reconstructed target value.

Description

Technical Field

[0001] The present invention relates to the field of synthetic aperture radar interferometry technology, and in particular to a multi-channel InSAR phase unwrapping method, system, terminal and storage medium. Background Art

[0002] Phase unwrapping is a key processing step in Interferometric Synthetic Aperture Radar (InSAR) and plays a vital role in InSAR data processing. Traditional phase unwrapping algorithms are mainly categorized into: single-baseline phase unwrapping, which requires the phase continuity assumption; multi-baseline phase unwrapping, which requires multiple interferograms to be unwrapped based on the Chinese Residue Theorem (CRT); however, small errors in the remainder can lead to significant deviations in the reconstructed target number; and deep learning-based phase unwrapping, which requires extensive training data to obtain phase noise characteristics to guide phase unwrapping.

[0003] Currently, the representative robust Chinese Remainder Theorem phase unwrapping algorithms mainly include: the search-based robust CRT phase unwrapping algorithm, which searches for the minimum value of all remainder deviation values ​​as a reliable congruence equation group for unwrapping, thereby improving the robustness of the algorithm. The disadvantage of this method is that each pixel point on the interference pattern needs to undergo a repeated search process, which makes its computational efficiency low; the closed robust CRT-based phase unwrapping algorithm, which performs phase unwrapping by performing a congruence equation group consisting of L-1 new congruence equations by taking the difference between L-1 congruence equations and the first reference congruence equation. The disadvantage of this method is that the selection of its reference congruence equation is not necessarily optimal, and improper selection may lead to unwrapping failure; the multi-stage robust CRT uses the modular grouping method based on the closed robust CRT algorithm to improve the value of its greatest common divisor M, thereby improving the remainder error bound, but does not fundamentally improve the range of the remainder error bound, which is still less than one-quarter of the greatest common divisor. The existing technology still has problems such as weak robustness and small residual error bound, and the unwrapping effect is still not ideal in some noisy areas. Summary of the Invention

[0004] The present invention provides a multi-channel InSAR phase unwrapping method, system, terminal and storage medium to solve the problem of weak robustness of phase unwrapping in the prior art.

[0005] In a first aspect, a multi-channel InSAR phase unwrapping method is provided, comprising the following steps:

[0006] S1: Construct a mathematical model of congruence equations based on the basic principles of interferometric synthetic aperture radar;

[0007] S2: Select multiple integer congruence equations from the existing K integer congruence equations corresponding to different moduli and divide them into multiple groups to obtain multiple integer congruence equations, and each integer congruence equation group satisfies that the moduli are mutually prime;

[0008] S3: Based on the acquired multi-channel InSAR data, calculate the deviation value set that may appear in each set of integer congruence equations during the search process, and find the intersection of multiple deviation value sets, and determine whether the intersection is {0, 1, 2}. If so, proceed to step S4; otherwise, proceed to step S5;

[0009] S4: Calculate a plurality of estimated value sets of terrain data based on the current plurality of integer congruence equations, obtain the final terrain data by intersecting and averaging the values, and then end;

[0010] S5: Select a new congruence equation with a different modulus to replace any congruence equation in the multiple sets of integer congruence equations, and then recalculate and determine whether the intersection of the deviation value set is {0,1,2}; if so, go to step S4; if not, integrate the multiple sets of integer congruence equations in step S2 and the current multiple sets of new integer congruence equations as the current multiple sets of integer congruence equations, and go to step S4.

[0011] Furthermore, in step S1, the mathematical model of the congruence equations is obtained by the following method:

[0012] The terrain data is represented by N, and the following real congruence equation model can be established; the terrain data is the terrain height or terrain shape variable;

[0013] N=n i *M i +r i ,i=1,2,....,L

[0014] Where n i is a fuzzy number, M i is the fuzzy height or fuzzy deformation, r i The interference phase is expressed as The height value or deformation value corresponding to the time; represents the winding phase; L represents the number of modules, i represents the i-th module or the i-th congruence equation;

[0015] If M i There is a common greatest common factor M, and the above formula can be further simplified as:

[0016]

[0017] make Subtract r from both sides of the above equation cAnd divided by M, we can get the following integer congruence equation model:

[0018] N0≡q i modΓ i

[0019] Where, Γ i represents the modulus of the i-th congruence equation;

[0020] The reconstruction formula of terrain data N is: N = M*N0+r c .

[0021] Furthermore, in step S2, each set of integer congruence equations is expressed as follows:

[0022]

[0023] Where, represents the estimated value of N0 in the Ath group of integer congruence equations, m AL , Γ AL ,q AL , Δq AL ,Δp AL represents the folding integer, modulus, remainder, remainder error, and complementary error of the Lth congruence equation in the Ath group of integer congruence equations, and Δq AL ∈{-1,0,1},Δp AL ∈{-1,0,1}.

[0024] Furthermore, in step S3, the estimated value of N0 is obtained by using the classical CRT for the Ath group of integer congruence equations:

[0025]

[0026] Where, is the deviation value of N0; Γ A =Γ A1 *Γ A2 *....*Γ AL ; γ Ai =Γ A / Γ Ai ;

[0027] When the intersection of the deviation value sets formed by the possible values ​​of ΔN0 of multiple sets of integer congruence equations during the search process is {0,1,2}, then enter step S4; if the intersection of the deviation value sets is {0,1,2,ΔN1}, ΔN1>>2, then enter step S5.

[0028] Furthermore, step S4 includes:

[0029] Based on the current multiple sets of integer congruence equations, multiple estimated value sets of N0 are calculated;

[0030] Obtaining a plurality of estimated value sets of terrain data according to a plurality of estimated value sets of N0 and a reconstruction formula of terrain data N;

[0031] The final terrain data is obtained by finding the intersection of the estimated value sets of multiple terrain data and taking the average value.

[0032] Furthermore, in step S5, the modulus of the selected new congruence equation satisfies the following conditions:

[0033] ΔN1=k 2L+1 *Γ 2L+1 +Δg 2L+1 ,Δg 2L+1 ∈{3,4,....,Γ 2L+1 -3}

[0034] Where k 2L+1 is the folding integer, Δg 2L+1 is the remainder.

[0035] Furthermore, in step S5, if the intersection of the recalculated deviation value sets is {0, 1, 2}, multiple estimated value sets of N0 are calculated based on the current multiple new integer congruence equations; multiple estimated value sets of terrain data are obtained based on the multiple estimated value sets of N0 and the reconstruction formula of the terrain data N; the final terrain data is obtained by calculating the intersection of the multiple estimated value sets of terrain data and taking the average value;

[0036] If the intersection of the recalculated deviation value sets is {0, 1, 2, ΔN2}, ΔN2>>2, multiple estimated value sets of N0 are calculated based on the multiple sets of integer congruence equations in step S2 and the current multiple new sets of integer congruence equations; multiple estimated value sets of terrain data are obtained based on the multiple estimated value sets of N0 and the reconstruction formula of the terrain data N; the final terrain data is obtained by calculating the intersection of the multiple estimated value sets of terrain data and taking the average value.

[0037] In a second aspect, a multi-channel InSAR phase unwrapping system is provided, comprising:

[0038] A model building module is used to build a mathematical model of a group of congruence equations based on the basic principles of interferometric synthetic aperture radar;

[0039] A congruence equation group construction module is used to select multiple integer congruence equation groups from the existing K integer congruence equation groups corresponding to different moduli and divide them into multiple groups to obtain multiple groups of integer congruence equation groups, where each group of integer congruence equation groups satisfies that the moduli are mutually prime.

[0040] The terrain data reconstruction module is used to calculate the deviation value set that may appear in the search process for each set of integer congruence equations based on the acquired multi-channel InSAR data, and to find the intersection of multiple deviation value sets and determine whether the intersection is {0, 1, 2};

[0041] If so, a plurality of estimated value sets of terrain data are calculated based on the current multiple sets of integer congruence equations, and the final terrain data is obtained by finding the intersection and taking the average value;

[0042] If not, select a new congruence equation with a different modulus to replace any one of the multiple new sets of integer congruence equations, and then recalculate and determine whether the intersection of the deviation value sets is {0,1,2}; if so, obtain a set of estimated values ​​of multiple terrain data based on the current multiple sets of new integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average; if not, obtain a set of estimated values ​​of multiple terrain data based on the multiple sets of integer congruence equations constructed by the congruence equation group construction module and the current multiple sets of new integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average.

[0043] In a third aspect, an electronic terminal is provided, comprising:

[0044] a memory having a computer program stored thereon;

[0045] The processor is configured to load and execute the computer program to implement the multi-channel InSAR phase unwrapping method as described above.

[0046] In a fourth aspect, a computer-readable storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, the multi-channel InSAR phase unwrapping method as described above is implemented.

[0047] The present invention proposes a multi-channel InSAR phase unwrapping method, system, terminal, and storage medium. By dividing a set of moduli into two or several groups of moduli combinations, the method searches for all possible estimated values ​​when the integer part of the remainder occurs in [-1, 0, 1], and then performs an intersection operation on them to reconstruct the accurate target value. This method improves the robustness of the Chinese remainder theorem by increasing the redundancy of the moduli. Compared with the existing methods, the range of the remainder error bound is less than or equal to M / 4 or M / 2. This method further improves the range of the remainder error bound and expands it to less than or equal to M, where M is the greatest common divisor under the same group of moduli combinations. The unwrapping success rate of the present invention is significantly improved compared to other methods. The root mean square error between the reconstructed target value and the true target value is also much smaller than that of existing algorithms, resulting in higher reconstruction accuracy. It is worth noting that the size of the deviation value is fixed after the modulus is given, so this method does not need to repeatedly calculate the size of the deviation value. Therefore, its computational complexity is similar to that of the classic CRT algorithm, and it has high computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0049] Figure 1 This is a flow chart of a multi-channel InSAR phase unwrapping method provided by an embodiment of the present invention;

[0050] Figure 2 This is a comparative analysis chart of TFR performance of four different robust CRT algorithms provided by an embodiment of the present invention;

[0051] Figure 3 This is the terrain height inversion of different robust CRT algorithms provided by the embodiments of the present invention under the residual error bound, wherein (a), (d), and (g) are 256×256 interference patterns simulated by the peak function, (b), (e), and (h) are the terrain height inversion results for the residual noise size within the residual error bound under the method of the present invention, closed robust CRT, and multi-stage robust CRT, respectively, and (c), (f), and (i) are the difference images between the estimated terrain height and the actual value inverted by the method of the present invention, closed robust CRT, and multi-stage robust CRT, respectively. DETAILED DESCRIPTION

[0052] To make the objectives, technical solutions, and advantages of the present invention more apparent, the technical solutions of the present invention will be described in detail below. Obviously, the embodiments described are only some of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other implementations obtained by those of ordinary skill in the art without inventive effort are within the scope of protection of the present invention.

[0053] The present invention mainly considers that with the development of synthetic aperture radar interferometry (InSAR) technology, the application of InSAR technology is becoming more and more extensive, such as long-term monitoring of surface elevation and deformation monitoring, landslide monitoring, urban and mining area ground subsidence and stability of artificial buildings. Phase unwrapping is one of the key steps in InSAR technology processing, and how to better phase unwrapping is a problem worth exploring. The present invention establishes a real congruence equation group model for the surface elevation h(s) or the surface deformation d(s) based on the multi-channel InSAR system parameters and the interferometric phase information obtained by measurement, and converts the model into an integer congruence equation group model by dividing by the greatest common factor M and rounding, thereby converting the problem of solving the real numbers h(s) and d(s) into the problem of solving the unknown integer N0 using the Chinese remainder theorem. In the process of solving N0, the present invention divides the congruence equations into two or more groups to obtain a set of estimated values ​​for the target integer N0, and takes the average of the intersection of the two or more estimated value sets as the final N0 estimate. In this case, the error of the estimated value is minimized, thereby ensuring that the reconstructed surface elevation h(s) or surface deformation d(s) is more accurate. The grouping of the present invention is based on the requirements of the classical Chinese remainder theorem for the modulus, and only requires that the modulus of each group meet the coprime requirement. The technical solution of the present invention is described in detail below with reference to specific embodiments.

[0054] The embodiment of the present invention provides a multi-channel InSAR phase unwrapping method, such as Figure 1 As shown, the following steps are included:

[0055] S1: Construct a mathematical model of congruence equations based on the basic principles of interferometric synthetic aperture radar.

[0056] Consider a system with a working wavelength of λ and an effective baseline length of B. i And a multi-channel InSAR system with a viewing angle of θ. The terrain height h(s) of the (i, j)th pixel s can be expressed as:

[0057] h(s)=((λ*r(s)*sinθ) / (4π*B i ))*φ i (s) (1)

[0058] Where r(s) is the slant distance between the target s and the main antenna corresponding to the (i, j) pixel, φi (s) is the absolute interferometric phase of the target in the i-th interferogram.

[0059] Similarly, after removing the terrain phase, the surface deformation variable d(s) of the (i, j)th pixel can be expressed as:

[0060] d(s)=(λ / 4π)*φ i (s) (2)

[0061] d(s) is defined as the magnitude of surface deformation under different time baselines.

[0062] The absolute phase and winding phase have the following relationship:

[0063]

[0064] Where n i (s) is an unknown integer, called a fuzzy number, is the winding phase.

[0065] Assuming that there is no deformation when reconstructing h(s) (or the terrain phase has been removed when reconstructing the deformation d(s)), we can know the height corresponding to the phase change of 2π in the i-th interference pattern - the fuzzy height (or fuzzy deformation).

[0066] Mi(s)=(λ*r(s)*sinθ) / (2*B i ) (or M i (s)=λ / 2) (4)

[0067] And the interference phase is The corresponding height value (or deformation value) is:

[0068]

[0069] The terrain height h(s) and terrain deformation size d(s) in (1) and (2) are uniformly expressed by N(s), and the following real congruence equation model can be established:

[0070] N(s)=n i (s)*M i (s)+r i (s),i=1,2,....,L (6)

[0071] Simplify the above formula:

[0072] N=n i *M i +r i ,i=1,2,....,L (7)

[0073] Where L represents the number of modules, i represents the i-th module or the i-th congruence equation;

[0074] If M i There is a common greatest common factor M, and the above formula can be further simplified as:

[0075]

[0076] make Subtract rc from both sides of equation (8) and divide by M to obtain the following integer congruence equation model:

[0077] N0≡q i modΓ i (9)

[0078] Where, Γ i represents the modulus of the i-th congruence equation.

[0079] So N can be reconstructed as follows:

[0080] N=M*N0+r c (10)

[0081] Where N0 can be determined by classic CRT as:

[0082]

[0083] Where, Γ=Γ1*Γ2*....*Γ L , satisfy γ i Defined as γ i =Γ / Γ i .

[0084] S2: Select multiple integer congruence equations from the existing K integer congruence equations corresponding to different moduli and divide them into multiple groups to obtain multiple integer congruence equations, and each integer congruence equation group satisfies that the moduli are mutually prime.

[0085] In the following embodiments, the system is described by taking the equations divided equally into two groups of integer congruence equations as an example.

[0086] Assume that the modules corresponding to the two integer congruence equations satisfying pairwise coprime are {Γ a1 ,Γ a2 ....Γ aL} and {Γ b1 ,Γ b2 ....Γ bL}, K>2L. From equation (9), we can see that the corresponding two integer congruence equations are:

[0087]

[0088] Among them, m ai represents the folded integer or fuzzy number in the i-th congruence equation in the a-th group of integer congruence equations.

[0089] When the remainder q ai and q bi Δq appears ai and Δq bi When the error is , the above integer congruence equations become:

[0090]

[0091] It can be seen that the estimated values ​​obtained by CRT inversion are:

[0092]

[0093] When Δq ai and Δq bi When the size of does not exceed 1, ΔN0(a) and ΔN0(b) each have 3 L Possible errors, and There are also 3 L possible results. Obviously, when the remainder q ai and q bi When there is no error, the value obtained by CRT inversion must satisfy Therefore, if Δq ai and Δq bi On the basis of each, a complementary error Δp ​​is added. ai ,Δp bi ∈{-1,0,1}, then when Δp ai =-Δq ai ,Δp bi =-Δq bi The result obtained by CRT inversion is the error-free N0. ai ,Δq bi ∈{-1,0,1}, and Δp ai ,Δp bi ∈{-1,0,1}, so when Δq ai and Δq bi When it is determined, it can be determined from Δp ai and Δp bi 3 of each L Find a set among the error sets that satisfies the condition Δp ai =-Δq ai and Δp bi =-Δq bi Therefore, the 3 obtained by adding the error to each group LThe estimated values ​​must contain the error-free N0. After adding the complementary error, the corresponding integer congruence equations become:

[0094]

[0095] The estimated value obtained using CRT at this time is:

[0096]

[0097] Without loss of generality, the group subscripts are omitted here.

[0098] S3: Based on the acquired multi-channel InSAR data, calculate the deviation value set that may appear in each set of integer congruence equations during the search process, and find the intersection of multiple deviation value sets, and determine whether the intersection is {0, 1, 2}. If so, proceed to step S4, otherwise proceed to step S5.

[0099] According to the following formula in the "Lectures on Elementary Number Theory":

[0100] (A+B)mod C=(A mod C+B mod C)mod C (18)

[0101] So formula (17) can be converted into the following formula:

[0102]

[0103] Right now:

[0104]

[0105] Among them, the deviation value of N0 And (Δp i +Δq i ) is known to have a range of {-2, -1, 0, 1, 2}. When a set of moduli is given, its γ i 、 Γ is also known; that is, when the modulus is given, the range of the difference between the target value and the estimated value obtained by inversion is known. Because the estimated values ​​under different modulus combinations all contain the target value, when two or more sets of congruence equations are used to solve the target value estimate set and find the intersection, the target value must be included in the intersection, in which case ΔN0 = 0. However, if ΔN0>>2 exists in the intersection, it will lead to the incorrect target value selection. Therefore, how to eliminate the situation where ΔN0>>2 exists in the intersection is the next step.

[0106] When 2L modules are selected from K modules and divided into two groups according to the condition that each module combination satisfies the coprime condition, its γ can be calculated based on the module information. ai ,γ bi , Γ a ,Γ b According to (19), the estimated deviation value ΔN under different modulus combinations is a0 ,ΔN b0 The value of is also determined. That is:

[0107]

[0108] If {ΔN0(a)}∩{ΔN0(b)}={0,1,2}, it belongs to the first case. This combination is the simplest and most expected modular combination. At this time, enter step S4. If {ΔN0(a)}∩{ΔN0(b)}={0,1,2,ΔN1}, ΔN1>>2, then if two sets of modular combinations are used for inversion to reconstruct the target value, the following possible values ​​may be obtained. At this time, it is impossible to accurately determine which value in the estimated value set is the correct target value, and the process goes to step S5.

[0109] S4: Based on the current multiple sets of integer congruence equations, a set of estimated values ​​of multiple terrain data is calculated, and the final terrain data is obtained by the method of finding the intersection and taking the average value, and then the process ends.

[0110] Specifically, step S4 includes:

[0111] Based on the current multiple sets of integer congruence equations, multiple estimated value sets of N0 are calculated;

[0112] Obtaining a plurality of estimated value sets of terrain data according to a plurality of estimated value sets of N0 and a reconstruction formula of terrain data N (terrain height or terrain shape);

[0113] The final terrain data is obtained by finding the intersection of the estimated value sets of multiple terrain data and taking the average value.

[0114] Of course, in other embodiments, it is also possible to obtain N0 by intersecting and averaging multiple estimated value sets of N0, and then calculate the final terrain data based on the obtained N0.

[0115] S5: Select a new congruence equation with a different modulus to replace any one of the multiple new integer congruence equations, then recalculate and determine whether the intersection of the deviation value sets is {0, 1, 2}; if so, proceed to step S4; if not, integrate the multiple integer congruence equations of step S2 and the current multiple new integer congruence equations as the current multiple new integer congruence equations (the integration here can be understood as combining the X groups of integer congruence equations of step S2 and the new integer congruence equations after replacing one congruence equation to form X+1 new integer congruence equations), and proceed to step S4.

[0116] Specifically, when {ΔN0(a)}∩{ΔN0(b)}={0,1,2,ΔN1}, ΔN1>>2 occurs, it is impossible to accurately determine which value in the estimated value set is the correct target value. Therefore, a modulus that satisfies the following conditions is selected from the remaining K-2L moduli:

[0117] ΔN1=k 2L+1 *Γ 2L+1 +Δg 2L+1 ,Δg 2L+1 ∈{3,4,....,Γ 2L+1 -3} (22)

[0118] Where k 2L+1 is the folding integer, Δg 2L+1 is the remainder.

[0119] The purpose of this selection is because ΔN1 is obtained by solving the first two or more groups of moduli, so it meets the following conditions:

[0120]

[0121] So we set Γ that satisfies (22) 2L+1 From [Γ a1 ,Γ a2 ....Γ aL ] and [Γ b1 ,Γ b2 ....Γ bL ] and replace one of them at random to obtain the intersection of the estimated deviation values ​​under the new modular combination. If {ΔN0(a)}∩{ΔN'0(b)}={0,1,2}, it belongs to the first case, which is also the most expected modular combination. This modular combination can be directly used for untangling, that is, the current congruence equation group after replacing one congruence equation is used for untangling. If {ΔN0(a)}∩{ΔN'0(b)}={0,1,2,ΔN2}, ΔN2>>2, similarly, according to formula (23), the size of its ΔN2 satisfies the following:

[0122] ΔN2=k' 2L+1 *Γ 2L+1 +Δq 2L+1 +Δp 2L+1 (twenty four)

[0123] And because Δg 2L+1 ≠Δq 2L+1 +Δp 2L+1 , so k 2L+1 and k' 2L+1No matter what value is taken, ΔN1≠ΔN2, so the final estimated deviation value is only {0,1}, so by replacing the new Γ 2L+1 This can solve the problem of ΔN>>2 in the estimated deviation value. In this way, the remainder q can be obtained by modulus grouping. i The accurate reconstructed target value can be obtained by inverting all target estimation value sets obtained when the error does not exceed 1, and then finding the intersection of all sets and taking the average.

[0124] This method increases the robustness of the target value inversion using the Chinese Remainder Theorem by increasing the modular redundancy. a ,Γ b ) is greater than the maximum target value, the fewer modules in each group, the better, reducing computational complexity. Therefore, compared to search-based robust CRT methods, this method can maintain algorithm reliability and improve computational efficiency with fewer modules. Compared to closed robust CRT and multi-stage robust CRT, this method further expands the range of the remainder error bound to |τ| ≤ M, where M is the greatest common divisor of the same set of modules. Experiments on simulated and real InSAR data show that, compared with existing methods, the multi-channel InSAR elevation inversion method based on the grouped Chinese remainder theorem can improve the target value reconstruction success rate by approximately 20% to 30%, while maintaining the same signal-to-noise ratio. The root mean square error between the reconstructed target value and the true target value is also significantly lower than that of existing algorithms, demonstrating the effectiveness of this method. Notably, the deviation value is fixed once the module is given, eliminating the need for repeated calculation of the deviation value. Consequently, its computational complexity is similar to that of the classic CRT algorithm, resulting in high computational efficiency.

[0125] An embodiment of the present invention further provides a multi-channel InSAR phase unwrapping system, comprising:

[0126] A model building module is used to build a mathematical model of a group of congruence equations based on the basic principles of interferometric synthetic aperture radar;

[0127] A congruence equation group construction module is used to select multiple equal groups from the existing integer congruence equation groups consisting of K integer congruence equations corresponding to different moduli, and obtain multiple new integer congruence equation groups, where each new integer congruence equation group satisfies that the moduli are mutually prime.

[0128] The terrain data reconstruction module is used to calculate the deviation value set that may appear in the search process for each new set of integer congruence equations based on the acquired multi-channel InSAR data, and to find the intersection of multiple deviation value sets and determine whether the intersection is {0, 1, 2};

[0129] If so, a plurality of estimated value sets of terrain data are calculated based on the current multiple sets of new integer congruence equations, and the final terrain data is obtained by finding the intersection and taking the average value;

[0130] If not, select a new congruence equation with a different modulus to replace any congruence equation in the multiple new sets of integer congruence equations, and then recalculate and determine whether the intersection of the deviation value sets is {0,1,2}; if so, obtain a set of estimated values ​​of multiple terrain data based on the current multiple new sets of integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average; if not, obtain a set of estimated values ​​of multiple terrain data based on the multiple new sets of integer congruence equations constructed by the congruence equation group construction module and the current multiple new sets of integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average.

[0131] It should be understood that the functional unit modules in various embodiments of the present invention can be concentrated in one processing unit, or each unit module can exist physically separately, or two or more unit modules can be integrated into one unit module, and can be implemented in the form of hardware or software.

[0132] An embodiment of the present invention further provides an electronic terminal, including:

[0133] a memory having a computer program stored thereon;

[0134] The processor is configured to load and execute the computer program to implement the multi-channel InSAR phase unwrapping method as described above.

[0135] An embodiment of the present invention further provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the multi-channel InSAR phase unwrapping method as described above is implemented.

[0136] It can be understood that the same or similar parts of the above embodiments can be referenced to each other, and the contents not described in detail in some embodiments can refer to the same or similar contents in other embodiments.

[0137] The solution of the present invention is further described below with reference to experiments.

[0138] In the first experiment, the parameters selected for the simulation are 3≤L≤10, M=40, where Γ1 to Γ 10They are {5, 7, 11, 13, 17, 19, 23, 29, 31, 37}, etc. The target integer N (terrain height or terrain shape) is a real number uniformly distributed in 0≤N≤M(Γ-1). 10,000 experiments were conducted on different classical algorithms such as search-based robust CRT, closed robust CRT (CFR-CRT), and multi-stage robust CRT (MP-CRT). Figure 2 Figure 2 compares the TFR (Test Failure Rate, defined as the percentage of pixels where the difference between the reconstructed target value and the actual target value exceeds the range [-M, M]) performance of four robust CRT algorithms under the same noise conditions. The noise level follows a normal distribution.

[0139] Experimental data show that under the same noise intensity, the method of the present invention (New-CRT) has a lower experimental failure rate. That is, for the terrain height estimate obtained by inverting any pixel, the higher the probability that the true target value is the same as the integer part of the terrain height estimate, the smaller the deviation between the inverted terrain and the actual terrain.

[0140] In the second experiment, the peak function provided by Matlab was used to simulate the interference pattern. The size of the interference pattern was 256×256. Under different noise conditions, the method was compared with classic algorithms such as closed robust CRT and multi-stage robust CRT to verify the untangling effect of the peak interference pattern. Figure 3 (a), (d), and (g) are 256×256 interference patterns simulated by the peak function, and (b), (e), and (h) are the terrain height inversion results for the residual noise size within the residual error bound under the method of the present invention, closed robust CRT, and multi-stage robust CRT, respectively. (c), (f), and (i) are the difference images between the estimated and actual terrain height values, respectively. Among them, (a), (b), and (c) are the methods discussed in the present invention to obtain terrain height inversion data. When the size of the noise added to N is positive or negative M, the size of the noise added to N0 is positive or negative 1; (d), (e), and (f) are the closed robust CRT algorithm to obtain terrain height inversion data. When the size of the noise added to N is positive or negative M / 4, the size of the noise added to N0 is positive or negative 1 / 4; (g), (h), and (i) are the multi-stage robust CRT algorithm to obtain terrain height inversion data. When the noise added to N is M max / 4, the noise added to N0 is plus or minus 1 / 4; where M max is the greatest common divisor after grouping.

[0141] The following is an actual example to illustrate the above process:

[0142] The first one is to assume that each selected modulus combination contains three coprime moduli. For example, the selected modulus combinations are [13, 17, 19] and [23, 25, 31] respectively. The above formulas (11) and (21) are used to calculate the two sets of congruence equations to obtain the sizes of all possible error deviation values ​​ΔN0(a) and ΔN0(b) in all search processes. ΔN0(a)∩ΔN0(b)={0,1,2} satisfies the first case. The two sets of congruence equations can be directly used to perform terrain data inversion.

[0143] The second type is that if the selected modulus combinations are [11, 17, 23] and [13, 19, 21] respectively, the above formulas (11) and (21) are used to calculate the two sets of congruence equations to obtain the size of all possible error deviation values ​​ΔN0(a) and ΔN0(b) in all search processes, where ΔN0(a)∩ΔN0(b)={0,1,2,988}. At this time, when these two sets of modulus combinations are used, In the case where it is impossible to determine which solution is the correct one, we can select Γ7 from the remaining K-2L modules to satisfy (22), which can determine the size of ΔN'0(b). We can obtain ΔN0(a)∩ΔN'0(b)={0,1,2,4369}, that is, {ΔN0(a)∩ΔN0(b)}∩{ΔN0(a)}∩{ΔN'0(b)}={0,1,2,4369}∩{0,1,2,988}={0,1,2}, which satisfies the conditions and can be used for terrain data inversion after adding a module Γ7. The third method is to replace the congruence equation to satisfy ΔN0(a)∩ΔN'0(b)={0,1,2}. The untangling process is the same as the first method and will not be elaborated on.

[0144] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0145] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0146] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0147] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0148] Although the embodiments of the present invention have been shown and described above, it will be understood that the above embodiments are illustrative and are not to be construed as limitations on the present invention. A person skilled in the art may change, modify, replace and modify the above embodiments within the scope of the present invention.

Claims

1. A multi-channel InSAR phase unwrapping method, characterized in that: The steps include: S1: Construct a mathematical model of congruence equations based on the basic principles of interferometric synthetic aperture radar; S2: Select multiple integer congruence equations from the existing K integer congruence equations corresponding to different moduli and divide them into multiple groups to obtain multiple integer congruence equations, and each integer congruence equation group satisfies that the moduli are mutually prime; S3: Based on the acquired multi-channel InSAR data, calculate the deviation value set that may appear in each set of integer congruence equations during the search process, and find the intersection of multiple deviation value sets, and determine whether the intersection is {0, 1, 2}. If so, proceed to step S4; otherwise, proceed to step S5; S4: Calculate a plurality of estimated value sets of terrain data based on the current plurality of integer congruence equations, obtain the final terrain data by intersecting and averaging the values, and then end; S5: Select a new congruence equation with a different modulus to replace any congruence equation in the multiple sets of integer congruence equations, and then recalculate and determine whether the intersection of the deviation value set is {0,1,2}; if so, go to step S4; if not, integrate the multiple sets of integer congruence equations in step S2 and the current multiple sets of new integer congruence equations as the current multiple sets of integer congruence equations, and go to step S4.

2. The multi-channel InSAR phase unwrapping method according to claim 1, characterized in that: In step S1, the mathematical model of the congruence equations is obtained by the following method: The terrain data is represented by N, and the following real congruence equation model can be established; the terrain data is the terrain height or terrain shape variable; N=n i *M i +r i ,i=1,2,....,L Where n i is a fuzzy number, M i is the fuzzy height or fuzzy deformation, r i The interference phase is expressed as The height value or deformation value corresponding to the time; represents the winding phase; L represents the number of modules, i represents the i-th module or the i-th congruence equation; If M i There is a common greatest common factor M, and the above formula can be further simplified as: make Subtract r from both sides of the above equation c And divided by M, we can get the following integer congruence equation model: N0≡q i modΓ i Where, Γ i represents the modulus of the i-th congruence equation; The reconstruction formula of terrain data N is: N = M*N0+r c .

3. The multi-channel InSAR phase unwrapping method according to claim 2, characterized in that: In step S2, each set of integer congruence equations is expressed as follows: Where, represents the estimated value of N0 in the Ath group of integer congruence equations, m AL , Γ AL ,q AL , Δq AL ,Δp AL represents the folding integer, modulus, remainder, remainder error, and complementary error of the Lth congruence equation in the Ath group of integer congruence equations, and Δq AL ∈{-1,0,1},Δp AL ∈{-1,0,1}.

4. The multi-channel InSAR phase unwrapping method according to claim 3, characterized in that: In step S3, the estimated value of N0 is obtained by using the classical CRT for the Ath group of integer congruence equations: In the formula, is the deviation value of N0; Γ A = Γ A1 * Γ A2 *....* Γ AL ; γ Ai = Γ A / Γ Ai ; When the intersection of the deviation value sets formed by the possible values ​​of ΔN0 of multiple sets of integer congruence equations during the search process is {0,1,2}, then enter step S4; if the intersection of the deviation value sets is {0,1,2,ΔN1}, ΔN1>>2, then enter step S5.

5. The multi-channel InSAR phase unwrapping method according to claim 4, characterized in that: Step S4 includes: Based on the current multiple sets of integer congruence equations, multiple estimated value sets of N0 are calculated; Obtaining a plurality of estimated value sets of terrain data according to a plurality of estimated value sets of N0 and a reconstruction formula of terrain data N; The final terrain data is obtained by finding the intersection of the estimated value sets of multiple terrain data and taking the average value.

6. The multi-channel InSAR phase unwrapping method according to claim 4, characterized in that: In step S5, the modulus of the new congruence equation selected satisfies the following conditions: ΔN1=k 2L+1 *C 2L+1 +Δg 2L+1 ,Δg 2L+1 ∈{3,4,....,Γ 2L+1 -3} Where k 2L+1 is the folding integer, Δg 2L+1 is the remainder.

7. The multi-channel InSAR phase unwrapping method according to claim 6, characterized in that: In step S5, if the intersection of the recalculated deviation value sets is {0, 1, 2}, multiple estimated value sets of N0 are calculated based on the current multiple new integer congruence equations; multiple estimated value sets of terrain data are obtained based on the multiple estimated value sets of N0 and the reconstruction formula of the terrain data N; the final terrain data is obtained by calculating the intersection of the multiple estimated value sets of terrain data and taking the average value; If the intersection of the recalculated deviation value sets is {0, 1, 2, ΔN2}, ΔN2>>2, multiple estimated value sets of N0 are calculated based on the multiple sets of integer congruence equations in step S2 and the current multiple new sets of integer congruence equations; multiple estimated value sets of terrain data are obtained based on the multiple estimated value sets of N0 and the reconstruction formula of the terrain data N; the final terrain data is obtained by calculating the intersection of the multiple estimated value sets of terrain data and taking the average value.

8. A multi-channel InSAR phase unwrapping system, characterized in that: include: A model building module is used to build a mathematical model of a group of congruence equations based on the basic principles of interferometric synthetic aperture radar; A congruence equation group construction module is used to select multiple integer congruence equation groups from the existing K integer congruence equation groups corresponding to different moduli and divide them into multiple groups to obtain multiple groups of integer congruence equation groups, where each group of integer congruence equation groups satisfies that the moduli are mutually prime. The terrain data reconstruction module is used to calculate the deviation value set that may appear in the search process for each set of integer congruence equations based on the acquired multi-channel InSAR data, and to find the intersection of multiple deviation value sets and determine whether the intersection is {0, 1, 2}; If so, a plurality of estimated value sets of terrain data are calculated based on the current multiple sets of integer congruence equations, and the final terrain data is obtained by finding the intersection and taking the average value; If not, select a new congruence equation with a different modulus to replace any one of the multiple sets of integer congruence equations, and then recalculate and determine whether the intersection of the deviation value sets is {0,1,2}; if so, calculate the estimated value sets of multiple terrain data based on the current multiple sets of new integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average; if not, calculate the estimated value sets of multiple terrain data based on the multiple sets of integer congruence equations constructed by the congruence equation group construction module and the current multiple sets of new integer congruence equations, and obtain the final terrain data by finding the intersection and taking the average.

9. An electronic terminal, characterized in that: include: a memory having a computer program stored thereon; A processor, configured to load and execute the computer program to implement the multi-channel InSAR phase unwrapping method according to any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the multi-channel InSAR phase unwrapping method according to any one of claims 1 to 7 is implemented.

Images