CN106097404A - 利用非线性矢量曲面构建InSAR相位图像模型的方法 - Google Patents

Abstract

本发明涉及雷达、光学干涉图像处理，尤其涉及合成孔径雷达干涉测量(InSAR)相位图像处理。首先对离散InSAR相位图像进行多层次小波分解，将相位表示为低频部分和高频部分和的形式；利用估计窗口内相位点高频分量的最大值确定主矢量的方向，根据信号和残差点的模极值变化不同，对估计窗口内相位进行平滑；在估计窗口内对各分解方向求取相位矢量和，构建非线性相位矢量曲面模型；利用非线性相位矢量曲面顾及测区地形因素起伏变化的非线性部分，得到的相位模型能预估测区相位大小变化趋势，帮助相位滤波和解缠，提高处理结果的精度。

Description

利用非线性矢量曲面构建InSAR相位图像模型的方法

技术领域

本发明涉及雷达、光学干涉图像处理，尤其涉及合成孔径雷达干涉测量(InSAR)相位图像处理。

背景技术

InSAR技术利用侧视雷达回波信号所携带的相位信息来计算地面的三维数据。从最初的雷达数据格式SLC(单视复图像)得到InSAR相位图像，再得到最终高程数据产品，需要经过两幅图像配准、复共轭相乘、相位滤波、相位解缠和基线估计等数据处理步骤才能得到数字高程模型。

InSAR相位图像是SAR数据处理的中间产品，也是图像滤波和相位解缠的基础数据，由于后续的步骤需要转换成数学模型进行处理，构建精确、有效的InSAR图像的相位模型是重点要解决的技术难题之一。

InSAR干涉图像可表示为二维离散相位信号形式，干涉相位Φ可写成线性项φ_linear和非线性项φ_non-linear和的形式，如公式(1-1)：

$\begin{array}{c}\mathrm{\Φ}={\mathrm{\φ}}_{linear}+{\mathrm{\φ}}_{non-linear}=\frac{4\mathrm{\π}}{\mathrm{\λ}}\left(\frac{B_{\⊥}}{R\sin \mathrm{\θ}}\mathrm{\Δ}h+v\·T\right)+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·D_{non-linear}\\ =\frac{4\mathrm{\π}}{\mathrm{\λ}}\left(D_{linear}+D_{non-linear}\right)\end{array}---\left(1-1\right)$

其中，D_linear、D_non-linear分别表示数据处理系统线性和非线性组成部分；v表示线性地表起伏变化的速率，Δh表示地表起伏的高度变化，T表示获取干涉对的雷达数据时间基线。公式(1-1)中等号右边的第一项主要是垂直基线、时间基线构成的线性关系；第二项主要是由数据处理系统和地形欠采样引起相位的非线性变化以及相位残余构成，非线性部分造成传统枝切线回路无法闭合以及滤波和解缠步骤误差传递等现象，直接用线性的方法无法获取高精度相位，因此利用相位模型必须要顾及非线性变化的相位，即构建非线性相位模型。

目前，InSAR相位模型可通过线性相位估值的方法获取，设干涉图的线性相位模型为l_φ，那么：

$l_{\mathrm{\φ}}(x,y)=\frac{\underset{x,y\∈w_s}{\Σ}C(x,y)e^{-j2\mathrm{\π}(f_x^l+f_y^l)}}{\left|\underset{x,y\∈w_s}{\Σ}C(x,y)e^{-j2\mathrm{\π}(f_x^l+f_y^l)}\right|}\·e^{j2\mathrm{\π}(f_x^l+f_y^l)}---(1-2)$

