CN107783079B - Use L0Two-dimensional phase rapid unwrapping method of norm cost function - Google Patents
Use L0Two-dimensional phase rapid unwrapping method of norm cost function Download PDFInfo
- Publication number
- CN107783079B CN107783079B CN201710899261.5A CN201710899261A CN107783079B CN 107783079 B CN107783079 B CN 107783079B CN 201710899261 A CN201710899261 A CN 201710899261A CN 107783079 B CN107783079 B CN 107783079B
- Authority
- CN
- China
- Prior art keywords
- phase
- unwrapping
- norm
- cost function
- equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
Landscapes
- Engineering & Computer Science (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
The invention discloses a use L0The two-dimensional phase fast unwrapping method of the norm cost function comprises the following steps: the numerical implementation of the differential equation solving algorithm generally includes converting differential relations into differential relations corresponding to discrete data, and performing iterative solution on the differential equations, wherein each step of processing is performed on the basis of a processing result of a previous step during iterative processing, and a difference delta phi between two processing results becomes an important content of the solution. The use of L0The invention provides a two-dimensional phase rapid unwrapping method of a norm cost function, aims at the problem that the efficiency of processing a large data block by the conventional phase unwrapping method is low, analyzes the characteristics of the cost function by aiming at the conventional minimum norm unwrapping method, and provides a method which accords with L0The two-dimensional phase unwrapping method of the norm cost function is combined with a data blocking scheme, and the problem of full convergence of low-frequency errors in a linear differential equation is solved, so that the unwrapping processing efficiency is improved.
Description
Technical Field
The invention relates to the technical field of interferometric synthetic aperture radar (InSAR) phase unwrapping, in particular to a method for unwrapping phase by using L0A two-dimensional phase fast unwrapping method of a norm cost function.
Background
Phase unwrapping is an important step in relevant data processing in the research fields of synthetic aperture radar interferometry (InSAR), synthetic aperture sonar interferometry (InSAS), medical Magnetic Resonance Imaging (MRI), and the like, and has important significance for the research of phase unwrapping methods.
The conventional unwrapping method can be basically divided into two types of path guidance and optimization approximation, the basic principle of the path guidance method is to obtain related information from the interior of data according to certain rules and work out a phase calculation route for avoiding discontinuous boundaries, and the method comprises a branch cutting method, a quality diagram guidance method, a Flynn minimum discontinuity method, a network flow-based method and the like, and the method belongs to a local method. Optimization approachThe approximate method is to realize the approximation to the optimal solution through the extremum optimization process of the set global objective function, and the method usually adopts a norm form to define a cost function, which is also called a minimum norm method and comprises a least square method and an L methodpNorm method, etc., belonging to the global method.
The minimum norm unwrapping method is to use the norm of the gradient difference between an unwrapping phase and an wrapping phase as a cost function to carry out global statistics and obtain a global optimal solution through the minimization process of an objective function. The most classical approach is to use L2The least square method of norm is finally converted into Poisson equation, has the advantage of convenient calculation, and also has L2The norm is excessively smooth, and the detailed data information cannot be recovered. On the basis of the least square method, the weighted least square method is improved by introducing weights. Ghiglia et al introduce a generalized LPAnd (5) carrying out induction on a norm framework and a norm type model.
The difficulty of the minimum norm unwrapping method in research is how to select a proper cost function model, and the problem of low efficiency of iterative processing in practical application is also solved, so that the minimum norm unwrapping method is not ideal in practical aspect. On the basis of the previous research, the invention analyzes the characteristics of the cost function and provides a method for using L0The fast two-dimensional phase unwrapping method of the norm cost function is combined with a data blocking scheme, and the problem of full convergence of low-frequency errors in a linear differential equation is solved, so that the unwrapping processing efficiency is improved.
Disclosure of Invention
Technical problem to be solved
Aiming at the defects of the prior art, the invention provides a method for using L0The fast two-dimensional phase unwrapping method of the norm cost function has the advantages of high unwrapping efficiency and the like, and solves the problem that the existing phase unwrapping method is low in efficiency when processing large data blocks.
(II) technical scheme
In order to achieve the purpose, the invention provides the following technical scheme: use L0The fast two-dimensional phase unwrapping method of norm cost function includes the following steps：
The numerical realization of solving the differential equation algorithm firstly converts the differential relation into the differential relation corresponding to the discrete data, and when the iterative solution of the differential equation is carried out, each step of processing is carried out on the basis of the processing result of the previous step, the difference delta phi between the processing results of the two steps becomes the important content of the solution, and the unwrapping results of the nth step and the (n + 1) th step meet the requirement of the equation
The equation to be solved can be rewritten as
At this time, the right term of the solution equation is the part of the original winding phase information with the known solution removed, and the difference delta phi is obtainednAnd updating the unwinding result, and determining whether the next iteration is needed or not after the detection result is identified.
The iteration termination condition adopts the detection of residual points in the unwrapping residual quantity of the wrapping phase as a judgment basis, and the unwrapping residual quantity is expressed as
δψn＝ψ-φn
When data delta psinWhen no residual point exists in the data, the iteration can be terminated, and delta psi is usednIntegral unwrapping of any path and merging into phinAnd forming a final unwrapping result, otherwise, continuing iteration.
In combination with the data blocking strategy, the basic process of phase unwrapping is represented in a pseudo-code manner as follows:
given the winding phase psi, the following steps are performed in sequence:
winding phase block psii,jInitial unwrapping phase
Computing
Calculating the right term of the equation;
the PCG method solves the equation to obtain a solution
Updating phase
Calculating winding phase residual
Detecting residual points, if not, stopping calculation, otherwise, returning to the step b;
the residual quantity of winding phase is unwound and added
Finishing the block processing;
after all the blocks are finished, the adjacent period offset calculation and the data splicing are carried out,
the parts b and h are independent of each other during the block processing, and can be organized to realize parallel processing.
Preferably, the cost function of the minimum norm unwrapping method comprises the following assuming that the two-dimensional phase ψ is the phase
After phase winding, i.e. as a result
Wherein W is the winding function and phase unwrapping is to achieve recovery from psi
The process of (2):
wherein
Phi is the result of the unwrapping process,
for the true value of the phase,
for domain definition, the function g describes the influence degree of different gradients in statistics, called cost function, and the constraint on the gradients keeps the unwrapping result consistent in phase change to some extent,
assuming that the function g is continuously cut to be first order derivable, the equation satisfied by minimization can be obtained from the Euler-Lagrange equation
Wherein div is a divergence operator, and psi and phi are calculated under the condition that the phase data satisfies the Nyquist sampling theorem
Satisfy
Equation (2) can be rewritten as
The functions taking the L of the gradient respectively0、L1And L2In norm form, different norm forms are used in solutionDifferent strategies are respectively corresponded to during the winding process, and L is in the case of discontinuous phase0The norm method discriminates between continuous and discontinuous phase cases, and L2The norm method has strong smoothing effect, L1The norm method is easy to have step effect, and the ideal method is to adopt L0Cost function in norm form, but L0The norm definition does not satisfy the requirement of continuous derivation of the cost function, and therefore, a method conforming to L is required to be found0The cost function expression mode in the form of norm can be AND-ed in the two-dimensional case
The associated expression is decomposed into the tangential component u of the u-contourTTAnd normal component uNNTwo parts, according to the literature, can be (3) rewritten as
When in use
When the value is small, i.e. corresponding to the absence of phase discontinuity, the tangential component uTTAnd normal component uNNThe corresponding constraints should play a role at the same time, and it is assumed that the two constraints are equivalent at this time, the coefficients have the same value, and the coefficients are within the range
The lower limit takes the same positive value, which requires that
When in use
When the value is larger, the situation of phase discontinuity may occur correspondingly, in order to reduce the phase discontinuityInfluence of the need to attenuate the normal component u of the contourNNSpecific gravity of, maintaining tangential constraint, i.e. conditional
However, these two conditions cannot be satisfied simultaneously for L0The norm being defined by modifying the above formula
I.e. the coefficients of both components approach zero when the gradient is large, but the tangential component uTTIs important over the normal component uNNAnd (4) partial.
Preferably, said is based on L0The cost function of norm criterion includes the following content, according to the requirement of formula (5) for cost function g, the design conforms to L0A function model of the norm, where the specific form of the cost function is selected as
Completely meets the requirement of the formula (5),
substituting the cost function of the formula (6) into the formulas (2) and (3) to obtain the unwrapping equation based on the new cost function
To eliminate the scaling problem of the phase range, the above equation can be modified
Where alpha is a positive constant, alpha can be slightly increased slightly when the phase psi is larger in scale, and vice versa,
it is not easy to directly solve the equation of the elliptic differential equation (7), so that
When solving, b is firstly taken as a constant, and an equation is solved
Obtaining a phase phi, updating b with phi, resolving again, and iterating until a final result meeting conditions is obtained, wherein the condition number of a linear equation set of the sparse coefficient matrix is large, the solution is carried out by adopting a Krysolv subspace algorithm, and a Preconditioned Conjugate Gradient (PCG) method is used.
Drawings
FIG. 1 is a diagram of three norm forms of the cost function of the present invention;
FIG. 2 is a schematic diagram of cost functions of different scale coefficients according to the present invention;
FIG. 3 is a block unwrapping diagram illustrating the present invention;
FIG. 4 is a schematic diagram of the composite winding phase (a), the unwrapping result (b) and its corresponding discontinuous boundary (c) according to the present invention;
FIG. 5 is a diagram of the first 4 and last equation coefficients (a-e) and corresponding phase residuals (f-j) according to the present invention;
FIG. 6 is a graph of the invention for a block derived unwrapping result and a discontinuity boundary profile;
FIG. 7 is a schematic diagram of 14 × 6 block unwrapping processing of interferometric phase data of InSAR in accordance with the present invention;
FIG. 8 is a statistical chart comparing the related indicators of the unwrapping result of the present invention with those of the prior art.
(III) advantageous effects
Compared with the prior art, the invention provides a method for using L0The fast two-dimensional phase unwrapping method of the norm cost function has the following beneficial effects:
1. the use of L0The invention provides a two-dimensional phase rapid unwrapping method of a norm cost function, aims at the problem that the efficiency of processing a large data block by the conventional phase unwrapping method is low, analyzes the characteristics of the cost function by aiming at the conventional minimum norm unwrapping method, and provides a method which accords with L0The two-dimensional phase unwrapping method of the norm cost function is combined with a data blocking scheme, and the problem of full convergence of low-frequency errors in a linear differential equation is solved, so that the unwrapping processing efficiency is improved.
2. The use of L0The fast two-dimensional phase unwrapping method of the norm cost function can obtain a better unwrapping result and a smaller discontinuous boundary, and can achieve higher processing efficiency by adopting a proper block scheme and a parallel processing technology.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The invention provides a technical scheme that: use L0The two-dimensional phase fast unwrapping method of the norm cost function comprises the following steps:
the numerical realization of solving the differential equation algorithm firstly converts the differential relation into the differential relation corresponding to the discrete data, and when the iterative solution of the differential equation is carried out, each step of processing is carried out on the basis of the processing result of the previous step, the difference delta phi between the processing results of the two steps becomes the important content of the solution, and the unwrapping results of the nth step and the (n + 1) th step meet the requirement of the equation
The equation to be solved can be rewritten as
At this time, the right term of the solution equation is the part of the original winding phase information with the known solution removed, and the difference delta phi is obtainednUpdating the unwinding result, determining whether the next iteration is needed or not after the detection result is identified,
the iteration termination condition adopts the detection of residual points in the unwrapping residual quantity of the wrapping phase as a judgment basis, and the unwrapping residual quantity is expressed as
δψn＝ψ-φn
When data delta psinWhen no residual point exists in the data, the iteration can be terminated, and delta psi is usednIntegral unwrapping of any path and merging into phinForming a final unwrapping result, otherwise continuing iteration,
when a linear equation is solved, if the data volume is large, the problem of slow attenuation of low-frequency errors can be solved, which can cause slow iterative convergence and time-consuming calculation, the invention adopts a strategy of data blocking processing, the unwrapping processing is only carried out on phase data in small areas of each block, and after the completion, data splicing is carried out, the process is as shown in figure 3, because data blocking is carried out, high-frequency errors in the blocking areas are quickly attenuated in the iterative processing, the iterative calculation speed is accelerated, the unwrapping processing efficiency is improved, low-frequency errors existing among the blocking areas can be compensated by calculating the offset between adjacent areas during data splicing, and thus, adverse effects caused by the low-frequency errors can be avoided;
in addition, the data blocking processing has no dependency relationship for the operation of each block data, and can be performed independently, in this case, a parallel computing mode can be adopted, the computing capability of the processor is fully exerted, the efficiency of the phase unwrapping processing for large data volume is improved, it is noted that the setting of the size of the blocking area should not destroy the integrity of the phase discontinuity boundary, otherwise the unwrapping result may have regional errors,
in combination with the data blocking strategy, the basic process of phase unwrapping is represented in a pseudo-code manner as follows:
given the winding phase psi, the following steps are performed in sequence:
winding phase block psii,jInitial unwrapping phase
Computing
Calculating the right term of the equation;
the PCG method solves the equation to obtain a solution
Updating phase
Calculating winding phase residual
Detecting residual points, if not, stopping calculation, otherwise, returning to the step b;
the residual quantity of winding phase is unwound and added
Finishing the block processing;
after all the blocks are finished, the adjacent period offset calculation and the data splicing are carried out,
wherein, the b and h parts are independent of each other in the blocking treatment and can be groupedThe cost function for the weave-implemented parallel processing, minimum norm unwrapping method contains the following, assuming that the two-dimensional phase ψ is the phase
After phase winding, i.e. as a result
Wherein W is the winding function and phase unwrapping is to achieve recovery from psi
When the unwrapping result and the true value meet a certain criterion, the unwrapping process can be considered to reach an ideal state, and the unwrapping process criterion is converted into an optimization problem related to gradient difference:
wherein
Phi is the result of the unwrapping process,
for the true value of the phase,
for domain definition, the function g describes the influence degree of different gradients in statistics, called cost function, and the constraint on the gradients keeps the unwrapping result consistent in phase change to some extent,
assuming that the function g is continuously cut to be first order derivable, the equation satisfied by minimization can be obtained from the Euler-Lagrange equation
Wherein div is a divergence operator, and psi and phi are calculated under the condition that the phase data satisfies the Nyquist sampling theorem
Satisfy
Equation (2) can be rewritten as
When the functions take the L of the gradient respectively0、L1And L2In the norm form, the patterns corresponding to the three norms are as shown in fig. 1, and different norm forms correspond to different strategies in the unwrapping process, and L is a phase discontinuity0The norm method discriminates between continuous and discontinuous phase cases, and L2The norm method has strong smoothing effect, L1The norm method is easy to have step effect, and the ideal method is to adopt L0Cost function in norm form, but L0The norm definition does not satisfy the requirement of continuous derivation of the cost function, and therefore, a method conforming to L is required to be found0The cost function expression mode in the form of norm can be AND-ed in the two-dimensional case
The associated expression is decomposed into the tangential component u of the u-contourTTAnd normal component uNNTwo parts, according to the literature, can be (3) rewritten as
When in use
When the value is small, i.e. corresponding to the absence of phase discontinuity, the tangential component uTTAnd normal component uNNThe corresponding constraints should play a role at the same time, and it is assumed that the two constraints are equivalent at this time, the coefficients have the same value, and the coefficients are within the range
The lower limit takes the same positive value, which requires that
When in use
When the value is larger, phase discontinuity may occur correspondingly, and in order to reduce the influence of phase discontinuity, the normal component u of the isoline needs to be weakenedNNSpecific gravity of, maintaining tangential constraint, i.e. conditional
However, these two conditions cannot be satisfied simultaneously for L0The norm being defined by modifying the above formula
I.e. the coefficients of both components approach zero when the gradient is large, but the tangential component uTTIs important over the normal component uNNMoiety based on L0The cost function of norm criterion includes the following content, according to the requirement of formula (5) for cost function g, the design conforms to L0A function model of the norm, where the specific form of the cost function is selected as
The function is in
Is zero and follows
Becomes close to constant, graph and L0Norm substantially identical, as in the curve of fig. 4
Shown in the specification, and
completely meets the requirement of the formula (5),
substituting the cost function of the formula (6) into the formulas (2) and (3) to obtain the unwrapping equation based on the new cost function
To eliminate the scaling problem of the phase range, the above equation can be modified
Where α is a positive constant, when the phase ψ is large in scale, α can be slightly increased slightly, and vice versa, as shown in fig. 2, the specific value can be determined as the case may be,
it is not easy to directly solve the equation of the elliptic differential equation (7), so that
When solving, b is firstly taken as a constant, and an equation is solved
Obtaining phase phi, updating b by phi, resolving again, and iterating until a final result meeting conditions is obtained, wherein the condition numbers of a linear equation set of the sparse coefficient matrix are all larger, a Krysolv subspace algorithm is adopted for solving, and a Preconditioned Conjugate Gradient (PCG) method is used in the invention.
In addition, the use of L0The information such as specific experimental data related to the two-dimensional phase rapid unwrapping method of the norm cost function is as follows:
a machine used for experimental data processing is provided with a 2.7GHz logic 8-core CPU and a dual-channel 1600MHz8GB DDR3 memory, a 64-bit operating system is adopted, a multi-core parallel processing mode of a shared memory is used, and a linear equation is used for resolving a conjugate gradient algorithm adopting a semi-coarsening multi-grid precondition with good stability.
The experimental data comprises synthetic winding phase data and real InSAR interference phase data, the size of the synthetic data is 513 multiplied by 513 pixels, a vertical upper half-period sine strip is superposed in the center of a plane background, and then the continuous sine strip is subjected to phase winding processing and is changed into 7 sections after 6 times of truncation. The difficulty of the unwrapping processing of the data lies in that the discontinuous phase difference exceeds 3 periods, and the method for controlling the unwrapping result by depending on the differential relation between phase data has certain difficulty. And (3) taking alpha to be 0.003, carrying out unwrapping treatment on the whole without partitioning, after 8 times of iterative calculation, successfully recovering the continuity of the sine strip, and completely detecting the discontinuous boundary between the strip data and the background. The two-dimensional display of the raw composite data and the unwrapped result and its corresponding discontinuous boundary distribution data is shown in fig. 4, with the data unit being the phase period (2 pi).
The iterative calculation process is also a process of performing attraction and stretching from the winding phase to the unwinding result according to the constraint condition of a given equation, and the action strength is completely determined by the coefficient b. The values (a-e) of the coefficients b for the first 4 times and the last time and the corresponding changes (f-j) of the winding phase residuals δ ψ for each time are given in fig. 5.
It can be seen from the figure that the coefficient b has a larger absolute value at the winding phase residual and a smaller value at the discontinuous boundary position, and the position with the smaller absolute value has a larger value. The original phase information has stronger attraction in a region with larger value, the solution is attracted to the region close to the original phase information, successive approach is finally stabilized near a final unwrapping result, the residual quantity at the moment does not contain residual points any more, and the residual quantity is obtained by direct path integration.
The data is processed in a blocking mode, so that the processing difficulty can be reduced, the data convergence speed is increased, and problems can be caused by improper use. The synthetic winding phase data is divided into 2 × 2 and 3 × 3 data blocks, and the data blocks are combined after the data blocks are separated and wound. The resulting unwrapping results and the corresponding discontinuous boundary distribution are shown in fig. 6. It can be seen that both of these blocking schemes split the discontinuous boundary, the boundary line recovered after unwrapping is no longer complete, and a break occurs, such as the circled position in (b), but the unwrapping result (b) of 2 × 2 blocking can still recover the continuity of the sinusoidal band, which is also an acceptable unwrapping result. Using the processing result (c) of 3 × 3 blocks, the large discontinuous phase difference of the band in the middle block area cannot be recovered, so that the positive period misalignment of the band in the middle block area occurs, and this misalignment cannot be compensated during block splicing, so that the sinusoidal band is broken in the unwinding result, and finally the unwinding processing fails.
Therefore, when data blocking is performed, attention must be paid to the setting of the block size, and although the smaller the size, the faster the settlement becomes, it is necessary to avoid the above-described case where the discontinuous boundary truncation causes the local data cycle in the block to be shifted. The phase discontinuity in the wrapped phase is often due to signal noise or anomalies in the local features of the observed target, and the wide range of discontinuity boundaries is rare, so the problem is not very severe.
The real InSAR interference phase is adopted as experimental data, the data size is 1398 × 622 pixels, the data has a gentle change trend, noise interference at several positions is serious, and some positions are located at data boundary positions, the influence generated by the noise can be eliminated only by utilizing the constraint of data at one side, and the iteration in a larger area is time-consuming. The data is now partitioned on a 200 pixel and 100 pixel scale, i.e., using 7 x 3 and 14 x 6 schemes, respectively, and then separately partitioned for unwrapping. Due to the difference of the phase data quality in the block region, iterative computation is not needed in blocks without residual points, and the blocks with more residual points may need many iterations to complete unwrapping, so that the task allocation of parallel computation is suitable for adopting a dynamic allocation mode so as to more fully utilize resources.
After the unwrapping processing of all the block data is completed, the unwrapping results of the block data need to be spliced and synthesized into a complete processing result. Before data splicing, low-frequency information contained in the winding phase cannot be presented in the block region, and is ignored during the unwinding process, but represents a full-period phase shift between data blocks during the block splicing, and therefore compensation of the offset value is required. Due to the phase continuity characteristic of the unwrapping result, the phase difference value between the blocks is counted at the adjacent positions of the blocks, discontinuous boundary data mixed by accident are easily submerged by continuous phase data, and the phase deviation value between the blocks can be obtained through simple statistical analysis. Selecting one block as a reference block, sequentially calculating deviation values according to the adjacency relation, compensating according to the deviation values, and performing data splicing to obtain a complete phase unwrapping result. In fig. 7, there are shown the InSAR interferometric phase data (a), the data 14 × 6 block map (b), the block unwrapping result (c), the block offset value table (d), the splicing result (e) after offset compensation, and the discontinuous boundary distribution map (f) corresponding to the unwrapping result, and it can be seen that even if the inter-data is decomposed into data blocks of 100 pixels in size, the unwrapping result can still maintain good integrity.
In order to compare the method of the present invention with other unwrapping methods, Goldstein branch cutting method, Quality map Guided method based on variance (Quality Guided), Flynn minimum discontinuous method, SNAPHU method and PUMA method based on network flow and the unwrapping method proposed by the present invention are used to perform unwrapping processing on interferometric phase of InSAR using various methods, respectively. The quality of the unwrapping result and the unwrapping efficiency are respectively counted by using quantitative indexes, including the average value and the error statistics of the difference value between the unwrapped result and the original wrapped phase data, the length of the phase discontinuity boundary and the running time length, and the statistical results are listed in table 1.
It should be noted that the data type of the experimental processing is 8-bit floating point type data, and the effective number of the data is not more than 7 bits, so that when the related precision index reaches more than 10-6, no error can be considered. It can be seen from the table that the mean statistics of the rewinding phase difference values all reach more than 10 < -6 >, and the statistical value of the medium error index is more prominent in the Goldstein bisection method because the bisection method may cause the situation of islands or holes in the area with poor data quality. The statistic result of the discontinuous boundary length index is that the Flynn method is optimal, which is the most important principle of the unwrapping method, and secondly, the invention accords with L0The norm criterion method can obtain a smaller discontinuous boundary while considering data continuity, other methods such as a branch cut method and a quality guidance method are slightly poor, the worst method is SNAPHU, and noise in data causes a discontinuous boundary generated by phase deviation on corners. Finally, the branch and cut method is fastest, the quality guide method is next to the quality guide method, and the PUMA method is slowest, but the method provided by the invention has slightly different required running time under two blocking schemes. The smaller the data block size, the faster the processing speed, but to avoid destroying the integrity of discontinuous boundaries in the unwrapping result, the time consumption of the 14 × 6 blocking scheme is better than that of the SNAPHU method, which is equivalent to that of the Flynn method, even if the 7 × 3 blocking scheme is slightly shorter than that of the SNAPHU method, and the efficiency is better than that of the PUMA unwrapping method as a whole.
By integrating the statistical indexes in the table (fig. 8), it can be seen that the phase unwrapping method provided by the present invention can obtain a better unwrapping result and a smaller discontinuous boundary, and can achieve a higher processing efficiency by using a suitable blocking scheme and a parallel processing technique, and the advantage of both quality and efficiency is not possessed by the classical minimum norm method.
The invention has the beneficial effects that: the use L0Norm ofThe invention provides a two-dimensional phase rapid unwrapping method of a cost function, aims at the problem that the efficiency of processing a large data block by the conventional phase unwrapping method is low, analyzes the characteristics of the cost function by aiming at the conventional minimum norm unwrapping method, and provides a method which is in line with L0The two-dimensional phase unwrapping method of the norm cost function is combined with a data blocking scheme to solve the problem of full convergence of low-frequency errors in a linear differential equation so as to improve the efficiency of unwrapping processing, obtain a better unwrapping result and a smaller discontinuous boundary, adopt a proper blocking scheme and a parallel processing technology, and use L0The fast two-dimensional phase unwrapping method of the norm cost function can achieve high processing efficiency, the quality and the efficiency have the advantages that the classical minimum norm method does not have, and the method has good expansibility, and specific indexes refer to experimental statistical data.
Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that changes, modifications, substitutions and alterations can be made in these embodiments without departing from the principles and spirit of the invention, the scope of which is defined in the appended claims and their equivalents.
Claims (3)
1. Use L0The fast two-dimensional phase unwrapping method of the norm cost function is characterized in that: the method comprises the following steps:
the numerical realization of solving the differential equation algorithm firstly converts the differential relation into the differential relation corresponding to the discrete data, and when the iterative solution of the differential equation is carried out, each step of processing is carried out on the basis of the processing result of the previous step, the difference delta phi between the processing results of the two steps becomes the important content of the solution, and the unwrapping results of the nth step and the (n + 1) th step meet the requirement of the equation
The equation to be solved can be rewritten as
At this time, the right term of the solution equation is the part of the original winding phase information with the known solution removed, and the difference delta phi is obtainednUpdating the unwinding result, determining whether the next iteration is needed or not after the detection result is identified,
the iteration termination condition adopts the detection of residual points in the unwrapping residual quantity of the wrapping phase as a judgment basis, and the unwrapping residual quantity is expressed as
δψn＝ψ-φn
When data delta psinWhen no residual point exists in the data, the iteration can be terminated, and delta psi is usednIntegral unwrapping of any path and merging into phinForming a final unwrapping result, otherwise continuing iteration,
in combination with the data blocking strategy, the basic process of phase unwrapping is represented in a pseudo-code manner as follows:
given the winding phase psi, the following steps are performed in sequence:
winding phase block psii,jInitial unwrapping phase
Computing
Calculating the right term of the equation;
the PCG method solves the equation to obtain a solution
Updating phase
Calculating winding phase residual
Detecting residual points, if not, stopping calculation, otherwise, returning to the step b;
the residual quantity of winding phase is unwound and added
Finishing the block processing;
after all the blocks are finished, the adjacent period offset calculation and the data splicing are carried out,
the parts b and h are independent of each other during the block processing, and can be organized to realize parallel processing.
2. Use of L according to claim 10The fast two-dimensional phase unwrapping method of the norm cost function is characterized in that: the cost function of the minimum norm unwrapping method includes the following, assuming that the two-dimensional phase ψ is the phase
After phase winding, i.e. as a result
Wherein W is the winding function and phase unwrapping is to achieve recovery from psi
The process of (2):
wherein
Phi is the result of the unwrapping process,
for the true value of the phase,
for domain definition, the function g describes the influence degree of different gradients in statistics, called cost function, and the constraint on the gradients keeps the unwrapping result consistent in phase change to some extent,
assuming that the function g is continuously cut to be first order derivable, the equation satisfied by minimization can be obtained from the Euler-Lagrange equation
Wherein div is a divergence operator, and psi and phi are calculated under the condition that the phase data satisfies the Nyquist sampling theorem
Satisfy
Equation (2) can be rewritten as
When the functions take the L of the gradient respectively0、L1And L2In the case of the norm form, different norm forms are respectively associated with different strategies in the case of unwrapping, and in the case of phase discontinuity, L is0The norm method discriminates between continuous and discontinuous phase cases, and L2The norm method has strong smoothing effectWith, L1The norm method is easy to have step effect, and the ideal method is to adopt L0Cost function in norm form, but L0The norm definition does not satisfy the requirement of continuous derivation of the cost function, and therefore, a method conforming to L is required to be found0The cost function expression mode in the form of norm can be AND-ed in the two-dimensional case
The associated expression is decomposed into the tangential component u of the u-contourTTAnd normal component uNNTwo parts, according to the literature, can be (3) rewritten as
When in use
When the value is small, i.e. corresponding to the absence of phase discontinuity, the tangential component uTTAnd normal component uNNThe corresponding constraints should play a role at the same time, and it is assumed that the two constraints are equivalent at this time, the coefficients have the same value, and the coefficients are within the range
The lower limit takes the same positive value, which requires that
When in use
When the value is larger, phase discontinuity may occur correspondingly, and in order to reduce the influence of phase discontinuity, the normal component u of the isoline needs to be weakenedNNSpecific gravity of, maintaining tangential constraint, i.e. conditional
However, these two conditions cannot be satisfied simultaneously for L0The norm being defined by modifying the above formula
I.e. the coefficients of both components approach zero when the gradient is large, but the tangential component uTTIs important over the normal component uNNAnd (4) partial.
3. Use of L according to claim 20The fast two-dimensional phase unwrapping method of the norm cost function is characterized in that: based on L0The cost function of norm criterion includes the following content, according to the requirement of formula (5) for cost function g, the design conforms to L0A function model of the norm, where the specific form of the cost function is selected as
Completely meets the requirement of the formula (5),
substituting the cost function of the formula (6) into the formulas (2) and (3) to obtain the unwrapping equation based on the new cost function
To eliminate the scaling problem of the phase range, the above equation can be modified
Where alpha is a positive constant, alpha can be slightly increased slightly when the phase psi is larger in scale, and vice versa,
it is not easy to directly solve the equation of the elliptic differential equation (7), so that
When solving, b is firstly taken as a constant, and an equation is solved
Obtaining a phase phi, updating b with phi, resolving again, and iterating until a final result meeting the conditions is obtained, wherein the linear equation set of the sparse coefficient matrix has a large condition number, and is solved by a Krysolv subspace algorithm, and a Preconditioned Conjugate Gradient (PCG) method is used.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710899261.5A CN107783079B (en) | 2017-09-28 | 2017-09-28 | Use L0Two-dimensional phase rapid unwrapping method of norm cost function |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710899261.5A CN107783079B (en) | 2017-09-28 | 2017-09-28 | Use L0Two-dimensional phase rapid unwrapping method of norm cost function |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107783079A CN107783079A (en) | 2018-03-09 |
| CN107783079B true CN107783079B (en) | 2021-10-26 |
Family
ID=61434231
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710899261.5A Active CN107783079B (en) | 2017-09-28 | 2017-09-28 | Use L0Two-dimensional phase rapid unwrapping method of norm cost function |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107783079B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108615264A (en) * | 2018-05-09 | 2018-10-02 | 国家计算机网络与信息安全管理中心 | A kind of three dimensional Phase unwrapping method based on Duality Decomposition |
| CN109541593B (en) * | 2018-10-30 | 2022-04-12 | 北京航空航天大学 | Improved minimum cost stream InSAR phase unwrapping method |
| CN110068823B (en) * | 2019-04-02 | 2021-03-09 | 中国人民解放军海军工程大学 | Block parallel quality guide rapid phase unwrapping algorithm |
| CN110824474B (en) * | 2019-11-25 | 2023-04-21 | 南京邮电大学 | Rapid phase unwrapping method based on minimum balance tree |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2007039772A1 (en) * | 2005-10-06 | 2007-04-12 | Roke Manor Research Limited | Phase unwrapping algorithm for array calibration with signals of opportunity |
| CN103942095A (en) * | 2014-03-18 | 2014-07-23 | 中国科学院软件研究所 | Two-dimensional phase position unwrapping method based on heterogeneous accelerating platform |
| CN104808203A (en) * | 2015-03-03 | 2015-07-29 | 电子科技大学 | Multi-baseline InSAR phase unwrapping method by iterating maximum likelihood estimation |
-
2017
- 2017-09-28 CN CN201710899261.5A patent/CN107783079B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2007039772A1 (en) * | 2005-10-06 | 2007-04-12 | Roke Manor Research Limited | Phase unwrapping algorithm for array calibration with signals of opportunity |
| CN103942095A (en) * | 2014-03-18 | 2014-07-23 | 中国科学院软件研究所 | Two-dimensional phase position unwrapping method based on heterogeneous accelerating platform |
| CN104808203A (en) * | 2015-03-03 | 2015-07-29 | 电子科技大学 | Multi-baseline InSAR phase unwrapping method by iterating maximum likelihood estimation |
Non-Patent Citations (4)
| Title |
|---|
| InSAR相位解缠算法的比较;羋宁轩;《测绘与空间地理信息》;20160331;第39卷(第3期);全文 * |
| Large-Scale L0-Norm and L1-Norm 2-D Phase Unwrapping;Hanwen Yu等;《IEEE TRANSACIONS ON GEOSCIENCE AND REMOTE SENSING》;20170831;第55卷(第8期);全文 * |
| Probabilistic Cost Functions for Network Flow Phase Unwrapping;Gabriel F.Carballo等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20000930;第38卷(第5期);全文 * |
| 水下SAS数据建模及干涉测量关键技术研究;高建;《中国博士学位论文全文数据库信息科技辑》;20121015;第84页第1段至第108页第4段 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107783079A (en) | 2018-03-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107783079B (en) | Use L0Two-dimensional phase rapid unwrapping method of norm cost function | |
| CN106804091B (en) | A kind of cyclotron isochronous magnetic field shimming method and system | |
| CN111143989B (en) | Frequency adjustment amount calculation method, module, system, storage medium, and device | |
| Hall | On the non-linear evolution of Görtler vortices in non-parallel boundary layers | |
| CN111191783A (en) | Self-adaptive quantization method, device, equipment and medium | |
| CN104198463B (en) | Raman spectra pretreatment method and system | |
| Bellettini et al. | Stable CMC integral varifolds of codimension $1 $: regularity and compactness | |
| Xie et al. | A fourth-order kernel-free boundary integral method for implicitly defined surfaces in three space dimensions | |
| CN114756974B (en) | Wall surface distance calculation method considering object surface normal information | |
| CN114065114A (en) | Method and system for predicting metering error of capacitive voltage transformer | |
| Matyjasek | Quasinormal modes of dirty black holes in the effective theory of gravity with a third order curvature term | |
| Mustafa et al. | A new 4-point C 3 quaternary approximating subdivision scheme | |
| CN109765611A (en) | Seismic data interpolation method and device | |
| Luskin | An approximation procedure for nonsymmetric, nonlinear hyperbolic systems with integral boundary conditions | |
| Lurie | Numerical design of composite radiofrequency pulses | |
| CN103942095A (en) | Two-dimensional phase position unwrapping method based on heterogeneous accelerating platform | |
| Vetter et al. | Optimized GPU histograms for multi-modal registration | |
| CN108229071A (en) | Cutting performance degradation assessment method and system based on AR models Yu SVDD algorithms | |
| Boldyrev et al. | An approach to multidimensional nonlinear optimization | |
| CN114170276A (en) | Magnetic resonance brain image hippocampus registration method | |
| CN114355439A (en) | Seismic inversion method and system based on Lagrange operator dynamic change | |
| Han et al. | Dynamic scaling for low-precision learning | |
| Zhou et al. | Monotonicity analysis and the reduced gradient method in constrained optimization | |
| Bickley et al. | On finite-difference methods for the numerical solution of boundary-value problems | |
| Seok et al. | Stochastic Learning Equation using Monotone Increasing Resolution of Quantization |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |