CN110824474A - Quick phase unwrapping method based on minimum balance tree - Google Patents
Quick phase unwrapping method based on minimum balance tree Download PDFInfo
- Publication number
- CN110824474A CN110824474A CN201911164010.8A CN201911164010A CN110824474A CN 110824474 A CN110824474 A CN 110824474A CN 201911164010 A CN201911164010 A CN 201911164010A CN 110824474 A CN110824474 A CN 110824474A
- Authority
- CN
- China
- Prior art keywords
- value
- block
- phase
- point
- data
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 55
- 238000004364 calculation method Methods 0.000 claims abstract description 32
- 238000012545 processing Methods 0.000 claims abstract description 22
- 238000009826 distribution Methods 0.000 claims abstract description 8
- 238000004804 winding Methods 0.000 claims description 23
- 230000008569 process Effects 0.000 claims description 17
- 238000012360 testing method Methods 0.000 claims description 11
- 230000000694 effects Effects 0.000 claims description 5
- 230000015654 memory Effects 0.000 claims description 5
- 230000005540 biological transmission Effects 0.000 claims description 4
- 230000006870 function Effects 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 4
- 230000001902 propagating effect Effects 0.000 claims description 2
- 230000002159 abnormal effect Effects 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 7
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 6
- 238000005070 sampling Methods 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002457 bidirectional effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 230000012447 hatching Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
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
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
The invention discloses a rapid phase unwrapping method based on a minimum balance tree, which adopts a parallel processing mode, obtains data quality parameters by utilizing wrapping difference values in the horizontal direction and the vertical direction, respectively obtains positive and negative reliability graphs through positive and negative residual point sources based on the forward propagation of a data quality graph, searches a shortest phase discontinuous boundary by combining the propagation distribution of an abnormal residual point source, and iteratively constructs the reliability quality index distribution of a phase data angular point; and then converting the reliability index based on the angular point into an index of a data point, carrying out non-incremental back propagation on the basis of the index, constructing a reliability index map of an absolute phase calculation path, and carrying out phase unwrapping calculation by stages by using a flood flooding method on the basis of the reliability index map, so that the global optimal unwrapping processing can be realized, and all the calculation is based on a parallel mode.
Description
Technical Field
The invention relates to the technical field of InSAR data processing, in particular to a quick phase unwrapping method based on a minimum balance tree.
Background
Synthetic Aperture Radar (SAR) is a device that synthesizes radar observation data in the distance direction and the direction to form high-resolution large-Aperture radar observation data, and interferometric Synthetic Aperture radar (interferometric Synthetic Aperture radar) is a device that performs interferometric measurement to obtain spatial three-dimensional information of an observation target based on the SAR technology. Because the adopted radar signal has the advantages of penetrating clouds and fog, the InSAR can be used in some extreme environments which cannot use conventional observation means, and plays an important role in investigation and rescue of disastrous accidents. With the progress of sensor technology and the improvement of an observation platform, the space-time resolution of SAR observation data is higher and higher, the processed data volume is larger and larger, and a new requirement is provided for the efficiency of data processing.
Phase unwrapping is a key step of interferometric processing, different SAR data obtain interferometric data in a dual product mode, and phase difference information of the SAR data corresponding to the phase data part of the interferometric data is a key for solving target three-dimensional information. However, the interference phase is an auxiliary value that is periodically wound, and the absolute phase difference cannot be directly reflected. In the case that the data satisfy the sampling theorem, the data can be derived from the adjacent absolute phase according to the difference value smaller than the half period. The difficulty of phase unwrapping is the problem of local discontinuity of absolute phases caused by noise, interference and special conditions of an observation target, the discontinuity cannot be directly detected, a discontinuity edge needs to be directly or indirectly identified through different methods, and finally an unfolded phase result is obtained through calculation.
The existing unwrapping method mainly comprises a quality diagram guiding method, a minimum norm method, a network flow and the like. The quality map guiding method mainly obtains the index data reflecting the data quality distribution by obtaining the quality auxiliary information of the phase, and then calculates the low-quality data from the area with high data quality according to the quality distribution to realize the unwrapping of the phase data. The minimum norm method designs a discontinuous minimum norm rule according to continuous characteristics among data, and solves an absolute phase which meets an optimal criterion through optimization of the minimum norm. The network flow method mainly forms a network through the connection of residual points, calculates the flow of the connected network according to the corresponding cost, and obtains the network reflecting the discontinuous phase distribution by minimizing the overall cost, thereby realizing the phase conversation processing. This is also a currently more practical application of an unwinding method. For phase unwrapping processing of large data blocks, most of simple schemes for parallel processing are based on block independent processing of data, and a practical mainstream parallel unwrapping method basically does not exist.
Disclosure of Invention
The purpose of the invention is as follows: in order to overcome the defects of the prior art, the invention provides a quick phase unwrapping method based on a minimum balance tree.
The technical scheme is as follows: the invention provides a quick phase unwrapping method based on a minimum balance tree, which comprises the following steps:
(1) reading original winding phase data and transmitting the data to an equipment memory;
(2) according to a batch processing mode, decomposing the read winding phase data into regular data blocks, wherein the data blocks are processing units for parallel calculation;
(3) detecting the position of a winding phase residual error point, and taking the winding phase residual error point as a forward propagation source point;
1) traversing coordinates (x, y) of each target position in a two-dimensional data definition domain omega, performing closed-loop winding phase gradient integration at the position, and marking the position as a residual point when the integration is a nonzero value;
wherein ,
in order to wind the phase data on,is a gradient operator, W is a winding function, and s is a length differential quantity on a closed loop route;
(4) calculating a phase data quality propagation cost indicator Q
1) Traversing the definition position (x, y) in a phase definition domain omega, calculating a statistical index T reflecting the distribution continuity of phase data in a given domain of the phase definition domain omega, and counting the variance value of winding phase gradient in a neighborhood and the maximum value of the absolute value of the winding phase gradient in the neighborhood;
2) according to the calculation content of T, converting the T into a continuously forward increasing cost index value Q;
(5) setting a residual point as a source point of forward propagation of the reliability, setting an initial reliability index value of the source point as zero, and defining the reliability at (x, y) as:
wherein l is a minimum path from S to (x, y) such that
Obtaining a minimum value, wherein Q is a cost index value function, S is a source point, R is a source point set, and omega is a domain boundary;
(6) searching data block b to which point source belongs0Marking the block as a movable block;
(7) traversing all the movable blocks in parallel, and iterating to perform block convergence processing operation until all the movable blocks disappear;
(8) traversing all the positive and negative parameter difference points in parallel, and searching different-number parameter difference point pairs which are the closest reliability relationship to each other to construct a balanced tree;
(9) repeating the processes of 5-8 until no source point pairs with the shortest opposite-sign residual errors exist;
(10) searching the highest position of the reliability value, setting the priority value of the highest position as the priority value of the position path, using the highest position as a source point of backward propagation, setting the block as a movable block, and propagating according to the following formula, specifically as described in step 11;
P(l)＝min{P(u，v)|(u，v)∈l}
(11) calculating a priority value in a definition domain, traversing all the movable blocks in parallel and iterating the following steps until all the movable blocks disappear;
1) traversing each data definition point in the block, and iterating the following operations until each definition position in the block has a constant integral path priority value:
1.1) traversing (x, y) external adjacent points, and selecting the path with the highest priority value;
1.2) calculating the path priority value at (x, y) according to the path priority value of the point;
2) setting the active block state as a convergence state;
3) testing the transmission of a calculation path priority value to a neighborhood fast intermediate inactive block through the boundary of the four adjacent domains of the block;
4) setting the block state of the neighborhood block to be an active state when the priority value of the neighborhood block is changed;
(12) setting the phase of the highest priority position of the reverse path as a reference phase value as a source point of the phase back propagation expansion calculation, and setting the block as a movable block;
(13) calculating the lowest priority threshold value of the current expansion calculation according to the grading setting parameters;
(14) traversing each movable block in parallel to iterate the following processes until all the activities disappear quickly;
(15) increasing grades and returning to step 13 until the threshold value covers all definition domains;
(16) and finishing the phase unwrapping data processing, transmitting the unwrapping result back to the host and outputting the result.
The step (7) comprises the following steps:
1) traversing each data definition position (x, y) in parallel in the block, and iterating the following operations until each definition position in the block has a constant priority value:
a) traversing the (x, y) external adjacent points, and if the known reliability value exists, resolving the reliability value at the (x, y) position according to an Eikonal equation;
b) if the reliability value exists at the position (x, y), comparing the calculation result with the reliability value, and taking a smaller value to update, otherwise, directly updating by using the calculation result;
2) setting the processed active block state as a convergence state;
3) testing the transmission reliability value of the field boundary of the block to the inactive block in the field block;
4) the reliability value of the domain block is changed, and the block state is set to be an active state.
And (4) iteratively performing the following operations:
1) setting a negative (or positive) parameter difference point as a starting position, and marking a discontinuous boundary;
2) according to the positive (or negative) reference difference point forward propagation process, tracking reversely from the current position and moving to the forward position thereof, and marking the flow to increase progressively until the positive or negative reference difference point position is tracked;
3) canceling out the different-sign residual point pair to propagate the point source identity, resetting the reliability value of the propagation to be undetermined, setting the non-zero flow boundary to be superconducting and setting the corresponding propagation cost to be zero.
And (3) in step 11, calculating the path priority value at (x, y) according to the path priority value at the point in step (1.2):
1.2.1 if no treatment has been done at (x, y), update directly with the calculated value;
1.2.2 otherwise, comparing the new and old calculation priority values: if the new priority value is large, the old value is updated, otherwise, the old value is maintained.
Step (14) the process comprises the steps of:
1) and traversing each data definition point in the block, and iteratively performing the following operations until each defined position in the block has a constant phase value:
a) traversing (x, y) external adjacent points, and if known phases exist, selecting the path with the highest priority value;
b) calculating the unwrapped phase value at (x, y) according to the unwrapped phase value and the path priority value at the point;
2) setting the active block state as a convergence state;
3) testing the phase value transferred to the inactive block in the neighborhood block through the boundary of the four neighborhood regions of the block;
4) the neighboring blocks, with varying phase, set their block state to the active state.
Step b) calculating the unwrapped phase value at (x, y) according to the unwrapped phase value and the path priority value at the point, and specifically comprises the following steps:
1) if the (x, y) is not processed, directly updating by using a calculated value;
2) otherwise, comparing the new and old unwrapped phase values:
2.1) if the phase value is different, updating the phase value;
2.2) if there is no difference in phase value, the old value is maintained.
Has the advantages that: the invention discloses a quick phase unwrapping method based on a minimum balance tree, which has the following advantages:
1. the invention adopts a parallel propagation mode to realize bidirectional propagation resolving processing of a phase-offset point source, can effectively improve the calculation efficiency of large-data-block interference phase unwrapping processing, and can keep the global consistency of propagation by adopting a movable-block interaction mode;
2. the method for detecting the discontinuous phase edge by using the nearest reliability different-sign residual error point pair realizes dynamic marking through a boundary flow value and is beneficial to realizing complex multi-cycle discontinuous boundary detection by using the superconducting property of the boundary;
3. in the phase expansion calculation process, the calculation path priority value is constructed by using the reliability index distribution map, high-quality data with high reliability can be considered preferentially, and the reliability of phase calculation is ensured by adopting step-by-step calculation from high to low in a step-by-step mode.
Drawings
FIG. 1 is a graph of InSAR interferometric phase data used in an embodiment of the present invention;
FIG. 2 is a clockwise difference operator for residual point detection according to the present invention, wherein (a) is a horizontal first-order difference template, (b) is a vertical first-order difference template, (c) is a second-order difference template based on the graph (a), and (d) is a second-order difference template based on the graph (b);
FIG. 3 is horizontal direction quality description indicator data for an embodiment of the invention;
FIG. 4 is a schematic diagram of data points and corners and data blocks implemented in accordance with the present invention, wherein (a) is a schematic diagram of data points and corners; (b) is a data block diagram;
FIG. 5 is a phase data signature for an embodiment of the present invention;
FIG. 6 is a graph of initial positive residual reliability for an embodiment of the present invention;
FIG. 7 is an initial positive residual impact region of an embodiment of the present invention;
FIG. 8 is a final positive residual reliability map for an embodiment of the present invention;
FIG. 9 is a final positive residual impact area for an embodiment of the present invention;
FIG. 10 is a minimum balanced tree profile constructed in accordance with an embodiment of the present invention;
FIG. 11 is a graph of integrated reliability for an embodiment of the present invention;
FIG. 12 is a path integral precedence graph according to an embodiment of the invention;
FIG. 13 is the final unwrapping result for an embodiment of the present invention;
fig. 14 is a flow chart of a fast phase unwrapping method in accordance with the present invention.
Detailed Description
The technical scheme of the invention is further explained in detail by combining the attached drawings:
according to the flow chart of the fast phase unwrapping method based on the minimum balanced tree shown in fig. 14, the technical solution of the present invention is further explained;
technical implementation equipment parameters: configuring 1 CPU (4 cores, 2.7GHz), 1 GPU (384 stream processors, 0.71GHz and 2GB device memories) and 3 host memories (single 4GB and 1.6GHZ), running a 64-bit operating system, and testing by using C + + language codes.
Test objects: the real InSAR interference winding phase data is from test data provided by ROI _ PAC (repeat orthogonal interferometry PACkage), the data size (horizontal multiplied by vertical) is 1398 multiplied by 622, most data quality is better, the sampling data is continuous and strong, but local data has stronger noise interference, and the inside and the edge are distributed. The winding phase data is shown in its entirety in fig. 1.
For wrap data, the data structure is designed in 4 x 4 block sizes, including:
original winding data psi (float, 1400X 624, actual effective 1398X 622)
Winding horizontal differential psix(float, 1400X 624, actual effective 1397X 622)
Vertical differential twist psiy(float, 1400X 624, actually there isEffect 1398X 621)
Horizontal propagation differential qx(float, 1400X 624, actual effective 1398X 621)
Vertical propagation differential qy(float, 1400X 624, actual effective 1397X 622)
Multiple flag data flag (float, 1400X 624)
Positive active block list l+(int array)
Negative active block list l-(int array)
Block Mark data Block (char, 350X 156)
Positive reliability data p+(float, 1400X 624, actual effective 1399X 623)
Negative reliability data p-(float, 1400X 624, actual effective 1399X 623)
Positive residual source propagation region s+(int, 1400X 624, actually effective 1399X 623)
Negative residual source propagation region s-(int, 1400X 624, actually effective 1399X 623)
Discontinuous cycle count flow (char, 1400X 624, actually effective 1399X 623-1)
Integrated reliability data p (float, 1400X 624, actual effectiveness 1399X 623)
Corner priority data fg(float, 1400X 624, actual effective 1399X 623)
Data point priority data fc(float, 1400X 624, actual effective 1398X 622)
Unwrapped phase data
(float, 1400X 624, actual effective 1398X 622)
The flag data flag, bit flag, includes 1 bit of positive residual, 1 bit of negative residual paired with 1 bit of residual, 1 bit of discontinuous boundary point, 2 bits in positive reliability direction, 2 bits in negative reliability direction, 1 bit of boundary point; when the priority value is calculated reversely, the position of the unwound mark 1 and the unwound reference mark 1 can be repeatedly used. Discontinuous boundary cycle count is associated with direction, defining the direction from positive to negative horizontal to right and verticalAnd the lower part is positive, otherwise, the lower part is negative, 4 bits in the flow count horizontal edges, and 4 bits count vertical edges, wherein the ranges are-8-7. The block mark data block comprises 1 bit of positive and negative convergence marks, 1 bit of positive and negative convergence marks are listed in an active list, 1 bit of unwinding finish marks, 1 bit of priority value full coverage marks, and 1 bit of unwinding/flooding list active list marks. To save memory, the portions of data with non-overlapping cycles can be shared for multiplexing, e.g., p and p+/-And the like. After the winding phase data is read in psi, the winding difference value psi between adjacent samples in the horizontal direction and the vertical direction is calculatedx and ψy，
wherein
W(x)＝x-[x/2π]×2π
[. cndot. ] is a rounding operator. The differential templates in the horizontal and vertical directions are shown in fig. 2 (a) and (b), in which the hatching corresponds to the reference position. Then, the phase residual point is cumulatively detected according to the difference value in the clockwise direction in a 2 multiplied by 2 window,
of the 4 elements participating in the calculation, 2 are from the horizontal differential values and 2 are from the vertical differential values, and the calculation window template is shown in fig. 2 (c) and (d), shaded as the reference position. When the calculated value of r is not zero, namely a residual point exists, the corresponding position in the flag is marked. For the probing of the test data, a total of 16212 residual points were found to be present in the data.
Then, according to the following formula, the difference value variance in the horizontal direction and the vertical direction is respectively counted according to a 3 multiplied by 3 window,
therein take psix and ψyCorresponding to the difference values of the winding phase in the horizontal and vertical directions, respectively. Taking inverse variance as data quality drawingThe index is shown by the following formula
Where ε is a small constant to avoid the situation where the divisor is zero. Quality index data q in the horizontal direction obtainedxAs shown in fig. 3.
The data is partitioned according to the size of 4 × 4, the data points correspond to sampling positions, corner points correspond to angular positions between the data points, a schematic diagram is shown as (a) in fig. 4, a partitioning mode is shown as (b) in fig. 4, and 54600 blocks are totally partitioned. Respectively marking the blocks containing the residual error points and the blocks positioned at the boundary as the active states according to the positive and negative residual errors and the positions of the residual error points marked in the flag, and storing the active states in the blocks+ and l-In (1). Blocks of other blocks are marked as unprocessed. The position and boundary reliability values p of the residual error point+ and p-Initialized to zero and the reliability values at other locations are labeled as very large limiting constants. All elements of the dataflow are initialized to zero.
Forward propagation calculations of positive and negative reliability values are performed within the active block,
when the flow of the corresponding direction is nonzero, q of the corresponding side in the formula is zero. When the reliability value changes, the reliability value is updated, and the propagation direction at s are recorded in the flag+ and s+And marking a source point. Repeating the calculation until the reliability value in the active block does not change any more, marking the block as a convergence state in the block, marking the adjacent inactive block as an active block, and storing the active block in the block+ and l-In (1). Then checking one by one+ and l-And calculating the reliability value of the boundary neighborhood of the newly marked active block, and if the reliability value is not changed, cancelling the state of the active block by the block. Updating l according to block tag+ and l-The elements in (1) are processed in an iterative and repeated mode untilAnd marking all blocks to converge until the blocks mark, obtaining a positive and negative priority graph, and marking a forward propagation path. Fig. 5 is a path label of positive propagation of positive residual, and the gray values from 0 to 3 correspond to four different propagation directions, i.e., right, left, lower, and upper. Fig. 6 is a graph of the forward resulting propagation reliability for positive residuals. Fig. 7 is a propagation impact region label for positive residuals.
Searching the minimum residual point pair and the residual propagation influence area s by combining the positive and negative residual point propagation influence area maps+ and s+The minimum point should satisfy the following condition
Wherein u and v are a pair of residual point pairs, sv +Is the corresponding positive residual error point source at v, su -Is the corresponding negative residual point source at u. And searching discontinuous boundaries according to the propagation path flag, and correspondingly counting in the flow according to the direction. And simultaneously marking the propagation point source pairing marks of the residual points in the flag, and recalculating the reliability value, the propagation direction and the point source influence area of the paired residual point propagation areas according to the process.
And repeating the iteration of the processing, expanding the point pair balance tree until no residual point pair exists. Obtaining a final minimum balance tree and a final positive and negative reliability graph p+ and p-And a comprehensive reliability map p, the comprehensive reliability p being calculated as follows
FIG. 8 Final Positive residual reliability map p+. FIG. 9 is a final positive residual point source propagation impact area label diagram extracted in flag. FIG. 10 is a diagram of the minimum balancing tree resulting from the probes extracted in flow. Fig. 11 is a final integrated reliability map p.
On the basis of completing forward propagation and minimum balanced tree construction, a position where a maximum value is searched for p in the comprehensive reliability graph is set as a reference position, and a block is marked as an activity in a flag markBlock, listed in the active list l+。
Setting a current water level threshold value, repeatedly carrying out reverse flood inundation in the movable block by utilizing the discontinuous edge corresponding to the minimum balance tree in the comprehensive reliability graph p and the flow to calculate the unwrapping priority value,
fg-i，j＝min(pi，j，max(fg-i±1，j，fg-i，j±1))
when the priority value in the activity block is not changed any more, the block is marked to be in a convergence state, and when the threshold value is fully covered, the block is marked together. The block is marked as an active block adjacent to the new convergence block, and is listed as l+In (1). Then, the test is performed one by one and newly listed in+Block boundary priority value in (1), if no update is needed, then the activity flag is cancelled in block and from l+Is deleted.
Repeating the above process until l+Where no more active blocks are present. Then after changing the water level value, these processes are repeated until all blocks have completed the calculation of the priority value f and convergedgThe calculation of (2) is completed.
For fg and fcOnly the difference between the corner points and the data points can be realized by simple operations, such as fg-i，j＝max(fg-i±1，j，fg-i，j±1)
Processing to obtain a priority map fcAs shown in fig. 12.
Based on the priority value map fcAnd discontinuous boundary flow and reference position phase
Performing phase unwrapping, wherein the block with the reference position is marked as a movable block+The reference position flag marks the unwinding state.
Repeating the calculation in the active block according to the priority value water level threshold
Until all active blocks converge, set block flag to receiveThe convergence state is checked to see if all the blocks meet the water level threshold and marked, and then the adjacent blocks are set as active blocks and listed in+。
The above process is repeated until all the water level threshold blocks are met to complete the process. The water level threshold is adjusted and the processes are repeated until the unwrapping process of all blocks is completed and converges. Obtaining a disentanglement result
As shown in fig. 13.
As can be seen from the figure, the unwrapping data basically restores the phase continuity relationship, and the noise region processing is also more ideal. The whole treatment process takes 0.586 second, and the treatment is efficient. The method for rapidly phase unwrapping based on the minimum balance tree (hereinafter referred to as the method) is compared with the conventional SNAPHU method, and the comparison parameters are listed in Table 1, so that the norm statistical index of the method in a discontinuous boundary or the norm statistical index of the method in the discontinuous boundary can be seen
TABLE 1 comparison of the method with the SNAPHU method
| L0Norm/pixel | L1Norm/pixel | Run time/sec | |
| SNAPHU | 12989 | 12989 | 22.258 |
| Method for producing a composite material | 12597 | 12597 | 0.586 |
The above description is only an embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art can understand that the modifications or substitutions within the technical scope of the present invention are included in the scope of the present invention, and therefore, the scope of the present invention should be subject to the protection scope of the claims.
Claims (6)
1. A fast phase unwrapping method based on a minimum balance tree is characterized in that: the method comprises the following steps:
(1) reading original winding phase data and transmitting the data to an equipment memory;
(2) according to a batch processing mode, decomposing the read winding phase data into regular data blocks, wherein the data blocks are processing units for parallel calculation;
(3) detecting the position of a winding phase residual error point, and taking the winding phase residual error point as a forward propagation source point;
1) traversing coordinates (x, y) of each target position in a two-dimensional data definition domain omega, performing closed-loop winding phase gradient integration at the position, and marking the position as a residual point when the integration is a nonzero value;
wherein ,
in order to wind the phase data on,
for gradient operators, W is the winding function and s is the closed-loop pathLength differential of upper;
(4) calculating a phase data quality propagation cost indicator Q
1) Traversing the definition position (x, y) in a phase definition domain omega, calculating a statistical index T reflecting the distribution continuity of phase data in a given domain of the phase definition domain omega, and counting the variance value of winding phase gradient in a neighborhood and the maximum value of the absolute value of the winding phase gradient in the neighborhood;
2) according to the calculation content of T, converting the T into a continuously forward increasing cost index value Q;
(5) setting a residual point as a source point of forward propagation of the reliability, setting an initial reliability index value of the source point as zero, and defining the reliability at (x, y) as:
wherein l is a minimum path from S to (x, y) such that
Obtaining a minimum value, wherein Q is a cost index value function, S is a source point, R is a source point set, and omega is a domain boundary;
(6) searching data block b to which point source belongs0Marking the block as a movable block;
(7) traversing all the movable blocks in parallel, and iterating to perform block convergence processing operation until all the movable blocks disappear;
(8) traversing all the positive and negative parameter difference points in parallel, and searching different-number parameter difference point pairs which are the closest reliability relationship to each other to construct a balanced tree;
(9) repeating the processes of 5-8 until no source point pairs with the shortest opposite-sign residual errors exist;
(10) searching the highest position of the reliability value, setting the priority value of the highest position as the priority value of the position path, using the highest position as a source point of backward propagation, setting the block as a movable block, and propagating according to the following formula, specifically as described in step 11;
P(l)＝min{P(u,v)|(u,v)∈l}
(11) calculating a priority value in a definition domain, traversing all the movable blocks in parallel and iterating the following steps until all the movable blocks disappear;
1) traversing each data definition point in the block, and iterating the following operations until each definition position in the block has a constant integral path priority value:
1.1) traversing (x, y) external adjacent points, and selecting the path with the highest priority value;
1.2) calculating the path priority value at (x, y) according to the path priority value of the point;
2) setting the active block state as a convergence state;
3) testing the transmission of a calculation path priority value to a neighborhood fast intermediate inactive block through the boundary of the four adjacent domains of the block;
4) setting the block state of the neighborhood block to be an active state when the priority value of the neighborhood block is changed;
(12) setting the phase of the highest priority position of the reverse path as a reference phase value as a source point of the phase back propagation expansion calculation, and setting the block as a movable block;
(13) calculating the lowest priority threshold value of the current expansion calculation according to the grading setting parameters;
(14) traversing each movable block in parallel to iterate the following processes until all the activities disappear quickly;
(15) increasing grades and returning to step 13 until the threshold value covers all definition domains;
(16) and finishing the phase unwrapping data processing, transmitting the unwrapping result back to the host and outputting the result.
2. The method of claim 1, wherein the fast phase unwrapping method comprises: the step (7) comprises the following steps:
1) traversing each data definition position (x, y) in parallel in the block, and iterating the following operations until each definition position in the block has a constant priority value:
a) traversing the (x, y) external adjacent points, and if the known reliability value exists, resolving the reliability value at the (x, y) position according to an Eikonal equation;
b) if the reliability value exists at the position (x, y), comparing the calculation result with the reliability value, and taking a smaller value to update, otherwise, directly updating by using the calculation result;
2) setting the processed active block state as a convergence state;
3) testing the transmission reliability value of the field boundary of the block to the inactive block in the field block;
4) the reliability value of the domain block is changed, and the block state is set to be an active state.
3. The method of claim 1, wherein the method comprises: and (4) iteratively performing the following operations:
1) setting a negative (or positive) parameter difference point as a starting position, and marking a discontinuous boundary;
2) according to the positive (or negative) reference difference point forward propagation process, tracking reversely from the current position and moving to the forward position thereof, and marking the flow to increase progressively until the positive or negative reference difference point position is tracked;
3) canceling out the different-sign residual point pair to propagate the point source identity, resetting the reliability value of the propagation to be undetermined, setting the non-zero flow boundary to be superconducting and setting the corresponding propagation cost to be zero.
4. The method of claim 1, wherein the fast phase unwrapping method comprises: and (3) in step 11, calculating the path priority value at (x, y) according to the path priority value at the point in step (1.2):
1.2.1 if no treatment has been done at (x, y), update directly with the calculated value;
1.2.2 otherwise, comparing the new and old calculation priority values: if the new priority value is large, the old value is updated, otherwise, the old value is maintained.
5. The method of claim 1, wherein the fast phase unwrapping method comprises: step (14) the process comprises the steps of:
1) and traversing each data definition point in the block, and iteratively performing the following operations until each defined position in the block has a constant phase value:
a) traversing (x, y) external adjacent points, and if known phases exist, selecting the path with the highest priority value;
b) calculating the unwrapped phase value at (x, y) according to the unwrapped phase value and the path priority value at the point;
2) setting the active block state as a convergence state;
3) testing the phase value transferred to the inactive block in the neighborhood block through the boundary of the four neighborhood regions of the block;
4) the neighboring blocks, with varying phase, set their block state to the active state.
6. The method of claim 5, wherein the fast phase unwrapping based on a minimum balanced tree comprises: step b) calculating the unwrapped phase value at (x, y) according to the unwrapped phase value and the path priority value at the point, and specifically comprises the following steps:
1) if the (x, y) is not processed, directly updating by using a calculated value;
2) otherwise, comparing the new and old unwrapped phase values:
2.1) if the phase value is different, updating the phase value;
2.2) if there is no difference in phase value, the old value is maintained.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911164010.8A CN110824474B (en) | 2019-11-25 | 2019-11-25 | Rapid phase unwrapping method based on minimum balance tree |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911164010.8A CN110824474B (en) | 2019-11-25 | 2019-11-25 | Rapid phase unwrapping method based on minimum balance tree |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110824474A true CN110824474A (en) | 2020-02-21 |
| CN110824474B CN110824474B (en) | 2023-04-21 |
Family
ID=69558847
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911164010.8A Active CN110824474B (en) | 2019-11-25 | 2019-11-25 | Rapid phase unwrapping method based on minimum balance tree |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110824474B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104316922A (en) * | 2014-10-11 | 2015-01-28 | 南京邮电大学 | Regional division based multi-strategy InSAR (Interferometric synthetic aperture) radar phase unwrapping method |
| CN105158761A (en) * | 2015-08-31 | 2015-12-16 | 西安电子科技大学 | Radar synthetic phase unwrapping method based on branch-cut method and surface fitting |
| CN107783079A (en) * | 2017-09-28 | 2018-03-09 | 淮海工学院 | One kind uses L0The quick unwrapping method of two-dimensional phase of norm cost function |
| CN110490827A (en) * | 2019-08-30 | 2019-11-22 | 三亚中科遥感研究所 | The quick real-time processing method and system of airborne InSAR data |
-
2019
- 2019-11-25 CN CN201911164010.8A patent/CN110824474B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104316922A (en) * | 2014-10-11 | 2015-01-28 | 南京邮电大学 | Regional division based multi-strategy InSAR (Interferometric synthetic aperture) radar phase unwrapping method |
| CN105158761A (en) * | 2015-08-31 | 2015-12-16 | 西安电子科技大学 | Radar synthetic phase unwrapping method based on branch-cut method and surface fitting |
| CN107783079A (en) * | 2017-09-28 | 2018-03-09 | 淮海工学院 | One kind uses L0The quick unwrapping method of two-dimensional phase of norm cost function |
| CN110490827A (en) * | 2019-08-30 | 2019-11-22 | 三亚中科遥感研究所 | The quick real-time processing method and system of airborne InSAR data |
Non-Patent Citations (1)
| Title |
|---|
| JIAN GAO ET AL: "Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data", 《JOURNAL OF APPLIED REMOTE SENSING》 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110824474B (en) | 2023-04-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113014906A (en) | Daily scene reconstruction engine | |
| CN113781667B (en) | Three-dimensional structure simplified reconstruction method and device, computer equipment and storage medium | |
| Aslanyan et al. | The Knotted Sky I: Planck constraints on the primordial power spectrum | |
| CN110147815A (en) | Multiframe point cloud fusion method and device based on K mean cluster | |
| CN108007401A (en) | A kind of river and lake storehouse bank deformation detecting device and method based on boat-carrying InSAR platforms | |
| CN105654483A (en) | Three-dimensional point cloud full-automatic registration method | |
| CN102521870B (en) | Coloring reuse method for micro-polygon ray tracing | |
| Zhou et al. | GPS/terrestrial 3D laser scanner combined monitoring technology for coal mining subsidence: a case study of a coal mining area in Hebei, China | |
| Yu et al. | Large-scale multibaseline phase unwrapping: Interferogram segmentation based on multibaseline envelope-sparsity theorem | |
| CN113311433B (en) | InSAR interferometric phase two-step unwrapping method combining quality map and minimum cost flow | |
| Baselice et al. | Statistical edge detection in urban areas exploiting SAR complex data | |
| Arpaia et al. | Well balanced residual distribution for the ALE spherical shallow water equations on moving adaptive meshes | |
| CN113393577B (en) | Oblique photography terrain reconstruction method | |
| CN108761458B (en) | Morphological refinement-based interference SAR water body digital elevation model correction method | |
| CN103473749B (en) | A kind of method based on full variation image co-registration and device | |
| Zou et al. | WinIPS: WiFi-based non-intrusive IPS for online radio map construction | |
| CN110824474A (en) | Quick phase unwrapping method based on minimum balance tree | |
| Xiao et al. | Model-independent search for anisotropies in stochastic gravitational-wave backgrounds and application to LIGO-Virgo’s first three observing runs | |
| Gao et al. | Phase unwrapping method based on parallel local minimum reliability dual expanding for large-scale data | |
| Sohn et al. | Sequential modelling of building rooftops by integrating airborne LiDAR data and optical imagery: preliminary results | |
| Kitahara et al. | Algebraic phase unwrapping based on two-dimensional spline smoothing over triangles | |
| CN111007565A (en) | Three-dimensional frequency domain full-acoustic wave imaging method and device | |
| Chen et al. | An efficient phase unwrapping method based on unscented Kalman filter | |
| CN112835041A (en) | Multi-baseline InSAR elevation reconstruction method combining UKF and AMPM | |
| Naik et al. | A nonlinear iterative reconstruction and analysis approach to shape-based approximate electromagnetic tomography |
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 | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Gao Jian Inventor after: Shen Peipan Inventor after: Cheng Rui Inventor before: Gao Jian Inventor before: Shen Peishan Inventor before: Cheng Rui |