其中，为通过二维Chirp-Z变换得到的二维空间频率分量，对应于距离向频谱偏移量Δω，其中F_s为系统的距离向采样率，C(x,y)为被截取局部干涉图。详见文献：Bamler R,Adaln N,et al.Noise-Induced Slope Distortion in 2-D PhaseUnwrapping by Linear Estimators with Application to SAR Interferometry，1998，中公开的技术方案。从(1-2)式可以看出线性相位模型可以用一阶近似相位平面和一个复常数部分构成。算法求解过程中采用二维离散傅里叶变换(FFT)方法，通过搜索窗口内二维信号频谱的峰值位置来估计干涉条纹频率，取估计窗口的中心相位频率作为信号的相位频率估计，通过最大化估计窗口相位频率的似然函数，可以估计出窗口中心(m,n)处的真实局部空间频率f_x和f_y。经过二维离散傅立叶变换(DFT)后，得到频谱峰值处的频率能够近似局部频率，详见文献：朱岱寅.利用线性相位模型提高干涉SAR图像相干性及改进相干性估计，2005，中公开的技术方案。

这是相位模型的线性构造方法，模型假设估计窗口内二维相位满足线性变化的特征，然而在实际数据处理中，线性相位模型仍然存在以下缺陷：

(1)相位残差点反映在频率域中是远离估计频率的干扰信号，此时使用估计窗口的方法估计的相位不再满足线性模型假设，容易导致线性相位估计模型出现偏差。

(2)线性模型平滑了相位中的非线性变化部分，使用频谱峰值的最大化似然函数使得窗口内相位样本局部统一，模糊了相位变化的细节。

(3)线性模型的特点导致估计窗口小且固定，计算时间长，不能有效抑制残差点导致低相干区域出现较大的模型误差。

目前，对于非线性相位使用二维离散傅里叶变换(FFT)方法，相位估计窗口可以根据均值的大小自适应变化，反映了部分地形变化的轮廓，详见文献：Wang Q S,Huang H F,An X Y,et al.Improving phase Unwrapping Techniques by the Use of NonlinearPhase Model，2011，中公开的技术方案。但FFT方法仍不能顾及地形和相位变化的细节，无法依据地形变化的视线向、距离向和方位向进行区分，而且使用估计窗口的均值估计相位不准确，也容易造成误差的传递。

发明内容

本发明对离散InSAR相位图像进行多层次小波分解，将相位表示为低频部分和高频部分和的形式；利用估计窗口内相位点高频分量的最大值确定主矢量的方向，根据信号和残差点的模极值变化不同，对估计窗口内相位进行平滑；在估计窗口内对各分解方向求取相位矢量和，构建非线性相位矢量曲面模型；利用非线性相位矢量曲面顾及测区地形因素起伏变化的非线性部分，得到的相位模型能预估测区相位大小变化趋势，帮助相位滤波和解缠，提高处理结果的精度。

具体方案流程为：

1.1将相位干涉图以复数的浮点形式读取到Matlab中，标记矩阵行、列大小为M、N，单位为像素。将干涉相位图划分大小为k×l的初始估计窗口，并计算值得信任的第一个估计窗口，将窗口内相位点(m,n)单频信号形式转换为三角函数弦分量的复数相位形式：φ(m,n)＝angle(exp[j(f_xm+f_yn)])＝cos(f_xm+f_yn)+j sin(f_xm+f_yn)，其中，f_x和f_y分别为相位点行向量和列向量的频率；

1.2在估计窗口内，利用离散小波变换对干涉图像正弦分量和余弦分量进行多层分解，将相位表示为低频部分和高频部分和的形式；

1.3在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一步搜索窗口的方向；

1.4在估计窗口内，计算相位主矢量，利用高频部分主矢量平滑窗口内的相位频率；

1.5在估计窗口内，根据各层分解分量的高频矢量和构建顾及相位细节变化的非线性相位曲面模型，将相位表示成非线性相位和线性相位和的形式Φ＝φ_linear+φ_non-linear，对非线性部分利用对j层分解的三个高频方向求矢量和，其中j＝1,2…J，分别表示对j层高频的水平、垂直和视线方向分量求矢量和，矢量求和保证了非线性曲面的特性；

