利用非线性矢量曲面构建InSAR相位图像模型的方法
【专利摘要】本发明涉及雷达、光学干涉图像处理，尤其涉及合成孔径雷达干涉测量(InSAR)相位图像处理。首先对离散InSAR相位图像进行多层次小波分解，将相位表示为低频部分和高频部分和的形式；利用估计窗口内相位点高频分量的最大值确定主矢量的方向，根据信号和残差点的模极值变化不同，对估计窗口内相位进行平滑；在估计窗口内对各分解方向求取相位矢量和，构建非线性相位矢量曲面模型；利用非线性相位矢量曲面顾及测区地形因素起伏变化的非线性部分，得到的相位模型能预估测区相位大小变化趋势，帮助相位滤波和解缠，提高处理结果的精度。
【专利说明】
利用非线性矢量曲面构建I nSAR相位图像模型的方法
技术领域
[0001]本发明涉及雷达、光学干涉图像处理，尤其涉及合成孔径雷达干涉测量（InSAR)相 位图像处理。
【背景技术】
[0002] InSAR技术利用侧视雷达回波信号所携带的相位信息来计算地面的三维数据。从 最初的雷达数据格式SLC(单视复图像)得到InSAR相位图像，再得到最终高程数据产品，需 要经过两幅图像配准、复共辄相乘、相位滤波、相位解缠和基线估计等数据处理步骤才能得 到数字高程模型。
[0003] InSAR相位图像是SAR数据处理的中间产品，也是图像滤波和相位解缠的基础数 据，由于后续的步骤需要转换成数学模型进行处理，构建精确、有效的InSAR图像的相位模 型是重点要解决的技术难题之一。
[0004] InSAR干涉图像可表示为二维离散相位信号形式，干涉相位Φ可写成线性项 Φ linear和非线性项Φ_-linear和的形式，如公式（1-1):
(1.-1)
[0006] 其中，Dlinear、D_-linear分别表示数据处理系统线性和非线性组成部分;v表示线性 地表起伏变化的速率，A h表示地表起伏的高度变化，T表示获取干涉对的雷达数据时间基 线。公式（1-1)中等号右边的第一项主要是垂直基线、时间基线构成的线性关系；第二项主 要是由数据处理系统和地形欠采样引起相位的非线性变化以及相位残余构成，非线性部分 造成传统枝切线回路无法闭合以及滤波和解缠步骤误差传递等现象，直接用线性的方法无 法获取高精度相位，因此利用相位模型必须要顾及非线性变化的相位，即构建非线性相位 模型。
[0007] 目前，InSAR相位模型可通过线性相位估值的方法获取，设干涉图的线性相位模型 为U，那么：
(1-2)
[0009]其中，尤，/;为通过二维Chirp-Z变换得到的二维空间频率分量，X?对应于距离 向频谱偏移量Α ω
?其中Fs为系统的距离向采样率，C(x，y)为被截取局部干涉图。 详见文南犬：Bamler R，Adaln N，et al .Noise-Induced Slope Distortion in 2_D Phase Unwrapping by Linear Estimators with Application to SAR Interferometry，1998， 中公开的技术方案。从（1-2)式可以看出线性相位模型可以用一阶近似相位平面和一个复 常数部分构成。算法求解过程中采用二维离散傅里叶变换(FFT)方法，通过搜索窗口内二维 信号频谱的峰值位置来估计干涉条纹频率，取估计窗口的中心相位频率作为信号的相位频 率估计，通过最大化估计窗口相位频率的似然函数，可以估计出窗口中心(m, η)处的真实局 部空间频率fx和fy。经过二维离散傅立叶变换(DFT)后，得到频谱峰值处的频率能够近似局 部频率，详见文献:朱岱寅.利用线性相位模型提高干涉SAR图像相干性及改进相干性估计， 2005,中公开的技术方案。
[0010]这是相位模型的线性构造方法，模型假设估计窗口内二维相位满足线性变化的特 征，然而在实际数据处理中，线性相位模型仍然存在以下缺陷：
[0011] (1)相位残差点反映在频率域中是远离估计频率的干扰信号，此时使用估计窗口 的方法估计的相位不再满足线性模型假设，容易导致线性相位估计模型出现偏差。
[0012] (2)线性模型平滑了相位中的非线性变化部分，使用频谱峰值的最大化似然函数 使得窗口内相位样本局部统一，模糊了相位变化的细节。
[0013] (3)线性模型的特点导致估计窗口小且固定，计算时间长，不能有效抑制残差点导 致低相干区域出现较大的模型误差。
[0014] 目前，对于非线性相位使用二维离散傅里叶变换(FFT)方法，相位估计窗口可以根 据均值的大小自适应变化，反映了部分地形变化的轮廓，详见文献:Wang Q S，Huang H F， An X Y,et al. Improving phase Unwrapping Techniques by the Use of Nonlinear Phase Model，2011，中公开的技术方案。但FFT方法仍不能顾及地形和相位变化的细节，无 法依据地形变化的视线向、距离向和方位向进行区分，而且使用估计窗口的均值估计相位 不准确，也容易造成误差的传递。
【发明内容】
[0015] 本发明对离散InSAR相位图像进行多层次小波分解，将相位表示为低频部分和高 频部分和的形式;利用估计窗口内相位点高频分量的最大值确定主矢量的方向，根据信号 和残差点的模极值变化不同，对估计窗口内相位进行平滑;在估计窗口内对各分解方向求 取相位矢量和，构建非线性相位矢量曲面模型；利用非线性相位矢量曲面顾及测区地形因 素起伏变化的非线性部分，得到的相位模型能预估测区相位大小变化趋势，帮助相位滤波 和解缠，提高处理结果的精度。
[0016]具体方案流程为：
[0017] 1.1将相位干涉图以复数的浮点形式读取到Matlab中，标记矩阵行、列大小为M、N， 单位为像素。将干涉相位图划分大小为kXl的初始估计窗口，并计算值得信任的第一个估 计窗口，将窗口内相位点(m，n)单频信号形式转换为三角函数弦分量的复数相位形式：Φ (m，n)=angle(exp[ j(fxm+fyn)]) = cos(fxm+fyn)+j sin(fxm+fyn)，其中，fx和fy分别为相位 点行向量和列向量的频率；
[0018] 1.2在估计窗口内，利用离散小波变换对干涉图像正弦分量和余弦分量进行多层 分解，将相位表示为低频部分和高频部分和的形式；
[0019] 1.3在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一 步搜索窗口的方向；
[0020] 1.4在估计窗口内，计算相位主矢量，利用高频部分主矢量平滑窗口内的相位频 率；
[0021 ] 1.5在估计窗口内，根据各层分解分量的高频矢量和构建顾及相位细节变化的非 线性相位曲面模型，将相位表示成非线性相位和线性相位和的形式Φ = Φ line3ar + J _ _ _ Φ ii_r，对非线性部分利用(爪，》) = arg(Z< )= m，g(l + 6十表;）对j层分解的三 Η 个高频方向求矢量和，其中j = 1，2···J，^.、分别表示对j层高频的水平、垂直和视 线方向分量求矢量和，矢量求和保证了非线性曲面的特性；
[0022] 1.6对小波分解相位进行逆变换得到相位估计模型，在估计窗口 k X 1内，对相位点 也,1 (m，η)，利用(电，《) = a取(紅_1 = arg (I + ΣI j可得到估计窗口内非线 性相位估计模型，其中，WTl ·)为小波逆变换，(^为小波分解低频部分，ΣΙ为小波分解 高频部分和。
[0023] 本发明技术方案带来的有益效果
[0024] (1)利用构造的非线性相位模型顾及了非线性相位的偏差，实现了相位最优估计；
[0025] (2)相位滤波是去除干涉图中的噪声，而应区分地形因素引起的相位相干斑噪声； 相位解缠是将缠绕相位（以2π为模)展开成实际相位的过程，但由于受到地形因素的影响使 得解缠路径有环路，不能正确解缠，若使用本发明估计相位模型来顾及测区地形因素，能预 估测区相位变化趋势，帮助相位滤波和解缠，提高处理结果的精度。
[0026] (3)非线性相位模型不再以局部线性相位频率估计代替局部相位，而是以反映地 形因素的非线性曲面来近似估计窗口相位。
【附图说明】
[0027] 图1为以相位垂直方向为主方向的相位矢量求和图；
[0028] 图2为以相位水平方向为主方向的相位矢量求和图；
[0029] 图3为以相位视线方向为主方向的相位矢量求和图；
[0030] 图4为非线性矢量曲面构建InSAR相位图像流程图。
[0031 ]图例说明，1-侧视雷达视线向，2-水平方向，3-垂直方向，V -侧视雷达视线向反方 向，2'-水平方向反方向，3'-垂直方向反方向。
【具体实施方式】
[0032]下面结合附图对本发明技术方案进一步说明。
[0033]由于星载侧视雷达以相同频率发射窄带脉冲，在斜距相同的条件下，不同地形条 件反射的雷达信号不同，反映在干涉图中就是或密集或稀疏的干涉条纹，那么相位条纹的 疏密程度也反映了干涉图中地形的特征，可以根据相位条纹频率建立相位模型。
[0034]为了方便本发明中技术描述，首先进行了以下定义：
[0035] 定义1:采样点地形坡度的相位差导数表示。
[0036] 设参数τ为地面采样点的地形坡度，τ与相位信号系统各参数之间的关系为：
P)
[0038]其中，c为雷达波传播速度;B丄为基线在垂直于视线方向的投影，即有效基线长度； λ为波长;R为斜距；Θ为入射角;?·φ为干涉图相位差频率。根据公式(2)可推出地面倾角与入 射角之间的关系公式(3):
(3)
[0040]由上式可知，可以根据距离向和方位向上的局部频率Af来估计干涉图相位模型。 设Δ fx、△ fy分别为距离向、方位向相位差频率，可得下式(4):
[0042] 联立相位差的公式： 可得：
u
[0046]本方法利用离散小波分层分解，获取非线性相位曲面模型来逼近估计窗口内非线 性相位，并对曲面模型进行逆变换得到非线性相位模型。具体步骤如下：
[0047]步骤1:准备基础数据的工作。
[0048] 获取两幅SAR图像51与52配准、滤波后的干涉图像，包括附加的系统参数信息（波 长、基线、下视角、斜距和如果可以获取的地形DEM信息），计算干涉图像相干系数γ和离散 干涉相位Φ。在实际图像处理过程中，
，式中1表示多视 数，取值为1:5或2:10，表示SAR图像共辄相乘。
[0049] 步骤2:读取干涉图得到相位矩阵并转换成相位频率的形式。
[0050] 将相位干涉图以复数的浮点形式读取到Matlab中，设矩阵行、列大小为Μ、Ν，单位 为像素。将干涉相位图划分大小为ΚXL的初始估计窗口，并将窗口内相位点(m，η)的单频信 号形式转换为三角函数弦分量的复数相位形式： i,kw,n、=-"ngle(CKp j{ f'iii + f'n] \
[0051 ] / 1 L , J/ ?? =cos ( /> + /;") + / sin (./:/" + />)
[0052] 其中，f4Pfy分别为相位点行向量和列向量的频率。
[0053] 另外，使用相位导数也就是相位差的形式得到相位频率，用下面的公式进行计算，
[0054] L = 〇Φ/^ = φη1ιυι + ~2φηιη
[0055] f, = 9φ/δν = 4,"+ι + ~ 2Φη,."
[0056] 上式中φω,η表示干涉图像ΜΧΝ中相位点(m，n)的相位值;χ和y为矩阵行向量和列 向量。
[0059] 在实际离散相位图像算法中，根据第1步的相位相干系数搜索干涉图像，计算值得 信任的第一个估计窗口，可使用下面的方法:首先搜索干涉图中局部估计窗口的平均最大 相干系数若这样的窗口数量不止一个，设数量为N，则根据地形坡度角与相位频率 V kxk J 公式(3 )，求取窗口内最小地形坡度窗口设数量为Μ，确定为初始估计窗口；否则，若f > 1， 则将窗口下标值最小的窗口确定为第一个估计窗口，从第一个估计窗口开始进行相位模型 的计算。
[0060] 步骤3:利用离散小波变换的方法对相位频率进行分解。
[0061]在估计窗口内，利用离散小波变换对干涉图像正弦分量和余弦分量进行多层分 解，将相位表示为低频部分和高频部分和的形式。对估计相位Φ在矢量域求和，分别用不同 频率的矢量形式表示(根据频率变化的特点，使用小波变换后的不同分辨率表示，低频部分 表示非线性相位;高频部分为相位的细节信息，地形的起伏变化反映在频率上就是频率域 中的高频部分）。
[0062]相位Φ k, 1 (m，η)的小波频谱变换为：
[0064] 其中，Φ 〇为窗口中心相位点，成=训沙(為(L/)exp(-./2；r(./;/c +//))); Δ φ邊示 Φ 〇的微小变化区域;J为相位分解的层次J = 1，2，…，n;Cj和Dh，l，r分别表亦小波分解中的低 频和高频部分，其中Dh，l，r表示各个分辨率层的相位细节信息，H，L，R表示窗口相位搜索的三 个方向（水平、垂直、视线)或其反方向，详细的相位搜索方向可参见附图1、2。
[0065]在实际离散相位图像算法中，对于相位窗口的第t次分解可以由第t-Ι次分解得 到，t-Ι次的分解窗口结构为：
[0067] X表示估计窗口中的相位频率。
[0068] 步骤4:利用分层分解结果获取高频部分主矢量的方向，确定估计窗口滑动的方 向。
[0069] 根据公式(8)对估计窗口内小波分解的结果，得到窗口内相位点的频率表示形式：
[0071] 其中，Φ?ρ为小波分解低频部分K为小波分解高频部分;Pi为取值范围为(〇，1) 的加权系数;在矢量加权求和中，系数Pi的大小决定了第i个信号分量在估计窗口中的权 重。
[0072] 在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一步 搜索窗口的方向。设当前窗口 j层的主矢量方向为M，标记为^^，是根据j-Ι层的主矢量方向 确定的：
[0075] 上式中表示H、L、R三个方向频率对主矢量的角度偏差，也就是三分量在 主矢量方向的矢量分解。
[0076] 在实际离散相位图像算法中，对矢量分解的角度偏差通过如下公式进行计算：
[0080]在上式中，若视线向为主矢量方向，水平、垂直方向矢量的角度偏差用A 01/31获 取;若水平、垂直为主矢量方向，视线向矢量的角度偏差用Α Θ2Λτ、△ θ3/3?获取。
[0081 ]步骤5:在估计窗口内利用高频部分主矢量平滑窗口内的相位频率。
[0082]由于窗口内的相位噪声表现为信号无关性，若当前估计窗口判断为相位噪声的相 位点频率值大于窗口内相位频率值的和，则认定为噪声信号，其非线性相位为无效估计，利 用已经构建相位模型的相邻窗口频率均值，与当前窗口高频部分的数值之和作为当前窗口 .A'xL 的相位模型，若物?，《)>，贝1J: ?=1
[0084] 根据干涉相位频率提取相位矢量，可以设置阈值(窗口频率变换主瓣值的均值）， 当估计窗口内提取得到的频率值大于该门限时，由于窗口中频谱变化最大值表示了相位变 化的主方向，发明中使用窗口内信号频谱最大值作为矢量阈值：Ize
[0085] 对于高频部分的处理为：= Σ U.A. +/d- + /"A) 'H，L，R表示三个方向（水 平、垂直、视线）的正方向高频gz和低频hz，在方向搜索可时同时顾及这三个方向的反方向。 若此时根据要控制的方向个数(2个、4个、6个或8个)对各方向求矢量和。
[0086] 在实际离散相位图像算法中，相位高频部分通过如下公式进行计算：
[0087] Σ (尽--'?--十 W)= Σ (2_, {^ !D,, (in,n)-D, \ym-[k-l)j2,n-(k-l)l2^Cj (>?,?)] j = zeH,L,R z.eH,L,R + (2-7 (2-人1?_ (辦，π)-為(m -(左-1)/2，'卜(免" 1)/2)) (67)(肌？7)) '2~i(2-iDz(nKn)-Dz(m-(k~l)/2^i-(k-l)/2)f Η- |2-?λ (?；, η)- (m-(k-1)/2 .η-(k-I )/2))
[0088] 上式中(m,n)为估计窗 口中相位点，02(111-(1^-1)/2，11-(1^-1)/2)表示(111，11)邻域的 相位点。
[0089] (k-Ι )/2保证相位变换点在估计窗口范围内。
[0090] 步骤6:在估计窗口内利用高频部分矢量和构造非线性相位曲面模型。
[0091] 由公式（2) Φ = Φ η_γ+Φ non-linear ? 根据各层分解分量高频矢量和构造顾及相位 细节变化的构建非线性相位曲面模型。
[0092] 设式（8)中
为信号Φ:分解的i层主矢量， 那么主矢量和的相位反映非线性相位巾_得到公式(10): J
[0093] 為―·(衍，= ) (1〇)
[0094] 同时对局部窗口相位高频矢量求和保证了非线性曲面的特性，如公式(11):
[0096] 其中，上式中右边g、慰、g表示对j层高频的水平、垂直和视线方向分量求矢 量和。砹1表示视线向为主矢量方向，水平、垂直方向矢量和用(/? 1获取;若水 平、垂直为主矢量方向，视线向矢量的矢量和分别用、GfA爲/;τ获取。
[0097] 步骤7:对小波分解相位进行逆变换，可求取非线性相位估计模型。
[0098] 在估计窗口（k，l)内可以得到非线性相位： = arg^T'1 {fm + ./>?})
[0099] f r Λ =ill'S Σ+ AnM," {m-fl)i-'xP(-/2^(./> + A11)) Π2) v k,i ' ' J
[0100] 其中，胃!^1!： ·)为小波逆变换，<Kp为小波分解低频部分为小波分解高频部 分在实际离散相位图像算法中，对于估计窗口的第t次逆变换由第t+i次分解得到，t+i次的 逆变换窗口结构为：
[0102] X-1表示Cj或Dh,l,r的逆变换。
[0103] 步骤8:判断是否搜索完整个图像，是结束流程;否则，窗口向相位极值方向滑动。
[0104] 以上借助具体实施例对本发明做了进一步描述，但是应该理解的是，这里具体的 描述，不应理解为对本发明的实质和范围的限定，本领域内的普通技术人员在阅读本说明 书后对上述实施例做出的各种修改，都属于本发明所保护的范围。
【主权项】
1.利用非线性矢量曲面构建InSAR相位图像模型的方法，其特征在于，包括w下步骤： 步骤一:将相位干设图W复数的浮点形式读取到Matlab中，标记矩阵行、列大小为M、N， 单位为像素;将干设相位图划分大小为kXl的初始估计窗口，并计算值得信任的第一个估 计窗口，将窗口内相位点(m，n)单频信号形式转换为Ξ角函数弦分量的复数相位形式：Φ (m,n) =angle(e邱[j(fxm+fyn)]) =cos(fxm+fyn)+jsin(fxm+fyn)，其中，fx和fy分别为相位 点行向量和列向量的频率； 步骤二:在估计窗口内，利用离散小波变换对干设图像正弦分量和余弦分量进行多层 分解，将相位表示为低频部分和高频部分和的形式； 步骤Ξ:在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一 步捜索窗口的方向； 步骤四：在估计窗口内，计算相位主矢量，利用高频部分主矢量平滑窗口内的相位频 率. 步骤五:在估计窗口内，根据各层分解分量的高频矢量和构建顾及相位细节变化的非 线性相位曲面模型，将相位表示成非线性相位和线性相位和的形式Φ = Φ linear + Φ ηαη-1 inear，对非线性部分利月对j层分解的Ξ 个高频方向求矢量和，其中j = l，2...J，ft、敵、ft分别表示对j层高频的水平、垂直和视 线方向分量求矢量和，矢量求和保证了非线性曲面的特性； 步骤六:对小波分解相位进行逆变换得到相位估计模型，在估计窗口 k X 1内，对相位点 4k,i(m，n)，利用可得到估计窗口内非线 性相位估计模型，其中，wri(.)为小波逆变换，(iHp为小波分解低频部分，Σ拓,,为小波分解 高频部分和。
【文档编号】G06T9/00GK106097404SQ201610363407
【公开日】2016年11月9日
【申请日】2016年5月27日
【发明人】不公告发明人
【申请人】山东科技大学