-
相位解缠是合成孔径雷达干涉(synthetic aperture radar interferometry, InSAR)和干涉合成孔径声纳(interferometric synthetic aperture sonar, INSAS)数字高程模型(digital elevation model,DEM)重构中极其关键的步骤之一,直接影响到所生成的DEM的速度和精度[1-2]。在所有相邻点相位差满足奈奎斯特采样定理的情况下,直接沿任意一路径进行积分就可以获得解缠相位。但由于InSAR和InSAS成像模式为斜视,会在成像结果中出现透视压缩、阴影和层叠,使得部分成像区域相干性低下,相位解缠过程变得非常困难[3-4]。随着缠绕相位图行列数的增加,相位解缠速度也成为影响干涉信号处理的一个重要因素,尤其是在实时干涉信号处理过程中[5-6]。
目前,相位解缠算法主要分为路径跟踪算法[7-8]、最小范数法[9]和网络流法[10-11]。路径跟踪的基本策略是通过显性或隐性的方式生成枝切线,阻止误差的传递。路径跟踪算法中枝切法效率较高,但枝切法解缠结果中容易存在“死区”。质量引导算法[12-13]解缠思想非常简单,完全依靠相位质量图来指导解缠路径,它隐含着质量图所指引的解缠路径不包含未平衡的残差点。最小范数算法[9]是将相位解缠问题转化为最小化全局函数优化问题,通过寻求一个全局的解缠相位来拟合观测的缠绕相位。网络流法[10-11]将相位解缠问题转化为求解最小费用流的网络优化问题,以达到降低相位解缠算法时间和空间复杂度的目的。在众多相位解缠算法中,质量引导算法是一种非常稳健的方法,它的核心包括两个方面,一是快速的高可靠相位质量图生成,二是快速的质量引导过程。在高可靠质量图计算方面,大量学者进行了深入研究[14-15],取得了较好效果。当前相位点质量值的计算都是以当前点为中心,在局部窗口进行的一系列数值运算。通常,可靠性越高的质量图,计算复杂程度越高,尤其当干涉相位图维数增大时,质量图的计算时间开销已成为影响相位解缠效率的重要因素。在快速质量引导过程优化方面,文献[7]提出了量化质量引导方法,显著加速了质量引导过程。但质量引导过程是一个典型的串行执行过程,影响着算法效率的进一步提升。
本文在多CPU多核平台上,充分利用共享内存环境下线程间高速通信的特点,对质量引导算法的并行化方法进行了研究,提出了共享内存环境下的质量引导相位解缠算法。首先,在满足减少重复计算和计算过程独立性的条件下,对质量图的计算方法进行优化;然后,在共享内存环境下实现了并行质量图求解。针对质量引导过程四邻域生长的独立性,给出了共享内存环境下的并行执行流程,在满足线程开启和关闭时间可忽略的情况下,可以优化相位解缠效率。最后,对InSAR和InSAS干涉数据进行解缠实验,验证了所提算法的高效性,能满足实时干涉信号处理需求。
-
相位解缠就是从缠绕相位数组中恢复真实相位差的过程。解缠前后缠绕相位ψ(x, y)和解缠相位φ(x, y)之间满足如下关系:
$$ \begin{array}{l} \psi \left( {x, y} \right) = W\left[{\varphi \left( {x, y} \right)} \right] = mod\left\{ {\varphi \left( {x, y} \right), 2\pi } \right\}\\ 0 \le \psi \left( {x, y} \right) < 2\pi \end{array} $$ 式中,W是缠绕运算;mod (·)表示求余操作。质量引导算法是一种路径跟踪相位算法,其解缠顺序完全取决于质量图,它是基于这样一个假设:质量图可以引导积分路径,这些积分路径所包围的区域不包括未平衡的残差点。质量引导算法的基本操作过程为:选择一个高质量的像素点作为初始点,检测它的4个邻近像素,将未解缠的邻近像素进行解缠并存入队列,依次从队列中移除最高质量像素,重复以上过程直至所有像素解缠完毕。为了提高质量引导算法效率,文献[7]提出了量化质量引导方法,取得了显著效果。但量化质量引导算法的解缠时间与干涉相位图的大小成正比,随着干涉相位图行列数的增加,解缠时间显著增加,影响解缠效率。
-
在进行相位解缠的时候,必须对干涉相位图的特征有充分认识,理想情况下,能够有一个二维标识数组来标识干涉相位图中相邻两点之间是否满足奈奎斯特采样定理,相位解缠时只需依据标识数组选择路径进行解缠。但在实际中,这种二维标识数组是无法获得的,而只能将配准后的相干系数近似地当作干涉相位两点之间满足奈奎斯特采样定理可能性的大小。在无法得到相干系数图时,直接从干涉相位图中提取伪相干系数图、相位梯度变化图或最大相位梯度图等来作为相位质量的度量标准。
相位质量图的计算过程如图 1所示,假设当前像素下标为(m, n),其质量值为qm, n,Am, n表示以点(m, n)为中心、k1和k2为半窗口大小的矩形数组,则当前点的质量值qm, n=f(Am, n)完全取决于局部窗口的缠绕相位值,f(·)表示质量值计算公式。如果对每个相位点采用同样计算方法,则可完成整个干涉图相位质量图的计算。计算质量图的过程中,通常需要计算行梯度dX和列梯度dY,计算一个点的相位质量值时,行梯度和列梯度的计算次数分别为2k1(2k2+1)和2k2(2k1+1)。对同一行的相邻相位点计算质量值时,需要重复计算的行梯度次数为(2k1-1)·(2k2+1),列梯度次数为2k2(2k1-1)。直接按照局部窗口内缠绕相位来计算相位质量值具有很好的并行性,但是存在大量冗余计算,而且冗余计算量与窗口大小成正比,严重影响相位质量图的计算效率。
在众多相位质量图中,相位梯度变化质量图是一种可靠相位质量图,其定义为:
$$ {\gamma _{m, n}} = \frac{{\sqrt {\sum\limits_{i = m-{k_1}}^{m + {k_1}} {\sum\limits_{j = n-{k_2}}^{n + {k_2}} {{{({\Delta ^x}_{i, j}-\bar \Delta _{m, n}^x)}^2}} } } + \sqrt {\sum\limits_{i = m - {k_1}}^{m + {k_1}} {\sum\limits_{j = n - {k_2}}^{n + {k_2}} {{{({\Delta ^y}_{i, j} - \bar \Delta _{m, n}^y)}^2}} } } }}{{(2{k_1} + 1) \times (2{k_2} + 1)}}{\rm{ }} $$ (1) 式中,Δi, jx=W[ψi+1, j-ψi, j]; Δi, jy=W[ψi, j+1-ψi, j]; ψi, j为缠绕相位; W为缠绕运算; ${\bar \Delta _{m, n}^x} $、${\bar \Delta _{m, n}^y} $分别表示在(2k1+1)×(2k2+1)窗口中Δi, jx、Δi, jy的均值。使用相位梯度变化质量图有一点需要注意的是,相位梯度变化质量图表征的是干涉相位数据坏的程度。
从图 1和式(1)中可以看出,相邻相位点在局部窗口计算质量图时仅在窗口边缘使用的相位梯度不同,并且这些不同梯度值对局部窗口的行均值${\bar \Delta _{m, n}^x} $和列均值$ {\bar \Delta _{m, n}^y}$产生影响。相位质量图可看作行梯度变化量与列梯度变化量之和,并且两者计算相互独立,为并行计算创造了条件。为避免行梯度与列梯度的重复计算,需开辟两块大小分别为M(N-1)和(M-1)N的空间分别存放行梯度和列梯度,行、列梯度的计算过程是完全独立的,直接使用parallel for指令对循环体内的内容实现并行计算。对于局部窗口中行梯度和列梯度变化量的计算,最重要的就是如何快速计算出窗口内的行均值和列均值。由于行梯度变化量与列梯度变化量两者的计算方式相同,本文只给出行梯度变化量的计算方法。
在图 2中,假设当前局部窗口内的行梯度变化量的和为Sr(m, n),其右相邻点的局部窗口内的行梯度变化量的和为Sr(m, n+1),则Sr(m, n+1)=Sr(m, n+1)-Sl+Sr,Sl表示最左一列行梯度的和,Sr表示最右一列行梯度的和,这样就消除了相邻窗口内相位梯度和的重复计算,进一步优化了行梯度变量的计算。完成局部窗口行梯度和的计算后,直接除以窗口大小(2k1+1)×(2k2+1)即可得到局部窗口的行梯度均值,然后利用行梯度变化公式,直接从内存中读取保存的临时行梯度,就可以完成行梯度变化量的计算。采用同样的计算方式完成列梯度变化量的计算,最后将两者相加即可得到最终的相位梯度变化质量值。
采用共享内存方法计算质量图时,为充分利用相邻相位点质量值计算的相同点以减小计算量,计算任务分配只在行方向进行。假设缠绕相位数组的行数为M,线程数为Tthread,当前行数为i的相位质量值采用线程mod (i, Tthread)计算。采用OpenMP实现质量图的并行计算, 主要步骤为:
1) 载入干涉相位图;
2) 分配行、列梯度存储空间;
3) 采用OpenMP中的parallel for指令并行执行循环体内函数,实现行、列梯度的并行计算;
4) 调用parallel for指令按照行进行计算任务分配,循环体内的内容分配到不同的线程执行,实现不同行相位质量值的并行计算。
-
完成质量图的并行计算后,下一步就是对质量引导过程进行优化。质量引导过程中, 每次是对当前相位点的四邻域进行操作,存在天然的并行性。在解缠算法引导过程中对邻域的操作主要包括计算解缠相位质量值和入队操作。为更好地满足并行性要求,使得四邻域的操作完全独立,可以采用4个优先队列,分别保存当前生长节点的上、下、左、右邻域节点,这样, 可以独立进行不同节点的入队操作。并行质量引导的基本流程如图 3所示,当前节点的四邻域节点分别解缠并存储到各优先队列时,各优先队列需要返回当前队列的最高质量值,4个队列完成入队操作后,比较各优先队列的最高质量值的大小,并记录最大质量值所对应的优先队列。如果4个优先队列都不为空,则让最高质量值所对应的节点出队列,并对其四邻域中未解缠的节点进行解缠并送入对应的优先队列。如果4个优先队列都为空,则解缠完毕。
在并行质量引导过程中,当前节点4个邻域中未解缠的节点并行完成解缠和入队操作后,需要进行同步来获取4个队列中的最高值节点,这是一个串行比较过程。如果4个优先队列都为空,则完成解缠; 否则, 获取当前4个队列中的最高质量节点,对其未解缠邻域重复进行解缠和入队操作,直至所有队列为空结束。在图 3所示的并行质量引导流程中,存在串、并行运算交替转换问题,这是质量引导过程中不可避免的,关键问题在于采用OpenMP进行串、并行环境切换时,需要重复进行并行环境的启动与释放,严重影响解缠效率。共享内存并行环境一次启动与释放的时间开销在毫秒级,并且并行环境启动与释放次数与总的相位点数成正比。因此, 只有在并行环境启动与释放时间可以忽略的情况下,才可以显著提升质量引导过程效率。
-
为验证共享内存环境下的质量引导相位解缠算法解缠效率,在如下环境对InSAR和InSAS干涉相位图进行了解缠实验:处理器Intel (R) Xeon (R) CPU X5650 2.67 G (2处理器12核), 内存48 G, 操作系统Windows 7专业版, 软件环境为VS2008。
实验选用了两幅干涉图,分别如图 4(a)和图 5(a)所示。图 4(a)是一幅InSAR干涉图,其大小为2 048×2 048像素。图 5(a)是一幅2010年7月在某内陆湖进行InSAS海试样机实验中获取的原始数据处理后的干涉图,其距离向点数为8 800,方位向点数为1 512,残差点为24 652,整块原始数据的获取时间为14 s。实验过程中, 相位质量图选用可靠性较高的相位梯度变化质量图,窗口大小设置为7×7像素,采用量化质量引导时,质量等级设置为1 000。测试程序运行时间时,均采用Release模式优化程序。实验过程中采用了原始质量引导算法[3]、量化质量引导算法[7]和共享内存质量引导3种质量引导方法,3种解缠方法差别主要体现在解缠效率,解缠结果是相同的。
图 4(a)所对应的相位质量图如图 4(b)所示,从质量图中可以看出干涉图中存在部分噪声区域,采用质量引导算法计算所得的解缠结果如图 4(c)所示,高质量区域解缠结果相位连续性很好,噪声区域存在一定的误差累积。图 5(a)所对应的残差点分布如图 5(b)所示,残差点主要分布在远距离处,给远距离的相位解缠带来了一定困难,最后解缠结果如图 5(c)所示,在远距离位置存在一定误差积累。低质量区域误差累积是质量引导算法无法避免的,但可以通过设计更可靠的相位质量图来对误差传播进行有效的抑制。
3种质量引导算法的效率比较如表 1所示。从表 1中可以看出,对InSAR和InSAS干涉图进行解缠,原始质量引导相位解缠算法需要时间最长, 分别为57 342 ms和647 640 ms,其时间主要消耗在质量引导过程。采用量化质量引导策略后,质量引导过程时间大幅度减少,分别从55 550 ms和642 320 ms减少到408 ms和1 547 ms,显著提高了解缠速度,总解缠时间分别为2 147 ms和6 929 ms,加速比分别达到26.4和93.5。采用量化质量引导后,解缠时间主要消耗在质量图计算步骤,尤其是对于InSAS干涉图,其质量图计算时间已达到5 382 ms,远超过其质量引导时间1 547 ms。为进一步提升相位解缠效率,采用共享内存方法对质量图计算进行优化后,两幅干涉图的质量图计算时间分别降为182 ms和689 ms,远低于其质量引导过程所消耗时间。采用共享内存方法后,两幅干涉图的解缠时间分别为603 ms和2 410 ms,加速比分别为95.1和268.7,提升了质量引导算法解缠效率。
表 1 不同质量引导算法解缠效率比较
Table 1. Efficiency Comparison of Different Quality-Guided Algorithms
干涉图 InSAR InSAS 数据大小 2 048×2 048像素 1 512×8 800像素 解缠方法 原始方法 量化质量方法 共享内存方法 原始方法 量化质量方法 共享内存方法 质量图计算/ms 1 792 1 766 182 5 320 5 382 689 质量引导/ms 55 550 408 421 642 320 1 547 1 721 总时间/ms 57 342 2 147 603 647 640 6 929 2 410 加速比 1 26.4 95.1 1 93.5 268.7 -
本文提出了一种共享内存环境下的快速质量引导算法,优化了质量引导相位解缠效率,满足了实时InSAS系统相位解缠需求。质量引导算法的解缠结果精度完全取决于质量图的好坏,通常质量图的精度与计算时间成正比,采用共享内存方法进行质量图计算优化后,可在有限时间内进行更复杂的高精度质量图计算,为在实时条件下提高相位解缠精度提供了条件。此外,本文分析并给出了共享内存环境下的质量引导过程优化方法,基本策略是采用4个优先队列同时对当前节点的四邻域节点进行生长,但目前受限于多线程启动与关闭的时间开销,体现不了并行引导优势。对InSAR和InSAS的相位解缠实验结果表明, 共享内存环境下的质量引导相位解缠方法提高了解缠效率,解缠时间远小于原始数据获取时间,满足了实时InSAS相位解缠需求。下一步主要在满足实时相位解缠的条件下对高质量相位质量图的生成方法进行研究,提高相位解缠精度。
-
摘要: 提出了一种共享内存环境下的并行质量引导相位解缠算法。首先,分析了相邻相位点质量值计算的内在关系,采用行、列数组保存中间结果消除梯度的重复计算;然后,按行进行计算任务分配,简化相邻窗口内的行、列均值计算;最后,通过OpenMP指令实现计算任务分配,完成质量图的并行计算。对于质量引导过程内在的并行性进行了分析,但受限于重复的线程启动和退出时间开销,难以体现速度优势。对合成孔径雷达干涉(synthetic aperture radar interferometry,InSAR)和干涉合成孔径声纳(interferometric synthetic aperture sonar,INSAS)干涉相位图进行解缠实验,结果表明,所提方法提高了相位解缠效率,为在实时条件下进一步提高相位解缠精度奠定了基础。Abstract: We proposed a parallel quality-guided phase unwrapping algorithm in shared memory environment. The intrinsic relationship between the neighboring phase points about quality value computing is analyzed firstly, and then the row and column arrays are used to store the intermediate results in order to eliminate the repeated computation of gradient. The computing task is allocated by row, which can simplify the computing of row and column gradient mean values. Finally, the allocation of computing task is realized using OpenMP instructions, and the quality map is computed in parallel. The inherent parallelism of quality guided process is also analyzed in-depth, but subject to the repetitive thread startup and exit overhead, it is difficult to reflect the advantage of speed. Unwrapping tests performed on InSAR and InSAS interferogram show that the proposed method greatly improves the efficiency of phase unwrapping, and provides a foundation for improving the solution precision under real-time conditions.
-
Key words:
- phase unwrapping /
- quality guided /
- parallel computing /
- shared memory
-
表 1 不同质量引导算法解缠效率比较
Table 1. Efficiency Comparison of Different Quality-Guided Algorithms
干涉图 InSAR InSAS 数据大小 2 048×2 048像素 1 512×8 800像素 解缠方法 原始方法 量化质量方法 共享内存方法 原始方法 量化质量方法 共享内存方法 质量图计算/ms 1 792 1 766 182 5 320 5 382 689 质量引导/ms 55 550 408 421 642 320 1 547 1 721 总时间/ms 57 342 2 147 603 647 640 6 929 2 410 加速比 1 26.4 95.1 1 93.5 268.7 -
[1] Carballo G F.Statistically-Based Multiresolution Network Flow Phase Unwrapping for SAR Interferometry[D].Sweden:Kungliga Tekniska Hogskolan, 2000 [2] Chen.Statistical-Cost Network-Flow Approaches to Two-Dimensional Phase Unwrapping for Radar Interferometry[D].California:Stanford University, 2001 [3] Ghiglia C D, Pritt D M.Two-Dimensional Phase Unwrapping:Theory, Algorithm, and Software[M].New York:John Wiley & Sons.Inc, 1998 [4] 许才军, 王华. InSAR相位解缠算法比较及误差分析[J].武汉大学学报·信息科学版, 2004, 29(1):67-71 http://ch.whu.edu.cn/CN/abstract/abstract4670.shtml Xu Caijun, Wang Hua. Comparison of InSAR Phase Unwrapping Algorithms and Error Analysis[J]. Geomatics and Information Science of Wuhan University, 2004, 29(1):67-71 http://ch.whu.edu.cn/CN/abstract/abstract4670.shtml [5] Eineder M, Hubig M, Milcke B.Unwrapping Large Interferograms Using the Minimum Cost Flow Algorithm[C]. International Geoscience and Remote Sensing Symposium Seattle, WA, USA, 1998 [6] Chen C W, Zebker H A.Phase Unwrapping for Large SAR Interferograms:Statistical Segmentation and Generalized Network Models[J].IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8):1709-1719 doi: 10.1109/TGRS.2002.802453 [7] 钟何平, 唐劲松, 张森, 等.利用量化质量图和优先队列的快速相位解缠算法[J].武汉大学学报·信息科学版, 2011, 36(3):342-345 http://ch.whu.edu.cn/CN/abstract/abstract499.shtml Zhong Heping, Tang Jinsong, Zhang Sen, et al. A Fast Phase Unwrapping Algorithm Based on Quantized Quality Map and Priority Queue[J]. Geomatics and Information Science of Wuhan University, 2011, 36(3):342-345 http://ch.whu.edu.cn/CN/abstract/abstract499.shtml [8] Zheng D L, Da F P.A Novel Algorithm for Branch Cut Phase Unwrapping[J].Optics and Lasers in Engineering, 2011, 49(5):609-617 doi: 10.1016/j.optlaseng.2011.01.017 [9] 陈强, 杨莹辉, 刘国祥, 等.基于边界探测的InSAR最小二乘整周相位解缠方法[J].测绘学报, 2012, 41(3):441-448 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201203023.htm Chen Qiang, Yang Yinghui, Liu Guoxiang, et al. InSAR Phase Unwrapping Using Least Squares Method with Integer Ambiguity Resolution and Edge Detection[J]. Acta Geodaetica et Cartographica Sinica. 2012, 41(3):441-448 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201203023.htm [10] Costantini M. A Phase Unwrapping Method Based on Network Programming[C]. Fringe'96 Workshop-ERS SAR interferometry, Zurich, 1996 [11] Carballo G F, Fieguth P W.Probabilistic Cost Functions for Network Flow Phase Unwrapping[J].IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(51):2192-2201 https://www.researchgate.net/publication/3202452_Probabilistic_cost_functions_for_network_flow_phase_unwrapping [12] 李芳芳, 占毅, 胡东辉, 等.一种基于质量指导的InSAR相位解缠快速实现方法[J].雷达学报, 2012, 1(2):196-202 doi: 10.3724/SP.J.1300.2012.20023 Li Fangfang, Zhan Yi, Hu Donghui, et al. A Fast Method for InSAR Phase Unwrapping Based on Quality Guide[J]. Journal of Radars, 2012, 1(2):196-202 doi: 10.3724/SP.J.1300.2012.20023 [13] 黄柏圣, 许家栋.一种基于新质量图引导的干涉相位快速解缠方法[J].系统仿真学报, 2010, 22(2):528-531 http://www.cnki.com.cn/Article/CJFDTOTAL-XTFZ201002059.htm Huang Baisheng, Xu Jiadong. Fast Phase Unwrapping Method for Interferometric Phase Based on New Quality-Guided[J]. Journal of System Simulation, 2010, 22(2):528-531 http://www.cnki.com.cn/Article/CJFDTOTAL-XTFZ201002059.htm [14] Cho B L, Kong Y K, Kim Y S.Quality Map Extraction for Radar Interferometry Using Weighted Window[J].Electronics Letters, 2004, 40(8):472-473 doi: 10.1049/el:20040294 [15] Quan C, Tay C J, Chen L, et al.Spatial-Fringe-Modulation-Based Quality Map for Phase Unwrapping[J].Applied Optics, 2003, 42(35):7060-7065 doi: 10.1364/AO.42.007060 -