1.6对小波分解相位进行逆变换得到相位估计模型，在估计窗口k×l内，对相位点ψ_k,l(m,n)，利用可得到估计窗口内非线性相位估计模型，其中，WT^-1(·)为小波逆变换，φ_Lp为小波分解低频部分，为小波分解高频部分和。

本发明技术方案带来的有益效果

(1)利用构造的非线性相位模型顾及了非线性相位的偏差，实现了相位最优估计；

(2)相位滤波是去除干涉图中的噪声，而应区分地形因素引起的相位相干斑噪声；相位解缠是将缠绕相位(以2π为模)展开成实际相位的过程，但由于受到地形因素的影响使得解缠路径有环路，不能正确解缠，若使用本发明估计相位模型来顾及测区地形因素，能预估测区相位变化趋势，帮助相位滤波和解缠，提高处理结果的精度。

(3)非线性相位模型不再以局部线性相位频率估计代替局部相位，而是以反映地形因素的非线性曲面来近似估计窗口相位。

附图说明

图1为以相位垂直方向为主方向的相位矢量求和图；

图2为以相位水平方向为主方向的相位矢量求和图；

图3为以相位视线方向为主方向的相位矢量求和图；

图4为非线性矢量曲面构建InSAR相位图像流程图。

图例说明，1-侧视雷达视线向，2-水平方向，3-垂直方向，1′-侧视雷达视线向反方向，2′-水平方向反方向，3′-垂直方向反方向。

具体实施方式

下面结合附图对本发明技术方案进一步说明。

由于星载侧视雷达以相同频率发射窄带脉冲，在斜距相同的条件下，不同地形条件反射的雷达信号不同，反映在干涉图中就是或密集或稀疏的干涉条纹，那么相位条纹的疏密程度也反映了干涉图中地形的特征，可以根据相位条纹频率建立相位模型。

为了方便本发明中技术描述，首先进行了以下定义：

定义1:采样点地形坡度的相位差导数表示。

设参数τ为地面采样点的地形坡度，τ与相位信号系统各参数之间的关系为：

$f\left(\mathrm{\φ}\right)=\frac{{\mathrm{cB}}_{\⊥}}{\mathrm{\λ}Rtan(\mathrm{\θ}-\mathrm{\τ})}---\left(2\right)$

其中，c为雷达波传播速度；B_⊥为基线在垂直于视线方向的投影，即有效基线长度；λ为波长；R为斜距；θ为入射角；f_φ为干涉图相位差频率。根据公式(2)可推出地面倾角与入射角之间的关系公式(3)：

$\mathrm{\θ}-\mathrm{\τ}=arctan\[\frac{{\mathrm{cB}}_{\⊥}}{{\mathrm{\λRf}}_{\mathrm{\φ}}}\]---\left(3\right)$

由上式可知，可以根据距离向和方位向上的局部频率Δf来估计干涉图相位模型。设Δf_x、Δf_y分别为距离向、方位向相位差频率，可得下式(4)：

$\mathrm{\Δ}f={\[{\mathrm{\Δf}}_x{\mathrm{\Δf}}_y\]}^T=\frac{1}{2\mathrm{\π}}{\[\frac{\∂\mathrm{\φ}}{\∂x}\frac{\∂\mathrm{\φ}}{\∂y}\]}^T=\frac{1}{2\mathrm{\π}}{\[\frac{\mathrm{\Δ}\mathrm{\φ}}{\mathrm{\Δ}x}\·\frac{\mathrm{\Δ}\mathrm{\φ}}{\mathrm{\Δ}y}\]}^T---\left(4\right)$

联立相位差的公式：可得：

$\begin{array}{c}{\mathrm{\Δf}}_y=\frac{\∂\mathrm{\φ}}{\∂y}=k_B\frac{1}{tan(\mathrm{\θ}-{\mathrm{\τ}}_y)}\\ =k_B\frac{1+{\mathrm{tan\θtan\τ}}_y}{\tan \mathrm{\θ}-{\mathrm{tan\τ}}_y}\end{array}---\left(5\right)$

