《测绘学报》
建设与学术的桥梁，拉近与权威的距离。
链接复制，《测绘学报》抖音(抖音)关注！
[《测量学报》的个人主页]复制了这条信息，按副本抖音(抖音)打开，TA的更多作品##7NsBSynuc88##[抖音(抖音)密码]
本文内容来自《测绘学报》 2020年第8期，审查编号GS(2020)4062期。
相位解缠的CKF局部多项式系数递推估计法
谢先明 ,孙玉铮 ,梁小星,曾庆宁,郑展恒
桂林电子科技大学信息与通信学院, 广西 桂林 541004
收稿日期：2019-09-18；修回日期：2020-06-04
基金项目：国家自然科学基金（41661092；61961009）；广西自然科学基金（2018GXNSFAA281196；2016GXNSFDA380018）；认知无线电与信息处理省部共建教育部重点实验室（CRK170108）；广西无线宽带通信与信号处理重点实验室基金（GXKL06180102）
摘要：为从含噪干涉相位数据中估计出解缠相位，本文提出相位解缠的CKF局部多项式系数递推估计法。利用基于修正矩阵束模型的局部相位梯度估计算法（AMPM）来获取多项式系数中的梯度信息，在此基础上获得局部多项式系数初始值（即状态变量初值），最后利用容积卡尔曼滤波（CKF）算法递推估计多项式系数状态估计值，从而获得解缠相位。可根据干涉图条纹密度以及相位噪声情况，分别采用逐行（或逐列）扫描方式或质量图引导策略引导容积卡尔曼滤波器解缠干涉图缠绕像元。模拟样例与实测数据试验结果表明，与其他同类方法相比，本文算法能从噪声干涉图中获得更高的解缠精度。
关键词：相位解缠 局部多项式近似 容积卡尔曼滤波 梯度估计
引文格式：谢先明, 孙玉铮, 梁小星, 等. 相位解缠的CKF局部多项式系数递推估计法. 测绘学报，2020，49(8)：1023-1031. DOI: 10.11947.
阅读全文：
全文概述
相位解缠在合成孔径雷达干涉测量(InSAR)、核磁共振成像(MRI)与数字全息干涉技术(DHI)中具有重要意义，是各种干涉测量的关键步骤[1-2]。通常而言，实际获得的干涉相位被缠绕在(-π, π]的区间内，与解缠相位之间存在2π的整数倍差异，使得相位解缠成为一项不可或缺的任务。在无噪声条件下，基于相邻像元相位之间差值的绝对值不大于π的假设，通过对缠绕梯度积分即可实现相位解缠[3-5]，但实际干涉测量中，由于存在各种噪声源对原始干涉数据的干扰，直接进行梯度积分容易造成误差传播，难以获得精确的解缠相位。因此，许多具有噪声抗差性的相位解缠算法相继提出。
传统的相位解缠算法主要包括基于路径跟踪策略的算法、最小范数类算法和网络规划类算法等。基于路径跟踪策略的算法通过干涉图先验信息(主要包括残差点分布和像元质量)来设置解缠路径，尽量避免误差蔓延；但当干涉图存在严重噪声时，则难以设置合适的解缠路径，这将间接影响算法的效率与精度，其代表性的算法有枝切法[6-7]、质量图法[8-10]。最小范数类算法通常利用解缠相位梯度与缠绕相位梯度之间的差值来建立代价函数，将相位解缠问题转化为最小范数框架下的全局优化问题。最小范数类算法在高效进行相位解缠的同时容易忽略干涉图细节信息，特别是在处理低信噪比干涉图时易出现曲面的过度拟合，其代表性算法包括最小二乘法[11-12]，L-p范数相位解缠方法[13]。基于网络规划的相位解缠算法着力于将解缠相位导数与缠绕相位导数之间差异的最小化问题转化为求解网络规划中的最小费用流问题，来提高干涉图相位解缠效率，但由于无偏的相干系数和其他权重因子难以确定，算法的精度难以保证，其代表性算法有最小费用流算法[14-15]。传统相位解缠算法在处理高信噪比干涉图时通常可以取得很好的效果，但当干涉图信噪比较低时，往往相位解缠结果较差。文献[16]将卡尔曼滤波技术应用到相位解缠中来，在相位解缠的同时进行滤波；文献[17—20]提出基于路径跟踪策略的无迹卡尔曼滤波(UKF)相位解缠算法；文献[21—22]提出基于二阶相位导数状态变量估计的卡尔曼滤波相位解缠算法。
近年来，文献[23—24]将局部多项式近似(LPA)技术应用到相位解缠中，通过求取局部窗口内多项式系数最优解来获得解缠相位(“PhaseLa”算法)。在局部多项式估计的框架下，将多项式系数作为状态变量，通过选取合适尺寸的局部窗口，文献[25—26]采用扩展卡尔曼滤波算法(LAEKF)和一种自适应局部窗口的线性卡尔曼滤波算法(LALKF)对多项式系数进行估计，在信噪比较高的干涉图相位解缠试验中获得了较好的效果。但是对于干涉图的高频率区域，LAEKF与LALKF需扩展局部窗口尺寸来弥补测量样本不足所带来的噪声敏感性，不仅增加了算法运行时间，而且扩大窗口容易破坏干涉相位梯度的局部稳定的假设，导致多项式系数估计失真。本文提出相位解缠的CKF局部多项式系数递推估计法，先用AMPM[18]对干涉缠绕相位进行处理，获取初始相位梯度，加快窗口内多项式系数收敛速度，从而可以选取较小窗口，提高相位解缠效率；再使用CKF算法对多项式系数进行状态估计，获得更精确的相位解缠结果。模拟样例与实测数据试验结果表明，与其他同类方法相比，本文算法能从噪声干涉图中获得更高的解缠精度。
1 局部多项式近似与AMPM算法
1.1 局部多项式近似
文献[23—24]采用局部多项式近似算法对干涉图解缠相位进行估计，设(m,n)像元为干涉图当前估计像元，在以(m,n)像元为中心的局部窗口内，相邻像元(ms, ns)处的相位值ϕL(ms, ns)可由待解缠像元相位值ϕ(m,n)与其对应的梯度近似表示
(1)
式中，L代表局部窗口尺寸，即窗口大小为L×L；ml和nl为局部窗口内相邻像元与当前估计像元的坐标差；c0=ϕ(m,n)，c1、c2分别为当前估计相位ϕ(m,n)行方向和列方向的一阶导数，也即干涉图(m,n)像元行方向和列方向的一阶相位梯度；Vl=[1 mlnl]T(T为转置操作)；C=[c0c1c2]T为干涉图局部多项式系数矢量。LAEKF[25]算法和LALKF[26]算法分别采用扩展卡尔曼滤波算法和基于自适应局部窗口的线性卡尔曼滤波算法对干涉图局部多项式系数矢量C进行估计。
1.2 修正矩阵束模型的局部相位梯度估计算法(AMPM)
在干涉图中，(m,n)像元的复干涉信号Z(m,n)可以表示为
(2)
假设以(m,n)像元为中心的局部窗口尺寸为(2κ+1)×(2κ+1)(本文中取κ=3)，在式(2)中，fy和fx分别表示局部窗口中行方向和列方向的局部频率，ℓ(m,n)表示附加在复干涉信号Z(m,n)上的噪声；ϕy和ϕx分别表示局部窗口内行方向和列方向的单位相位梯度。因此，估计干涉图中像元(m1,n1)与像元(m2,n2)之间的相位梯度方法如下。
把以(m1,n1)像元为中心的局部窗口内所有复干涉像元看作为一个复数矩阵Z(m1,n1)
(3)
式中，Z(i,j)κ=Z(i+m1,j+n1)，对Z(m1,n1)进行矩阵奇异值分解表示如下
(4)
式中，矩阵ΩZ中的σi(i=1, 2, …, 2κ+1)为计算得到的奇异值；UZ和VZ分别表示Z(m1,n1)的左奇异矩阵和右奇异矩阵。
对矩阵奇异值分解后，舍去较小奇异值可达到去噪的效果，据此构造降噪之后的复干涉相位矩阵Z(m1,n1)*如下[18]
(5)
式中，[0]2κ×(2κ+1)为2κ×(2κ+1)的零矩阵；σ*为奇异值的主要部分。
从Z(m1,n1)*中创建子矩阵Y0、Y1和Y2
(6)
式中，Y0表示取Z(m1,n1)*前2κ行和列；Y1表示取Z(m1,n1)*后2κ行和前2κ列；Y2表示取Z(m1,n1)*前2κ行和后2κ列，即Y0和Y1以及Y2等3子矩阵尺寸为2κ×2κ。对Y0做奇异值分解，并构建Y0*、Y1*和Y2*
(7)
式中，U00和V00分别表示Y0的左奇异矢量和右奇异矢量，包含了奇异矩阵的主要元素；U01和V01则为余下的奇异矩阵。
计算相位梯度估计值g(m1,n1)|(m2,n2)
(8)
式中，y1+和y2+分别代表Y1*和Y2*的伪逆矩阵。conj(·)表示取复共轭运算；Arg[·]代表取复数相角；g(m1,n1)|(m2,n2)表示像元(m1,n1)与(m2,n2)之间的梯度信息；ϕy－和ϕx－分别为局部窗口内行方向和列方向的单位相位梯度估计值，即式(1)局部多项式系数矢量C中的梯度c1和c2。
2 相位解缠的CKF局部多项式系数递推估计法
2.1 CKF递推估计多项式系数矢量
为了有效评估干涉图局部多项式系数矢量C，特别当干涉图信噪比较低时，LAEKF与LALKF通常需通过增大当前待估计像元局部窗口尺寸来提高算法的噪声抗差性，不仅增加了算法运行时间，而且扩大窗口容易破坏干涉相位梯度局部稳定的假设，导致相位解缠精度严重下降。本文把AMPM算法与CKF算法结合起来，提出相位解缠的CKF局部多项式系数递推估计法，先用AMPM算法获取干涉图初始相位梯度，加快窗口内多项式系数收敛速度，从而可以选取较小窗口，提高相位解缠效率；再使用CKF算法对多项式系数矢量C进行状态估计，获得更精确的相位解缠效果。
设(m,n)像元为当前估计像元，利用当前估计像元干涉相位与其对应一阶导数近似来表示相邻像元干涉相位ϕL(ml, nl)[23-24]，则C的观测值可表示为
(9)
式中，(ml, nl)像元为以干涉图(m,n)像元为中心的局部窗口内的像元；j(ml, nl)为状态变量的观测值；v1,l和v2,l分别为附加在复干涉信号虚部和实部的噪声。本文算法中系数矢量C中的相位梯度已经通过AMPM算法进行估计获得初始值，使C更快且精确收敛，据此本文算法采用的局部窗口尺寸L=3，提高算法运行效率的同时，又能维持相位梯度的局部稳定性，准确恢复干涉相位。
窗口内中心像元的系数矢量和估计误差方差初始化如下
(10)
式中，E[·]表示求取期望。在局部窗口内对系数矢量递推估计L2次，根据l(l依次取1, 2, …,L2)循环使用CKF算法估计系数矢量C，主要步骤为：首先通过状态转换矩阵F估计当前系数矢量，然后生成容积点Ci,l|l－1*以及预测值ci,l|l－1并计算对应的误差方差矩阵Pl|l－1
(11)
式中，F为3×3的单位矩阵； k为系数矢量维度的2倍；i[1]表示单位矩阵的第i列向量；Ql代表由于相位梯度估计不准确造成的过程噪声误差；根据文献[25]，P0|0和Ql的值通常可分别设置为单位矩阵I和5I。
利用Levenberg-Marquardt方法优化误差方差Pl|l－1并生成新的容积点Ci,l|l－1+
(12)
式中，优化误差方差时为避免直接对矩阵求逆带来的耗时问题，Pl|l－1－表示取Pl|l－1的对角元素组成的误差向量；Pl|l－1*表示优化后的误差向量；Pl|l－1+为根据Pll-1*构造的对角矩阵。算子diag(·)可以将输入向量生成方阵，方阵的对角线元素为输入向量的元素。μ表示优化Pl|l－1-的调节参数[19]。根据容积点计算相位观测值jl|l－1、像元的观测方差PJJ,l|l－1以及互相关协方差PCJ,l|l－1
(13)
式中，Rl为噪声引起的误差方差矩阵，估计卡尔曼增益并更新多项式系数矢量和误差方差
(14)
式(10)—式(14)为CKF状态估计主要过程。其中，jl*代表相邻像元干涉相位的观测值，Kl|l为卡尔曼增益矩阵。局部窗口内逐点采用CKF对系数矢量递推估计，局部窗口内所有像元扫描完成后，系数矢量C=CL2|L2包含最终解缠相位及修正后的相位梯度。
当前像元相位解缠完成后，按行(或列)方向扫描选取下一个解缠像元，采用CKF递推估计其多项式系数矢量，直至所有像元缠绕相位均被解缠。
2.2 质量引导策略
对于干涉条纹非常密集的缠绕相位图，即使含有很少的噪声，也容易造成条纹密集区域的解缠相位产生突变，引起整个解缠相位图的误差变大。若仍选用逐行(或列)选取像元策略极易造成误差传播，难以准确恢复解缠相位图。为此，针对条纹非常密集的干涉相位图，选用质量图引导策略[5]确定相位解缠路径，使本文算法能够更好地恢复解缠相位，本文所选质量图为相位导数方差图，具体定义为
(15)
式中，Δi,jx和Δi,jy为大小τ×τ(本文中取τ=3)窗口内的干涉相位偏导数； γm,n即为表征(m,n)像元的质量值，且γ越小代表当前像元处相位质量越好。
3 试验与分析
3.1 模拟多山地形干涉图
为分析算法在处理干涉相位变化复杂的干涉图数据时的相位解缠效果，模拟试验采用相位变化较为复杂的多山地形，采用Matlab中生成仿真干涉图，尺寸为300×300，原始解缠相位3D视图与原始解缠相位图如图 1(a)与图 1(b)所示，图 1(c)为信噪比为3.01 dB的含噪干涉图，图 2为PhaseLa、LAEKF(L=5)、LALKF和本文算法(采用逐行扫描像元策略)对图 1(c)干涉图解缠结果。
图 1 仿真干涉Fig. 1 Simulated interferogram
图选项
图 2 算法解缠结果对比Fig. 2 Comparison of algorithm results
图选项
表 1表示各算法处理不同信噪比干涉图时的均方根误差。图 3为不同算法处理8 dB至-1 dB信噪比干涉图时的解缠误差对比图，图 3忽略了PhaseLa、LAEKF(L=3)和LALKF算法的解缠误差变化曲线(算法在低信噪比条件下误差非常大，可认为该算法解缠失败)。从表 1和图 2可以看出，随着干涉图信噪比的降低，PhaseLa逐渐失效；LAEKF和LALKF把多项式系数当作整体状态变量进行估计，其中LALKF处理相位变化复杂且信噪比较低的干涉图时，解缠误差较大最终解缠失败，而LAEKF在低信噪比情况下仍能完整恢复解缠相位图，但与本文算法解缠结果相比，误差较大。通过图 3可以看出，与其他方法相比，本文算法在处理不同信噪比的干涉图时解缠误差均较小，尤其在处理低信噪比干涉图时仍保持极低的解缠误差。从表 2可以看出，本文算法耗时与LALKF、LAEKF(L=3)相比较长，与PhaseLa、LAEKF(L> 3)相当，即本文算法在精度，稳定性和效率方面均能以较好的性能处理低信噪比下的干涉图。
表 1 算法均方根误差对比Tab. 1 Comparison of algorithm root mean square error
| SNR/dB | PhaseLa算法 | LAEKF算法 | LALKF算法 | 本文算法 | |||
| L=3 | L=5 | L=7 | L=9 | ||||
| 3.01 | 0.13 | 0.15 | 0.10 | 0.13 | 0.13 | 0.45 | 0.10 |
| 1.42 | 3.82 | 0.18 | 0.12 | 0.16 | 0.15 | 1.53 | 0.12 |
| 0 | 8.01 | 1.95 | 0.30 | 0.18 | 0.18 | 2.83 | 0.15 |
| -1.00 | 14.17 | 4.15 | 0.33 | 0.21 | 0.23 | 3.99 | 0.18 |
表选项
图 3 算法均方根误差对比Fig. 3 Comparison of algorithm root mean square error
图选项
表 2 算法运行时间对比Tab. 2 Comparison of algorithm running time
| 干涉图尺寸 | PhaseLa算法 | LAEKF算法 | LALKF算法 | 本文算法 | |||
| L=3 | L=5 | L=7 | L=9 | ||||
| 300×300 | 28.9 | 9.9 | 29.7 | 32.6 | 41.9 | 8.6 | 29.4 |
表选项
3.2 条纹密集的干涉图模型
针对干涉条纹非常密集的缠绕相位图，逐行(或列)选取像元策略容易造成相位解缠误差传播，因此利用质量图引导策略来确定相位解缠的CKF局部多项式系数递推估计法路径，以避免相位误差的蔓延。图 4(a)为Long’s Peak干涉相位图，Long’s Peak干涉数据来源于文献[1]中的试验数据, 该干涉图是根据其数字高程模型产生的，大小为458×152像素，图 4(b)表示用来掩去受损区域的掩膜图，图 4(c)为本文算法在质量图引导策略下对干涉图 4(a)处理得到解缠相位图，图 4(d)为图 4(c)的重缠绕结果。由图 4可知，本文算法解缠结果较为平滑，且重缠绕相位图与原干涉图条纹基本一致，这表明本文算法有效解缠了上述条纹密集的干涉图。
图 4 Long’s Peak相位Fig. 4 Long's Peak interferogram
图选项
4 实测干涉图数据解缠结果
实测数据采用欧洲航天局哥白尼计划中的哨兵1号(Sentinel-1)2019年6月28号与2019年9月14号对伊朗Bam地区观测的两幅卫星图像生成的部分干涉图数据。图 5(a)为实测干涉相位图，图 5(b)为干涉图残差点。
图 5 实测干涉Fig. 5 Real interferogram
图选项
各算法针对图 5(a)的处理结果如图 6所示，其中，图 6(a)—图 6(d)为各算法对实测干涉图的解缠相位图，图 6(e)—图 6(h)为各算法解缠相位的重缠绕结果。PhaseLa算法利用高斯-牛顿方法计算得到权重矩阵，并引入梯度下降计算最优参数向量，但是当实测数据存在大量噪声时，难以迭代获取最优参数向量使损失函数达到最小值，恢复的解缠相位易出现大片断裂，其解缠相位严重受到噪声干扰，故其解缠结果较差，见图 6(a)和图 6(e)。LAEKF算法不能对噪声区域产生抑制以避免误差传播至其他区域，易形成条带状的误差蔓延，难以有效地解缠噪声实测干涉图，见图 6(b)和图 6(f)。尽管LALKF算法采用了质量图引导线性卡尔曼滤波器沿高质量区域到低质量区域的路径解缠干涉图缠绕像元，但难以准确估计低信噪比区域的梯度信息，故其在干涉图噪声严重的局部区域相位一致性较差，见图 6(c)和图 6(g)。图 6(d)为本文算法在质量图引导策略的指导下解缠图 5(a)的结果，其重缠绕结果见图 6(h)。从图 6(d)和图 6(h)可以看出本文算法解缠相位光滑与连续，不仅重缠绕相位图条纹与原始干涉图条纹一致，且有效去除原始干涉图中相位噪声，即本文算法处理实测数据时也取得很好的效果。
图 6 各算法相位解缠结果Fig. 6 Phase unwrapping results
图选项
5 结论
本文提出相位解缠的CKF局部多项式系数递推估计法，将AMPM梯度估计方法与局部多项式近似技术相结合，来获取多项式系数矢量中的梯度信息初始值，极大地加速了局部窗口内多项式系数状态向量的收敛速度，有利于提高算法相位解缠效率，最后利用容积卡尔曼滤波算法递推估计多项式系数，同时完成相位解缠与滤波。试验结果表明，与PhaseLa、LAEKF和LALKF等方法相比，在效率相当的情况下，本文算法具有明显的精度优势，尤其在处理低信噪比干涉图时，本文算法误差仍然较小。在处理噪声较大的实测InSAR干涉相位数据时，原始干涉图残差点统计个数共19 338个，占比为3.38%，采用质量图引导策略引导容积卡尔曼滤波器解缠干涉图缠绕像元，其解缠相位含有残差点共597个，占比减少至0.100%，表明本文算法在处理实测数据时具有良好的性能。
作者简介
第一作者简介：谢先明(1979-)，男，博士，研究员，研究方向为雷达信号处理，Email:xxmxgm@163.com
第二作者（通信作者）简介：孙玉铮（1994-），男，硕士生，研究方向为InSAR相位解缠算法研究，Email:814824985@qq.com
第三作者简介：梁小星（1993-），男，硕士生，研究方向为多基线InSAR高程反演研究，Email: 592154331@qq.com
第四作者简介：曾庆宁（1963-），男，博士，教授，研究方向为语音与图像信号处理、瞬变电磁反演成像，Email：qingningzeng@126.com
第五作者简介：郑展恒（1978-），男，硕士，高级实验师，信号与信息处理，glzzh@guet.edu.cn
团队简介
本论文作者依托桂林电子科技大学“现代信号与信息处理研究室”开展科学研究，目前本团队主要研究方向包括：（1）INSAR以及多通道INSAR技术应用；（2）语音信号处理及应用；（3）TEM探地成像技术应用。本研究室近几年以来培养已毕业研究生20余名，本研究室发表学术论文60余篇，申请发明专利20余项，其中受权10余项。近五年，本研究室主持省部级以上科研课题10余项，其中国家自然科学基金项目5项、广西自然科学基金项目5项。
《测绘学报（英文版）》（JGGS）专刊征稿：LiDAR数据处理
重磅 | 《测绘学报》主编杨元喜院士获“钱学森杰出贡献奖”
资讯 | 2020年国基金资助率创历年新低！相比2019年降幅或达10%！
资讯 | 《测绘学名词》（第四版）正式公布
权威 | 专业 | 学术 | 前沿
微信、抖音小视频投稿邮箱 | song_qi_fan@163.com
欢迎加入《测绘学报》作者QQ群： 751717395
进群请备注：姓名+单位+稿件编号
