一种基于相位差统计模型的InSAR图像相位解缠方法
【专利摘要】本发明涉及雷达技术领域，尤其涉及合成孔径雷达干涉测量(InSAR)图像处理中的相位解缠技术领域。本发明构造了一种基于相位差统计模型的InSAR图像相位解缠方法，首先利用相位差导数函数和地形高度变化量关系式，推导出InSAR相位差统计函数模型，利用相位差统计函数估计相位的变化顾及系统去相干因素和地形起伏变化对相位解缠的影响；在相位解缠中，利用相位差统计模型作为最小二乘相位解缠的非线性约束条件，并在模型解算中将解缠前后相位差的极值问题替换为解缠前后相位差的整数倍数差值问题，通过求解相位解缠前后整数倍数差最小值进行相位解缠。
【专利说明】
一种基于相位差统计模型的I nSAR图像相位解缠方法
技术领域
[0001] 本发明涉及雷达技术领域，尤其涉及合成孔径雷达干涉测量（InSAR)图像处理中 的
[0002] 相位解缠方法。
【背景技术】
[0003] 相位解缠是InSAR数据处理的关键步骤之一，其数据处理结果的精度直接影响到 最终数据产品的精度，然而相位解缠被称为NP问题，它受到两类相位去相干因素（相位噪 声）的影响：一类是处理系统本身引起的，例如热噪声、时间去相干、基线去相干（由于天线 视线不完全平行），以及图像配准过程中插值算法或图像滤波算法而引入的噪声；另一类是 地形的起伏变化造成相位欠采样引起几何相位畸变等，引起相干斑噪声。
[0004] 将两幅雷达图像处理后获取的干涉相位Φ的函数模型写成线性项Φ linear和非线 性项Φ ncm-lirmr和的形式，如公式(1 ):
[0006] (1)式中，Dii_r、D_-ii_r分别表示地形起伏变化引起的线性和非线性部分，v表 示线性地表形变速率，T表示干涉对获取的雷达数据时间基线，λ为雷达波长，R为斜距长，Δ h为待测点两次观测高度差，Β±为基线在垂直于视线方向的投影，即有效基线长度。公式（1) 中等号右边的前两项是线性部分，与垂直基线、时间基线成线性关系；第三项为非线性部 分，非线性部分造成传统枝切线解缠回路无法闭合以及解缠误差传递等现象，直接用线性 的方法无法获取高精度的解缠结果。
[0007] 根据是否顾及相位非线性，相位解缠方法分为线性和非线性相位解缠方法。目前， 线性相位解缠方法中最具代表性的是最小二乘(Least Squares，LS)解缠方法和加权最小 二乘(Weighted Least Squares，WLS)相位解缠方法。详见文献：Lu Y G，Wang X ZH,Zhang X P.Weighted least-squares phase unwrapping algorithm based on derivative variance correlation map. 2007。线性方法基于一个假设前提:缠绕相位的离散偏导数 (即邻近像元的相位差)呈现递增或递减的线性变化趋势，相位差的绝对值都小于I对离散 的偏导数沿任意方向迭代积分就可求得解缠后的相位。线性方法可通过相位加权来降低去 相干因素引起的相位质量降低对解缠的影响，在一定程度上，避免误差在全局相位图中的 传递，通过迭代的方法提高模型收敛的速度。但是，线性相位解缠算法不区分识别相位噪声 的来源，容易导致全局误差传递，不能避免由于地形因素引起的干涉条纹丢失或叠掩的情 况，部分相位得不到充分解缠，求解结果的精度不高;加权相位解缠单纯避免了相位残差点 对解缠区域的影响，不能避免由于地形因素引起的干涉条纹丢失或叠掩的情况，原因是影 响相位质量的去相干因素如相位几何畸变、相位噪声等使得相位解缠结果出现的偏差实质 上是解缠模型没有顾及非线性变化的相位而产生的偏差。
[0008] 目前，对于非线性的优化方法，实际上是构造非线性模型并求解的过程，也就是非 线性优化估计问题。基于优化方法的相位解缠在进行迭代求解时，要求预先估计解缠相位 的初值，再通过线性近似的方法获取算法或模型的最优解。详见技术文献l:Wang Q S， Huang H F,An X Y,et al. Improving phase Unwrapping Techniques by the Use of Nonlinear Phase Model 2011，2:Eineder M,Adam N.A maximum likelihood estimator to simultaneously unwrap 2005,中公开的技术内容。但是，非线性解缠模型在求解大像 素干涉图像时，使用线性逼近的迭代方法，无法提高运算的效率和解缠结果的精度，甚至有 时无法求解，并且在滤波或解缠中，没有直接顾及地形变化对相位差的影响，得到的结果并 不精确。
【发明内容】
[0009] 本发明构造了基于相位差统计模型的InSAR图像相位解缠方法，首先利用相位差 导数函数和地形高度变化量关系式，推导出InSAR相位差统计函数模型，利用相位差统计函 数估计相位的变化顾及系统去相干因素和地形起伏变化对相位解缠的影响；在相位解缠 中，利用相位差统计模型作为最小二乘相位解缠的非线性约束条件，并在模型解算中将解 缠前后相位差的极值问题替换为解缠前后相位差的整数倍数差值问题，通过求解相位解缠 前后整数倍数差最小值进行相位解缠。
[0010] 为实现上述
【发明内容】
使用了如下具体技术方案：
[0011] 步骤1、假定待处理干涉图像行、列像素大小为MXN，将其划分为kXk大小的局部 估计窗口，计算相位在行列两个方向的差值Α Φ，并对相位差在方位向和距离向求导，得到 相位相位差导数(离散相位下的相位频率）
[0012]
[0013]其中，Δ Χ?,」= χ?+1,」+Χ?-hjJxi,」，Δ ytFyijn+yi,」-：L-2yi,j 分别表示干涉图沿X、y 方向的相位差值;五(&)，忍(~)分别是kXk窗口内的x、y方向的相位差均值
是以窗口 中心的相位差导数标准偏差值。
[0014] 步骤2、根据上步相位差导数值，计算局部估计窗口内相位差统计模型，并估计局 部窗口内地形高度变化量的分布；
[0015] 步骤3、根据相位差频率与地形坡度角的关系式估计地形高度变化量；
[0016] 步骤4、获取局部估计窗口内相位差统计模型；
[0017] 步骤5、利用估计窗口内地形变化量分布函数对窗口内相位差进行平滑，在估计窗 口内将距离向和方位向相位差值乘以相位差统计模型，可得到窗口内顾及地形变化量的相 位差估计模型。
[0018] 步骤6、对解缠前后距离向和方位向相位差的整数倍数差值△ ky和△ kx进行估计并 积分，就可以得到解缠后的相位差整数，将求解相位差的问题转换为求解整数倍数差值的 问题，构造整数倍数参数的函数方程式为:Δ\. = [Δ么-Δ%]/2;r和Δ<=[Δ么-Δ%]/2;τ， 其中△ Φ4ΡΔ φχ表示相位点缠绕相位差，Δ%.和Ar表示解缠后相位差；
[0019]步骤7、将最小二乘相位解缠约束模型替换为整数倍数差的约束最小化模型，如公 式：
[0022] 其中，Ρ( Δ φχ)、Ρ( Δ φγ)表示相位点行向量、列向量相位差统计值，Δ Ρ与Δ k1^1 表示两次迭代的整数倍数差，ξ为极小值;约束条件s . t表示在局部窗口内，第t次和第t_l次 迭代求解时的解缠整数倍数差值满足差值最小的要求，K表示估计窗口的大小；
[0023] 步骤8、利用离散最小二乘迭代方法获取解缠相位。
[0024] 本发明技术方案带来的有益效果
[0025] (1)基于相位差统计模型的InSAR相位解缠方法，利用干涉相位差统计模型顾及地 形因素，在解缠中能够对地形变化起到预估计的目的。对于山区地形和地形突变的地形以 及地形相位中有沉降时，形成的干涉条纹包含混叠相位的干涉数据来说，基于相位差统计 特征探测地形引起的噪声，获取的解缠图像平滑，重缠绕结果更加接近原始干涉图，解缠精 度更高。
[0026] (2)本发明使用相位差统计值代替相位差值作为最小二乘相位解缠的附加约束条 件，附加约束条件能够对受噪声影响的相位差进行筛选，在保证相位平滑的同时消除噪声 对解缠积分路径的影响。
[0027] (3)在解缠模型求解中，针对处理复杂的Hessian矩阵和大矩阵收敛的问题，而且 要兼顾如何改变约束条件，在保证精度的条件下以求达到最优的解缠效果，求解比较困难， 效果不理想的问题，本发明将模型转换为附加约束条件的最小二乘相位解缠模型等同公 式，将求解相位差估计的问题转换为求解整数倍数估计的问题，利用相位差统计特征函数 估计每次迭代的结果，每一次迭代都体现了离散的最优化的步骤。
【附图说明】
[0028] 图1为不同入射角对相同坡度的影响图；
[0029] 图2为InSAR干涉相位差统计特征模型获取流程。
[0030] 图例说明，图1中①、②、③分别表示不同的侧视雷达入射角度，分别表示 测区地形坡度倾角;shadow为无雷达反射信号的阴影区。
【具体实施方式】
[0031] 下面结合附图对本发明技术方案进一步说明。
[0032] 为了方便本发明中技术描述，首先进行了以下定义：
[0033] 定义1:设干涉雷达两次观测斜距差为Δ R，两次观测相位差为Δ φ公式，两者的公 式如下，
[0036] 其中，λ为雷达波长，Ro、!^分别为两次斜距长，h为待测点高度;B±为基线在垂直于 视线方向的投影，即有效基线长度;B//为水平基线；（ XQ，yQ，Z())和(x，y，z)分别为两次观测的 坐标;T X、Ty为地面观测点在沿方位向和距离向坡度角，θ〇为两次观测雷达下视角的平均值。
[0037] 首先在附图1中：在坡度起伏变化较大的区域，入射角大于θ2时雷达图像将产生阴 影，随着入射角减小，阴影区也在收缩;在入射角小到θ 3时将出现顶底倒置的现象造成雷达 影像的几何失真。当入射角相同时不同的地形坡度对相位差的影响也是相同的。
[0038]定义2:地形高度变化量公式
[0039]对于观测点的地形起伏变化可利用其所在点高度变化量进行估计，而地形高度变 化量可认为其服从三个参数影响的对数函数分布：
[0041]其中，|d|表示地形高度变化量d在雷达视线向方向的投影(侧视雷达对视线向的 相位变化最敏感）：d= | d | cos0(0表示观测点切线与雷达视线向间的夹角），| d | = [tanT 为三个函数控制参数，针对不同的地形通过改变三个参数可以对 地形高度变化进行较好的拟合。
[0042] 对参数的说明以及取倌可参见下表：
[0044]具体技术方案步骤如下：
[0045]步骤1、基础数据的准备工作。
[0046]获取两幅SAR图像配准、滤波后的干涉图像，计算相干系数γ和离散干涉相位Φ， 包括附加的系统参数信息(波长、基线、下视角、斜距和如果可以获取的地形DEM信息）。在实 际图像处理过程中，
进行干涉相干值的估计， 5；式表示滤波后两幅同一地区的SAR图像51与&共辄相乘。式中1表示多视数。并将干涉图读 取为浮点形式的相位矩阵，并标记行列大小为M、N，单位为像素。
[0047] 步骤2、计算干涉图行、列方向局部估计窗口内相位差导数值。
[0048] 将干涉图像划分为kXk的估计窗口，计算相位在行列两个方向的差值△ Φ，并对 相位差在方位向和距离向求导，得到相位差导数(离散相位下的相位差频率）。在离散图像 中，相位导数可用相位差分来表示：
[0049]
[0050]其中，Δ Χ?,」= χ?+1,」+Χ?-1;厂2xi,j，Δ ytfyid+i+yi,」-Hyi.j分别表示干涉图沿行列 方向的相位差值;名(&)，￡(△·>')分别是kXk窗口内的行列方向的相位差均值； 口中心的相位差导数标准偏差值。
[0051] 步骤3、根据上述步骤2得出的相位差导数值，计算局部估计窗口内相位差统计模 型，并估计局部窗口内地形高度变化量的分布。
[0052] 根据地形高度变化量公式（1-3)获取测绘区地形高度变化量的分布。对由三个参 数模拟的地形高度值变化量d取对数并除他们的雅可比行列式，设相位估计窗口的范围区 间大小为kXk，d x、dy为地形高度变化量d在距离向和方位向两个方向的投影，那么局部估计 窗口内地形高度变化量的二维概率密度函数公式为：
[0055]通过对估计窗口内地形高度变化量密度积分可以得到其分布函数，设
，那么干涉相位图中地形高度变化量的离散分布函数Fd为变量△ t与d的函 数，公式如下：
[0057] 步骤4、根据相位差导数与地形坡度角的关系式估计地形高度变化量。
[0058] 在上式地形变化量分布函数中存在一个描述地形高度的变量d，由于地形高度的 变化与地形坡度角的变化是成正比例的关系，那么可以由相位差导数与地形坡度角的关系 式表示这个变量。即地形坡度变化与相位导数的关系，如下式(4)、（5):
[0061 ]
:，x、y分别表示干涉相位的距离向和方位向，Δ fx、Δ fy分别表示 两个方向的相位差导数;C为雷达波传播速度;B丄为基线在垂直于视线方向的投影，即有效 基线长度;λ为波长;R为斜距；Θ为入射角；τ为瞬时地面采样点的坡度角。将(4)、（5)式进行 逆变换得到公式(6):
[0063]其中，W=cos0(dy-tan0)/kB，上式说明相位差导数是以相位频率、地形坡度角和入 射角为参数的函数。在公式中除地形高度变化量的分量dx、dy两个未知量外，其余参数均为 系统已知参数，或通过前面的步骤求取的变量值。
[0064] 步骤5、获取局部估计窗口内相位差统计模型。
[0065] 由步骤3对测区地形高度变化量的二维概率密度函数公式(2)以及地形坡度角与 相位差导数的关系函数式(6)，推导出两个方向的相位差统计模型如下式：
[0068] 其中，1是比例常数，在估计相位模型中受估计窗口大小的影响，实际估计窗口越 小影响越小，可将估计窗口设置为3 X 3大小。Tee为估计窗口内采样比，这是一个受入射角Θ 和地形坡度角arctandy同时控制的量，两者差值越大说明干涉效果越差需要提高采样率， 使用一个反映两者变化的三角函数sin(0-arctand y)表示。
[0069] 步骤6、利用局部估计窗口内地形变化量分布函数对窗口内相位差进行平滑。
[0070] 根据相位差统计函数式(7)、（8)，在估计窗口内将距离向和方位向相位差值乘以 相位差统计模型可得到窗口内顾及地形变化量的相位差估计模型，如(9)、（ 10):
[0073]
表示估计窗口 k X k内相位差的估计值，且它们之间的关系满足下式的
条件：
[0075] 由相位差统计特征与相位差值的变换关系可知，可以利用相位差统计特征估计相 位差的变化趋势，识别由于地形造成的相位残差点。
[0076] 步骤7、构造基于相位差统计模型的InSAR图像最小二乘相位解缠模型。
[0077] 利用相位解缠前后离散相位差的平方和最小构建最小二乘相位解缠模型，具体如 公式（12):
[0080]其中，Μ和N分别为干涉图相位矩阵的列数和行数，Δ 是相位点（i，j)处的缠绕 相位差，Φυ是相位点（i，j)处的解缠后相位;ξ为极小值;乂^,,.,,(~7>^尺,)表示估计窗 口中相位差的最大估值，求解中使用迭代求解的方法，约束条件s.t保证在估计窗口内，第t 次和第t-1次迭代求解时的解缠相位相位差估值满足差值最小的要求，K表示估计窗口的大 小，要保证中心相位点（i，j)处差值最小，K窗口的方向和大小可按照如下的规则进行自适 应调整：当差值变大时再分别搜索（i+n，j)、（ i，j+n)、（ i+n，j+n)三个方向，η的取值保证搜 索窗口为一个k X k的窗口；若改变搜索方向不能满足条件，则改变窗口的大小：1 =[鮮2」， [^/2」表示上取整;当差值变小时K = 2 X K，表示将窗口扩大一倍继续搜索。
[0081 ]步骤8、将相位解缠模型转变为整数倍数差的估计问题。
[0082]结合InSAR离散相位矩阵的特性，在相位解缠中，若只考虑相位解缠前后相位差的 实数部分则相位解缠与相位缠绕是两个互逆的过程：
[0084]其中，Φ是缠绕相位，Φ是解缠后绝对相位，k为一个整数倍数值，那么相位(i，j)的 缠绕相位差和绝对相位差定义为：
[0086] 对于一幅MX N的图像：
[0087] (i, j)eSii = 0,l,... ,M-2； j = 0,l,... ,Ν-1
[0088] (i, j)es2i = 0,l,. . . ,M-1； j = 0,l,. . . ,Ν-2
[0089] 绝对相位差与缠绕相位差相差2π整数倍：
[0091]由上式可以求解缠绕相位差与解缠后相位差的离散导数偏差(除以2π后），构造解 缠前后距离向和方位向相位差的整数倍数&%和4^的函数方程式为：
[0093] 只要对Δ&',·和进行估计并积分就可以得到解缠后的相位，将求解相位差的问 题转换为求解整数倍数差的问题。
[0094] 利用预先计算得到的相位差统计特征函数，重新构造解缠整数倍数差的函数模型 如式(13):
[0096] 其中，为利用相位差统计模型估计后的重新计算得到的相位差整数倍 数，此估值顾及了相位差的变化情况，也就是顾及了地形因素的变化，得到的解缠模型更加 准确、可靠。
[0097] 步骤9、利用解缠前后整数倍数差值最小构建解缠模型。
[0098] 将最小二乘相位约束解缠模型替换为整数倍数差的约束最小化模型，如公式 (14)：
[0100]其中，约束条件s.t保证在估计窗口内，第t次和第t-l次迭代求解时的解缠整数倍 数差值满足差值最小的要求，K表示估计窗口的大小；在离散干涉矩阵中，
?这一解缠模型表不最 小化相位差统计特征值，即用解缠相位差与预估缠绕相位差的偏差应尽可能为零。
[0101]通过公式（14)得到离散的相位差偏差量是无旋的绝对相位的相位差场，再沿路径 积分即可以得到解缠相位。积分公式如（15):
[0103] 其中，ent( ·)表示取整运算。
[0104] 步骤10、模型的解算方法。
[0105] 利用离散最小二乘迭代的方法求解解缠公式：
[0106] Pk；；； = (Pk；+lJ + Pk：_Xj +Pk：^ +PklJ_l-Pi)lA
[0107] 上式表示第t+1次相位值是由第t次的相位变化得到的，设（i，j)为估计窗口中心 像素，PU可以由下式得到：
[0109] 在每一次迭代过程中都检测相位约束条件，若则证 明为相位畸变，利用估计窗口中的最大概率相位差替代f +270进行相 位替换，保证两次搜索的估计窗口相位差统计值最小，否则按照步骤7中的方法自动改变窗 口的大小和方向。
[0110] 步骤11、当判断窗口搜索完整个干涉图则解缠结束，得到解缠结果。
[0111]以上借助具体实施例对本发明做了进一步描述，但是应该理解的是，这里具体的 描述，不应理解为对本发明的实质和范围的限定，本领域内的普通技术人员在阅读本说明 书后对上述实施例做出的各种修改，都属于本发明所保护的范围。
【主权项】
1. 一种基于相位差统计模型的InSAR图像相位解缠方法，其特征在于，包括w下步骤： 步骤一:假定待处理干设图像行、列像素大小为MXN，将其划分为kXk大小的局部估计 窗口，计算相位在行列两个方向的差值Δ Φ，并对相位差在方位向和距离向求导，得到相位 相位差导数(离散相位下的相位频率）其中，Δxi,j = Xi+l,j+Xi-l,j-2xi,j，Δyi,j = yi,j+l+yi,j-1-2yi,j分别表示干设图沿x、y方向的相 位差值 分别是kXk窗口内的x、y方向的相位差均值；^是W窗口中屯、的 啤说- 相位差导数标准偏差值； 步骤二:根据上步相位差导数值，计算局部估计窗口内相位差统计模型，并估计局部窗 口内地形高度变化量的分布； 步骤Ξ:根据相位差频率与地形坡度角的关系式估计地形高度变化量； 步骤四：获取局部估计窗口内相位差统计模型； 步骤五:利用估计窗口内地形变化量分布函数对窗口内相位差进行平滑，在估计窗口 内将距离向和方位向相位差值乘W相位差统计模型得到窗口内顾及地形变化量的相位差 估计模型； 步骤六:构造基于相位差统计模型的InSAR图像最小二乘相位解缠模型； 步骤屯:对解缠前后距离向和方位向相位差的整数倍数差值Aky和Akx进行估计并积 分，就可W得到解缠后的相位差整数，将求解相位差的问题转换为求解整数倍数差值的问 题，构造整数倍数参数的函数方程式为：其 中Δ Φ y和Δ Φ X表示相位点缠绕相位差，Af,和Δ從表示解缠后相位差； 步骤八:将最小二乘相位解缠约束模型替换为整数倍数差的约束最小化模型，如公式：其中，P(ΔΦx)、P(ΔΦy)表示相位点行向量、列向量相位差统计值，Λkt与ΔkW表示 两次迭代的整数倍数差，ξ为极小值;约束条件S. t表示在局部窗口内，第t次和第t-1次迭代 求解时的解缠整数倍数差值满足差值最小的要求，K表示估计窗口的大小； 步骤九:利用离散最小二乘迭代方法获取解缠相位。
【文档编号】G01S13/90GK106093939SQ201610362293
【公开日】2016年11月9日
【申请日】2016年5月27日
【发明人】不公告发明人
【申请人】山东科技大学