$\begin{array}{c}{\mathrm{\Δf}}_x=\frac{\∂\mathrm{\φ}}{\∂x}=k_B\frac{{\mathrm{tan\τ}}_x{\mathrm{cos\τ}}_y}{sin(\mathrm{\θ}-{\mathrm{\τ}}_y)}\\ =k_B\frac{{\mathrm{tan\τ}}_x}{\sin \mathrm{\θ}-{\mathrm{cos\θtan\τ}}_y}\end{array}---\left(6\right)$

其中：

本方法利用离散小波分层分解，获取非线性相位曲面模型来逼近估计窗口内非线性相位，并对曲面模型进行逆变换得到非线性相位模型。具体步骤如下：

步骤1：准备基础数据的工作。

获取两幅SAR图像S₁与S₂配准、滤波后的干涉图像，包括附加的系统参数信息(波长、基线、下视角、斜距和如果可以获取的地形DEM信息)，计算干涉图像相干系数γ和离散干涉相位φ。在实际图像处理过程中，式中l表示多视数，取值为1:5或2:10，表示SAR图像共轭相乘。

步骤2：读取干涉图得到相位矩阵并转换成相位频率的形式。

将相位干涉图以复数的浮点形式读取到Matlab中，设矩阵行、列大小为M、N，单位为像素。将干涉相位图划分大小为K×L的初始估计窗口，并将窗口内相位点(m,n)的单频信号形式转换为三角函数弦分量的复数相位形式：

$\begin{array}{c}\mathrm{\φ}\left(m,n\right)=angle\left(\exp \[j\left(f_{x}m+f_{y}n\right)\]\right)\\ =\cos \left(f_{x}m+f_{y}n\right)+j\sin \left(f_{x}m+f_{y}n\right)\end{array}---\left(7\right)$

其中，f_x和f_y分别为相位点行向量和列向量的频率。

另外，使用相位导数也就是相位差的形式得到相位频率，用下面的公式进行计算，

$f_x=\∂\mathrm{\φ}/\∂x={\mathrm{\φ}}_{m+1,n}+{\mathrm{\φ}}_{m-1,n}-2{\mathrm{\φ}}_{m,n}$

$f_y=\∂\mathrm{\φ}/\∂y={\mathrm{\φ}}_{m,n+1}+{\mathrm{\φ}}_{m,n-1}-2{\mathrm{\φ}}_{m,n}$

上式中φ_m,n表示干涉图像M×N中相位点(m,n)的相位值；x和y为矩阵行向量和列向量。

$cos(f_{x}m+f_{y}n)=\frac{\underset{m=1}{\overset{M-1}{\Σ}}\underset{n=1}{\overset{N}{\Σ}}({\mathrm{\φ}}_{m+1,n}+{\mathrm{\φ}}_{m-1,n}-2{\mathrm{\φ}}_{m,n})}{{\left(\underset{m=1}{\overset{M-1}{\Σ}}\underset{n=1}{\overset{N-1}{\Σ}}\right({({\mathrm{\φ}}_{m+1,n}+{\mathrm{\φ}}_{m-1,n}-2{\mathrm{\φ}}_{m,n})}^2+{({\mathrm{\φ}}_{m,n+1}+{\mathrm{\φ}}_{m,n-1}-2{\mathrm{\φ}}_{m,n})}^2\left)\right)}^{1/2}}$

