引文格式：谢先明, 孙玉铮, 梁小星, 等. 相位解缠的CKF局部多项式系数递推估计法. 测绘学报，2020，49(8)：1023-1031. DOI: 10.11947/j.AGCS.2020.20190385.
阅读全文：http://xb.sinomaps.com/article/2020/1001-1595/2020-8-1023.htm
全文概述
相位解缠在合成孔径雷达干涉测量(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)进行矩阵奇异值分解表示如下