$sin(f_{x}m+f_{y}n)=\frac{\underset{m=1}{\overset{M}{\Σ}}\underset{n=1}{\overset{N-1}{\Σ}}({\mathrm{\φ}}_{m,n+1}+{\mathrm{\φ}}_{m,n-1}-2{\mathrm{\φ}}_{m,n})}{{\left(\underset{m=1}{\overset{M-1}{\Σ}}\underset{n=1}{\overset{N-1}{\Σ}}\right({({\mathrm{\φ}}_{m+1,n}+{\mathrm{\φ}}_{m-1,n}-2{\mathrm{\φ}}_{m,n})}^2+{({\mathrm{\φ}}_{m,n+1}+{\mathrm{\φ}}_{m,n-1}-2{\mathrm{\φ}}_{m,n})}^2\left)\right)}^{1/2}}$

在实际离散相位图像算法中，根据第1步的相位相干系数搜索干涉图像，计算值得信任的第一个估计窗口，可使用下面的方法：首先搜索干涉图中局部估计窗口的平均最大相干系数若这样的窗口数量不止一个，设数量为N，则根据地形坡度角与相位频率公式(3)，求取窗口内最小地形坡度窗口设数量为N′，确定为初始估计窗口；否则，若N′>1，则将窗口下标值最小的窗口确定为第一个估计窗口，从第一个估计窗口开始进行相位模型的计算。

步骤3：利用离散小波变换的方法对相位频率进行分解。

在估计窗口内，利用离散小波变换对干涉图像正弦分量和余弦分量进行多层分解，将相位表示为低频部分和高频部分和的形式。对估计相位φ在矢量域求和，分别用不同频率的矢量形式表示(根据频率变化的特点，使用小波变换后的不同分辨率表示，低频部分表示非线性相位；高频部分为相位的细节信息，地形的起伏变化反映在频率上就是频率域中的高频部分)。

相位φ_k,l(m,n)的小波频谱变换为：

$\begin{array}{cc}WT\left({\mathrm{\φ}}_{k,l}\left(m,n\right)\right)=\underset{j}{\Σ}\left(C_j{J}^{j/J}\mathrm{\ψ}\left(\frac{{\mathrm{\φ}}_0-{\mathrm{\Δ\φ}}_t}{J}\right)+\underset{j}{\Σ}{D}_{H,L,R}{J}^{j/J}\mathrm{\ψ}\left(\frac{{\mathrm{\φ}}_0-{\mathrm{\Δ\φ}}_t}{J}\right)\right)& j=1,2,\mathrm{...}J\end{array}---\left(8\right)$

其中，φ₀为窗口中心相位点，Δφ_t表示φ₀的微小变化区域；J为相位分解的层次J＝1,2,…,n；C_J和D_H,L,R分别表示小波分解中的低频和高频部分，其中D_H,L,R表示各个分辨率层的相位细节信息，H,L,R表示窗口相位搜索的三个方向(水平、垂直、视线)或其反方向，详细的相位搜索方向可参见附图1、2。

在实际离散相位图像算法中，对于相位窗口的第t次分解可以由第t-1次分解得到，t-1次的分解窗口结构为：

${\left[\begin{array}{ccccc}{X}_{t+1}\left(0\right)& X_{t-1}\left(1\right)& \mathrm{...}& X_{t-1}\left(1\right)& 0\\ 0& X_{t-1}\left(0\right)& X_{t-1}(p-1)& X_{t-1}(p-2)& \mathrm{...}\\ 0& 0& X_{t-1}\left(0\right)& \mathrm{...}& X_{t-1}\left(p-1\right)\\ X_{t-1}(p-1)& \mathrm{...}& 0& 0& X_{t-1}\left(0\right)\end{array}\right]}_{0<p\&le;\frac{K-1}{2}\&times;\frac{L-1}{2}}$

X表示估计窗口中的相位频率。

步骤4：利用分层分解结果获取高频部分主矢量的方向，确定估计窗口滑动的方向。

根据公式(8)对估计窗口内小波分解的结果，得到窗口内相位点的频率表示形式：

$\begin{array}{c}\underset{i=1}{\overset{K\&times;L}{\&Sigma;}}\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_i^{\&prime;}}\left(m,n\right)=\underset{i=1}{\overset{K\&times;L}{\&Sigma;}}\arg \left({\mathrm{\&phi;}}_{Lp}+\&Sigma;{\mathrm{\&phi;}}_{G_M}\right)\\ \begin{array}{cc}=\underset{i=1}{\overset{K\&times;L}{\&Sigma;}}{P}_i\left(\underset{j}{\&Sigma;}\&lsqb;C_j{J}^{j/J}\mathrm{\&psi;}\left(\frac{{\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t}{J}\right)+\underset{j}{\&Sigma;}{D}_{H,L,R}{J}^{j/J}\mathrm{\&psi;}\left(\frac{{\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t}{J}\right)\&rsqb;\right)& j=0,1,\mathrm{2...}J\end{array}\end{array}---\left(9\right)$

其中，φ_Lp为小波分解低频部分，为小波分解高频部分；P_i为取值范围为(0,1)的加权系数；在矢量加权求和中，系数P_i的大小决定了第i个信号分量在估计窗口中的权重。

在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一步搜索窗口的方向。设当前窗口j层的主矢量方向为M，标记为是根据j-1层的主矢量方向确定的：

$G_M^{j-1}=MAX\left(D_{H,L,R}{J}^{i/J}\mathrm{\&psi;}\left(\frac{{\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t}{J}\right)\right)$

$G_M^j={\mathrm{\&Sigma;D}}_{H,L,R}{J}^{i/J}\mathrm{\&psi;}\left(\frac{{\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t}{J}\right)(G_M^{j-1}\mathrm{\&Delta;}\mathrm{\&theta;}/\mathrm{\&pi;})$

上式中表示H、L、R三个方向频率对主矢量的角度偏差，也就是三分量在主矢量方向的矢量分解。

在实际离散相位图像算法中，对矢量分解的角度偏差通过如下公式进行计算：

${\mathrm{\&Delta;\&theta;}}_1/\mathrm{\&pi;}=\left(1+\frac{{\mathrm{\&pi;G}}_R^{j-1}\cos \left(\mathrm{\&theta;}\right)}{G_H^{j-1}}+\frac{{\mathrm{\&pi;G}}_R^{j-1}\cos \left(\mathrm{\&theta;}\right)}{G_L^{j-1}}\right)\cos \left(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&theta;}\right)$

${\mathrm{\&Delta;\&theta;}}_2/\mathrm{\&pi;}=cos\left(arctan\left({\mathrm{\&pi;G}}_R^{j-1}cos\left(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&theta;}\right)/G_H^{j-1}\right)\right)$

${\mathrm{\&Delta;\&theta;}}_3/\mathrm{\&pi;}=cos\left(arctan\left({\mathrm{\&pi;G}}_R^{j-1}cos\left(\frac{\mathrm{\&pi;}}{2}-\mathrm{\&theta;}\right)/G_L^{j-1}\right)\right)$

在上式中，若视线向为主矢量方向，水平、垂直方向矢量的角度偏差用Δθ₁/π获取；若水平、垂直为主矢量方向，视线向矢量的角度偏差用Δθ₂/π、Δθ₃/π获取。

步骤5：在估计窗口内利用高频部分主矢量平滑窗口内的相位频率。

由于窗口内的相位噪声表现为信号无关性，若当前估计窗口判断为相位噪声的相位点频率值大于窗口内相位频率值的和，则认定为噪声信号，其非线性相位为无效估计，利用已经构建相位模型的相邻窗口频率均值，与当前窗口高频部分的数值之和作为当前窗口的相位模型，若则：

$\begin{array}{c}\begin{array}{cc}\underset{i}{\overset{k\&times;k}{\&Sigma;}}\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_i^{\&prime;}}\left(m,n\right)=\underset{k\&times;k}{\&Sigma;}{P}_M\&lsqb;\left(C_j\mathrm{\&psi;}\left(\left({\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t\right)/j\right)+\underset{i=1}{\overset{M}{\&Sigma;}}{D}_{H,L,R}\mathrm{\&psi;}\left({\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t/J\right)\right)/k\&times;k\&rsqb;& j=1,2,\mathrm{3...}J\end{array}\\ =\underset{k\&times;k}{\&Sigma;}{P}_M\&lsqb;\mathrm{\&psi;}\left(\left({\mathrm{\&phi;}}_0-{\mathrm{\&Delta;\&phi;}}_t\right)/j\right)\left(C_j+\underset{z\&Element;H,L,R}{\&Sigma;}\left(g_z{h}_z+h_z{g}_z+h_z{h}_z\right)C_j\right)/k\&times;k\&rsqb;\end{array}$

根据干涉相位频率提取相位矢量，可以设置阈值(窗口频率变换主瓣值的均值)，当估计窗口内提取得到的频率值大于该门限时，由于窗口中频谱变化最大值表示了相位变化的主方向，发明中使用窗口内信号频谱最大值作为矢量阈值：

对于高频部分的处理为：H,L,R表示三个方向(水平、垂直、视线)的正方向高频g_z和低频h_z，在方向搜索可时同时顾及这三个方向的反方向。若此时根据要控制的方向个数(2个、4个、6个或8个)对各方向求矢量和。

在实际离散相位图像算法中，相位高频部分通过如下公式进行计算：

$\begin{array}{c}\begin{array}{cc}\underset{z\&Element;H,L,R}{\&Sigma;}(g_z{h}_z+h_z{g}_z+h_z{h}_z)=\underset{z\&Element;H,L,R}{\&Sigma;}\left(2^{-j}\left(2^{-j}{D}_z\left(m,n\right)-D_z\left(m-\left(k-1\right)/2,n-\left(k-1\right)/2\right)\right)C_j\left(=m,n\right)\right)& j=1,2,\mathrm{3...}J\end{array}\\ +\left(2^{-j}\left(2^{-j}{D}_z\left(m,n\right)-D_z\left(m-\left(k-1\right)/2,n-\left(k-1\right)/2\right)\right)C_j\left(m,n\right)\right)\\ +\left(\begin{array}{c}{2}^{-j}\left(2^{-j}{D}_z\left(m,n\right)-D_z\left(m-\left(k-1\right)/2,n-\left(k-1\right)/2\right)\right)\\ \left(2^{-j}{D}_z\left(m,n\right)-D_z\left(m-\left(k-1\right)/2,n-\left(k-1\right)/2\right)\right)\end{array}\right)\end{array}$

上式中(m,n)为估计窗口中相位点，D_z(m-(k-1)/2,n-(k-1)/2)表示(m,n)邻域的相位点。

(k-1)/2保证相位变换点在估计窗口范围内。

步骤6：在估计窗口内利用高频部分矢量和构造非线性相位曲面模型。

由公式(2)Φ＝φ_linear+φ_non-linear，根据各层分解分量高频矢量和构造顾及相位细节变化的构建非线性相位曲面模型。

设式(8)中为信号φ_i分解的i层主矢量，那么主矢量和的相位反映非线性相位φ_non得到公式(10)：

${\mathrm{\&phi;}}_{non-linear}(m,n)=\arg \left(\underset{j=1}{\overset{J}{\&Sigma;}}{\mathrm{\&phi;}}_{G_M}\right)---\left(10\right)$

同时对局部窗口相位高频矢量求和保证了非线性曲面的特性，如公式(11)：

其中，上式中右边表示对j层高频的水平、垂直和视线方向分量求矢量和。表示视线向为主矢量方向，水平、垂直方向矢量和用获取；若水平、垂直为主矢量方向，视线向矢量的矢量和分别用获取。

步骤7：对小波分解相位进行逆变换，可求取非线性相位估计模型。

在估计窗口(k,l)内可以得到非线性相位：

$\begin{array}{c}{\mathrm{\&psi;}}_{k,l}\left(m,n\right)=\arg \left({\mathrm{WT}}^{-1}\left(f_{x}m+f_{y}n\right)\right)\\ =\arg \left(\underset{k,l}{\&Sigma;}{\mathrm{WT}}^{-1}\left(f_{x}m+f_{y}n\right){\mathrm{\&phi;}}_{G_0}\left(m,n\right)\exp \left(-j2\mathrm{\&pi;}\left(f_{x}m+f_{y}n\right)\right)\right)\end{array}---\left(12\right)$

其中，WT^-1(·)为小波逆变换，φ_Lp为小波分解低频部分，为小波分解高频部分在实际离散相位图像算法中，对于估计窗口的第t次逆变换由第t+1次分解得到，t+1次的逆变换窗口结构为：

${\left[\begin{array}{ccccc}{X}_{t+1}^{-1}\left(0\right)& X_{t+1}^{-1}\left(1\right)& \mathrm{...}& X_{t+1}^{-1}(p-1)& 0\\ 0& X_{t+1}^{-1}(p-1)& X_{t+1}^{-1}\left(0\right)& X_{t+1}^{-1}\left(1\right)& \mathrm{...}\\ 0& X_{t+1}^{-1}(p-1)& \mathrm{...}& X_{t+1}^{-1}\left(1\right)& X_{t+1}^{-1}\left(0\right)\\ X_{t+1}^{-1}(p-1)& X_{t+1}^{-1}\left(1\right)& X_{t+1}^{-1}\left(0\right)& 0& \mathrm{...}\end{array}\right]}_{0<p\&le;\frac{K-1}{2}\&times;\frac{L-1}{2}}$

X^-1表示C_J或D_H,L,R的逆变换。

步骤8：判断是否搜索完整个图像，是结束流程；否则，窗口向相位极值方向滑动。

以上借助具体实施例对本发明做了进一步描述，但是应该理解的是，这里具体的描述，不应理解为对本发明的实质和范围的限定，本领域内的普通技术人员在阅读本说明书后对上述实施例做出的各种修改，都属于本发明所保护的范围。

Claims (1)

1.利用非线性矢量曲面构建InSAR相位图像模型的方法，其特征在于，包括以下步骤：

步骤一：将相位干涉图以复数的浮点形式读取到Matlab中，标记矩阵行、列大小为M、N，单位为像素；将干涉相位图划分大小为k×l的初始估计窗口，并计算值得信任的第一个估计窗口，将窗口内相位点(m,n)单频信号形式转换为三角函数弦分量的复数相位形式：φ(m,n)＝angle(exp[j(f_xm+f_yn)])＝cos(f_xm+f_yn)+jsin(f_xm+f_yn)，其中，f_x和f_y分别为相位点行向量和列向量的频率；

步骤二：在估计窗口内，利用离散小波变换对干涉图像正弦分量和余弦分量进行多层分解，将相位表示为低频部分和高频部分和的形式；

步骤三：在估计窗口内，利用相位点高频分量的最大值确定主矢量的方向，也即是下一步搜索窗口的方向；

步骤四：在估计窗口内，计算相位主矢量，利用高频部分主矢量平滑窗口内的相位频率；

步骤五：在估计窗口内，根据各层分解分量的高频矢量和构建顾及相位细节变化的非线性相位曲面模型，将相位表示成非线性相位和线性相位和的形式Φ＝φ_linear+φ_non-linear，对非线性部分利用对j层分解的三个高频方向求矢量和，其中j＝1,2…J，分别表示对j层高频的水平、垂直和视线方向分量求矢量和，矢量求和保证了非线性曲面的特性；

步骤六：对小波分解相位进行逆变换得到相位估计模型，在估计窗口k×l内，对相位点ψ_k,l(m,n)，利用可得到估计窗口内非线性相位估计模型，其中，WT^-1(·)为小波逆变换，φ_Lp为小波分解低频部分，为小波分解高频部分和。