CN103645476A - 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 - Google Patents

Abstract

一种合成孔径雷达差分干涉图序列的时空同质滤波方法，它有八大步骤：步骤一：主副图像的自适应精密配准及副图像的重采样；步骤二：生成每张雷达图像的强度图；步骤三：生成序列差分干涉图；步骤四：确定滤波搜索窗口和统计邻域窗口的大小；步骤五：计算差分干涉图中每个像元的强度累积分布函数；步骤六：确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元；步骤七：相位同质的度量；步骤八：差分干涉图的空间域滤波。本发明能有效地同时滤除差分干涉图中的乘性斑点噪声和加性白噪声，并能很好地保留雷达图像的原始分辨率，从而可以大大地提高差分干涉图时序分析的有效性和精度。

Description

一种合成孔径雷达差分干涉图序列的时空同质滤波方法

技术领域

本发明涉及一种合成孔径雷达差分干涉图序列的时空同质滤波方法，属于合成孔径雷达干涉技术领域。它可以有效地同时滤除差分干涉图中的乘性斑点噪声和加性白噪声，并能很好地保留雷达图像的原始分辨率，从而可以大大地提高差分干涉图时序分析的有效性和精度。

背景技术

同所有的其它图像数据一样，合成孔径雷达差分干涉图也会受到各种噪声的污染。为了消除噪声，从被噪声污染了的图像中识别出真实图像，需要用滤波器对图像进行滤波。人们常用低通滤波器对图像进行滤波，在时域里，其实就是事先定义一个适当大小的小窗口（boxcar），让小窗口在图像上滑动，对这个小窗口中所有像元的数值进行平均运算，将平均值赋予该窗口所在位置的中心点像元。用数学表达式来表示的话，可以写成：

$\hat{u}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{u}_i---\left(1\right)$

上式中，表示在窗口内所有（N个）像元（值）组成的N维向量（u₁，u₂,……，u_N）^T的均值。显然，这种滤波器只适用于常数图像的加性噪声的滤波，其期望值就等于其均值，也即图像常数。但在实际情况下，图像都是非常数的。这时，方差可以被减小了N倍：

$\mathrm{Var}\left(U_k\right)=\frac{{\mathrm{\σ}}^2}{N}---\left(2\right)$

上式中，用随机变量U_k表示窗口内序号为k的像元，Var表示方差。

方差虽然被减小了N倍，但均值却是有偏的，且可以是非常大：

$E\left(U_k\right)-{\hat{u}}_k=\frac{1}{N}\underset{l\≠k}{\mathrm{\Σ}}({\hat{u}}_l-{\hat{u}}_k)---\left(3\right)$

这种有偏估计显然会极大地破坏图像的原始分辨率。

雷达差分干涉图可以表示为：

$\sqrt{I_m{I}_s}\exp \left(\mathrm{j\φ}\right)---\left(14\right)$

上式中，I_m和I_s分别表示主副图像的雷达反射强度（Intensity），φ表示主副图像去除参考面及高程相位后的相位差，即差分干涉相位。如果把雷达反射强度和差分干涉相位都看成是随机变量V和φ的话，那么雷达图像的强度将服从伽玛分布（GammaDistribution）：

$p\left(v\right|\hat{u})=\frac{L^L{v}^{L-1}{e}^{-L\frac{v}{\hat{u}}}}{\mathrm{\Γ}\left(L\right){\hat{u}}^L}---\left(5\right)$

上式中，v是雷达图像强度的观测值，

是的V均值，L是雷达图像的多视数，Γ是伽玛函数。雷达图像强度的数学期望为：

$E\left(V\right)=\hat{u}---\left(6\right)$

方差为：

$\mathrm{Var}\left(V\right)=\frac{{\hat{u}}^2}{L}---\left(7\right)$

显然，雷达强度信号的噪声是一种与信号相关的乘性噪声，要去除这类噪声，线性滤波器是无能为力的。

雷达差分干涉图的相位服从高斯分布：

上式中，

是相位φ的观测值，是相位均值，σ是其方差。

为了将雷达差分干涉图中由强度乘性噪声，雷达专业中称之为斑点噪声（speckle），和相位高斯噪声消除，本发明设计出了一种基于噪声特性的滤波方法。在雷达差分干涉技术的实际运用中，常常是需要对一序列差分干涉图进行时序分析，因此，本发明利用分布在时间域上的序列雷达图像的强度信息去除斑点噪声，再利用空间域上的相位信息去除高斯噪声，故称之为雷达差分干涉图的时空同质滤波。所谓图像滤波，就是对图像中每一个像元要用一个新的值去替代它原来的值，这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口，如图（1）中粗黑线框所示。对搜索窗口中心像元s滤波，即是用搜索窗口内所有像元的加权平均来替代其原来的值：

${\hat{u}}_s=\frac{\underset{t\∈\mathrm{\Ω}}{\mathrm{\Σ}}\mathrm{weight}(s,t)u_t}{\underset{t\∈\mathrm{\Ω}}{\mathrm{\Σ}}\mathrm{weight}(s,t)}---\left(9\right)$

上式中，Ω是以像元s为中心的一个规则窗口，

为差分干涉图像在像元s处滤波后的值，u_t表示像元t处的值，weight(s,t)表示像元t对于像元s所具有的权。本发明的亮点之处在于定义了判定像元t与像元s是同质的两个条件，即它们的强度分布和差分相位分布必须是相同的。也就是说，在滤波搜索窗口内，只有同质点才有资格参与滤波。在本发明中，强度分布的检验是在时域中进行的，而差分相位分布的检验是在空间域，即在每一个像元周边一个相对较小的规则窗口内进行的，这个较小的规则窗口称为统计邻域，或简称为邻域，见图1中阴影所示。

发明内容

1.目的：本发明的目的是提供一种合成孔径雷达差分干涉图序列的时空同质滤波方法。它克服了现有滤波技术的不足，能在保持住原始雷达图像分辨率的情况下，滤除雷达斑点噪声和差分干涉相位中的高斯噪声，从而可获取高分辨率的雷达差分干涉相位，可以大大提高后续差分相位时序分析的有效性和精度

2.技术方案：本发明是一种合成孔径雷达差分干涉图序列的时空同质滤波方法，该方法具体步骤如下：

步骤一：主副图像的自适应精密配准及副图像的重采样

在主副图像进行配准时，不仅要考虑雷达干涉系统（成像几何，干涉几何，雷达运行轨道等等）引起的地面同一目标在主副图像上的位置的偏移，还要考虑地面目标因较大的高程或较大的形变位移引起的主副图像上的位置偏移。在本发明中，副图像利用Sinc函数进行重采样，Sinc函数的位移量由下式决定：

Offset(i,j)=(△i,△j)=((△i,△j)_i,j)_System+((△i,△j)_i,j)_Object   (10)上式中，offset(i,j)表示了副图像中位于第i行第j列的像元相对于主图像上同一像元位置的偏移量，它是系统偏移和目标本身因高程或形变引起的偏移量的和。offset(i,j)由副图像s中位于(i,j)的像元的邻域与相应的主图像p中像元的邻域的协方差函数取最大值时的位移量决定，见下面的(11)和(12)式

$c_{\mathrm{ps}}{(n,m)}_{i.j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\Σ}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\Σ}}}(p(k,l)-\tilde{p})\·(s(k+n,l+m)-\tilde{s})---\left(11\right)$

$\left\{\begin{array}{c}\mathrm{\Δi}=n_{\max }\\ \mathrm{\Δj}=m_{\max }\end{array}\right.---\left(12\right)$

上两式中，p和s分别表示主副图像，c表示在主副图像中以第i行第j列为中心的N行M列的两个窗口的协方差函数，n和m分别表示行和列偏移量，k和l表示窗口中参与相关计算像元的行列位置，(△i，△j)表示当c取得最大值时的行列号，即（nmax，mmax）。

步骤二：生成每张雷达图像的强度图

单视复数雷达图像的每一个像元都是复数，它包含了两个信息，即雷达目标的后向散射强度，简称强度信息，和表征雷达与目标距离的相位信息。为了计算像元强度的统计特征，要先计算出每一个像元的强度：

$I_{i,j}={\mathrm{re}}_{i,j}^2+{\mathrm{im}}_{i,j}^2---\left(13\right)$

上式中，I_i,j表示雷达图像中位于第i行第j列像元的强度值，re表示像元的实部值，im表示像元的虚部值。

步骤三：生成序列差分干涉图

$\begin{array}{c}{f}_{\mathrm{ms}}=(f_m\·f_s^{*})\·\exp (-{\mathrm{\φ}}_{\mathrm{sim}})\\ =\sqrt{I_m{I}_s}\exp \left(\mathrm{j\φ}\right)\end{array}---\left(14\right)$

上式中，f_m，f_s和f_ms分别表示主图像，副图像和它们进行差分干涉后生成的差分干涉图，φ_sim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位，I_m和I_s分别表示主副图像的强度，φ表示差分干涉相位。当雷达图像较多，多于30张时，采用一张主图像的模式，否则采用多主图像模式，也就是所谓的短基线模式。

步骤四：确定滤波搜索窗口和统计邻域窗口的大小

所谓图像滤波，就是对图像中每一个像元要用一个新的值去替代它原来的值，这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口，通常取21行21列，或31行31列，如图1.中黑线框所示。为了计算差分相位的统计特性，需要定义一个像元的统计邻域窗口，如图1.中阴影所示，通常取5x5,或7x7个像元。

步骤五：计算差分干涉图中每个像元的强度累积分布函数

首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈，如图2.所示。对每一个像元构造一个强度向量:

U_s=(u_s1,u_s2,……u_sN)^T   (15)上式中，U_s表示像元s的强度向量，u_s1,u_s2,……u_sN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序。然后计算每一个像元的强度累积分布函数：

$F_s\left(x\right)=\left\{\begin{array}{ccc}0& & x<u_{s1}\\ \frac{k}{N}& \mathrm{if}& u_{\mathrm{sk}}\&le;x<u_{s(k+1)}\\ 1& & x\&GreaterEqual;u_{\mathrm{sN}}\end{array}\right.---\left(16\right)$

上式中，F_s(x)表示像元s在N张雷达强度图中的累积分布函数。理论上，雷达强度服从Gamma分布，但在大数据应用中，这种参数分布，即条件似然函数，计算复杂，耗时长，故本发明中，用无参数累积分布函数来近似替代其理论分布函数。

步骤六：确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元这一步，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元，换句话说，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元具有相同，至少是相似的强度累积分布函数。如果F_s(x)和F_t(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数，它们之间的最大距离定义为这两个分布的相似性的量度，见下式，并如图3.所示：

$D=\underset{s,t\&Element;\mathrm{\&Omega;},\mathrm{allx}}{\max }|F_s\left(x\right)-F_t\left(x\right)|---\left(17\right)$

由统计检验理论得知，两个累积分布的距离D服从如下式所示的KS分布：

$\begin{array}{c}{Q}_{\mathrm{ks}}\left(\mathrm{\&lambda;}\right)=2\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{j-1}\&CenterDot;e^{-2j^2{\mathrm{\&lambda;}}^2}\\ \mathrm{\&lambda;}=(\sqrt{N/2}+0.12+0.11/\sqrt{N/2})\&CenterDot;D\end{array}---\left(18\right)$

KS分布Q_ks(λ)具有如下的性质：

Q_ks(0)=1及Q_ks(∞)=0   （19）可见，Q_ks(λ)是一个单调函数，如果两个累积分布函数的距离D对应的Q_ks(λ)大于某一个阈值，即D是足够的小，则认为被比较的累积分布是相似的，用到两个像元的强度的比较上，则认为这两个像元在强度上是同质的。图4.是搜索强度同质像元后其结果的一个示意图，在以像元i为中心的搜索窗口内有j，k，l，m，n，o，p等6个像元与它是强度意义下的同质像元，这些像元构成搜索窗口内以像元i相关联的强度同质区域，它是搜索窗口Ω的一个子集Ω’，这个子集可以是不规则排列的。每一个像元都具有它自己的强度同质子集，并且对堆栈中所有差分干涉图都有效。

步骤七：相位同质的度量

通过前面的步骤，差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集，斑点噪声在这个子集内被去除了。步骤七，是要确定，在这个子集内，在差分相位概率分布的意义下，哪些像元与搜索窗口中心像元是同质的呢？在强度同质子集内，差分干涉相位服从高斯分布，它的噪声是一个均值为零的高斯噪声。居于此，本发明选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量，用W来表示：

$W_{s,t}=\exp [-\frac{\underset{k\&Element;\mathrm{\&omega;}}{\mathrm{\&Sigma;}}{({\mathrm{\&phi;}}_{s,k}-{\mathrm{\&phi;}}_{t,k})}^2}{{2\mathrm{\&sigma;}}^2}]---\left(20\right)$

上式中，W_s,t可以看作是像元s和t之间差分相位分布的相似性，ω是每一个像元的统计邻域，如图4.中同质像元周边的阴影所示，通常是一个5x5或7x7的小窗口。φ_s,k和φ_t,k分别是像元s和t在其各自邻域ω内对应像元k的差分相位，σ是邻域内噪声的方差。可见：

0<W_s,t≤1   （21）

步骤八：差分干涉图的空间域滤波

上一步骤中，像元间差分相位分布相似的程度可以用（20）式中的W_s,t来度量，亦可以理解为它们在差分相位意义上同质性的度量，W_s,t还具有（21）式的性质，因此，本发明利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权（weight）：

$u_s=\frac{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)\&CenterDot;u_t}{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)}---\left(22\right)$

上式中，u_s是差分干涉图中像元s的复数值，Ω'_s是属于像元s的强度同质子集，u_t是像元s的强度同质子集Ω'_s内的像元t的复数值，Weight(s,t)是像元t对于像元s所具有的权，它等于（20）式中的W_s,t。

作为总结，可以说，本发明提出的一个序列的差分干涉图的滤波同其它图像滤波方法的基础是一样的，都是利用了被滤波的像元的周边像元的加权平均。但本发明的特点在于，它没有将这个周边的形状固定，而是自适应地由周边那些在强度和差分相位概率分布意义下的同质像元所确定。因此，本发明提出的滤波方法，不仅可以同时消除差分干涉图中的乘性和加性噪声，同时还在很大程度上保留了雷达图像的原始分辨率。

附图说明

图1.搜索窗口和统计邻域窗口示意图，每一个小方格表示一个像元。搜索窗口是定义在以某一个像元，比如s，为中心的一个规则区域（M行，N列），统计邻域窗口是指在搜索窗口内某一像元，比如k，其周边一个较小的规则区域（m行，n列），如阴影所示。

图2.差分干涉图的堆栈及像元强度向量的结构示意图。

图3.累积分布函数间的距离定义示意图。

图4.在搜索窗口内的强度同质子集分布示意图。

图5.范家坪滑坡TerraSAR-X序列雷达图像强度均值图

图6.a范家坪滑坡20120107-20120129原始差分干涉图

图6.b范家坪滑坡20120107-20120129时空滤波后的差分干涉图

图7.本发明数据处理流程框图。

图中符号说明如下：

图1.中：s表示搜索窗口中心像元，k表示位于搜索窗口中的像元，它将被检验是否与像元s强度同质。

图2.中：N表示构成堆栈的差分干涉图数量，time表示按时间顺序堆栈。

图3.中：D表示两个累积分布函数之间的最大距离。

图4.中：i表示搜索窗口中的中心像元，j，k，l，m，n，o，及p表示搜索窗口内与i强度同质的像元。阴影表示强度同质像元的为下一步估计相位同质统计的邻域窗口。

具体实施方式

以长江三峡大坝地区范家坪滑坡地区INSAR形变监测中的差分干涉图的滤波为例，说明本发明在实际工程应用中的具体操作步骤。见图7，本发明是一种合成孔径雷达差分干涉图序列的时空同质滤波方法，该方法具体步骤如下：

步骤一：主副图像的自适应精密配准及副图像的重采样

为利用INSAR技术监测长江三峡大坝地区范家坪滑坡的非稳定特性，搜集了德国TerraSAR-X卫星X波段3米分辨率的21景单视复数（Single Look Complex）数据，时间跨距是从2012年1月至2013年2月。首先要选取一幅图像为主图像，作为参考图像，原则上它应该位于时域(时间基线集)和空间域(空间基线集)的中心。为了后续时序分析方便，这里选取第一张图像，即2012年1月7号的数据为主图像。以(1)式为模型基础，利用(2)式和(3)式表示的协方差方法计算其余各副图像中所有像元相对于主图像的相对偏移向量，用移位的Sinc函数对各副图像进行自适应二维重采样，以达到精密配准的目的，精度为十分之一个像元。

Offset(i,j)=(△i,△j)=((△i,△j)_i,j)_System+((△i,△j)_i,j)_Object   (1)

$c_{\mathrm{ps}}{(n,m)}_{i.j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\&Sigma;}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\&Sigma;}}}(p(k,l)-\tilde{p})\&CenterDot;(s(k+n,l+m)-\tilde{s})---\left(2\right)$

$\left\{\begin{array}{c}\mathrm{\&Delta;i}=n_{\max }\\ \mathrm{\&Delta;j}=m_{\max }\end{array}\right.---\left(3\right)$

式（1）中的第一项表示主图像中位于第i行第j列的像元与副图像中的同名点的位置的偏差量，该偏差量是一个二维向量，即式(1)中第二项。偏差向量是由系统产生的分量和主图像中该像元对应的地面目标本身形变产生的分量构成，分别为式（1）中第三和第四项所示。

在卫星轨道数和雷达成像时间等参数已知的情况下，可以先计算出系统产生的偏移分量。以该系统分量为基础，可以大致确定副图像中对应的像元位置。这样，对于主图像的中每一个像元（i,j），皆可以找到一个以该像元为中心的一个适当大小的窗口，比如512行和512列，同时亦可以确定副图像中对应的该窗口的位置。计算这两个窗口的协方差函数，如式（2）所示。M，N为窗口的行列数，p和s分别表示主副图像。使得协方差函数取得最大值时的m和n即为主图像中像元（i,j）在副图像中的偏移量，如式（3）所示。求出了所有副图像所有像元相对与主图像的偏移向量后再利用移位Sinc函数对所有幅图像进行重采样，使得所有幅图像皆与主图像精密配准。在该例中，感兴趣的地区仅是范家坪滑坡，故将所有21张雷达图像剪裁出以范家坪滑坡为中心的小区域图像，其包含2040（列）x1100（行）个像元。

步骤二：生成每张雷达图像的强度图

单视复数雷达图像的每一个像元都是复数，它包含了两个信息，即雷达目标的后向散射强度，简称强度信息，和表征雷达与目标距离的相位信息。为了计算像元强度的统计特征，要计算出每张图像中每一个像元的强度：

$I_{i,j}={\mathrm{re}}_{i,j}^2+{\mathrm{im}}_{i,j}^2---\left(4\right)$

上式中，I_i,j表示雷达图像中，位于第i行第j列像元的强度值，re表示像元的实部值，im表示像元的虚部值。为了检查21张雷达图像配准后的效果，可以计算21张强度图像的均值图像，如图5.，图像清晰程度可以大致显示配准效果。

步骤三：以2012年1月7日图像为主图像，其余为副图像，生成序列差分干涉图

$\begin{array}{c}{f}_{\mathrm{ms}}=(f_m\&CenterDot;f_s^{*})\&CenterDot;\exp (-{\mathrm{\&phi;}}_{\mathrm{sim}})\\ =\sqrt{I_m{I}_s}\exp \left(\mathrm{j\&phi;}\right)\end{array}---\left(5\right)$

上式中，f_m，f_s和f_ms分别表示主图像，副图像和它们进行差分干涉后生成的差分干涉图，φ_sim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位。利用TerraSAR-X卫星的WGS84精密轨道数据和地球WGS84椭球体（其表面即为INSAR系统参考面）参数，根据雷达成像时间及成像参数，可以计算出参考面及参考面上地形的（DEM）引起的干涉相位φ_sim。根据（5）式即可计算出差分干涉图。I_m和I_s分别表示主副图像的强度，φ表示差分干涉相位。

步骤四：确定滤波搜索窗口和统计邻域窗口的大小

所谓图像滤波，就是对图像中每一个像元要用一个新的值去替代它原来的值，这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口，通常取21行21列，或31行31列。现实采用的TerraSAR-X数据像元大小为3米x3米，搜索窗口取21x21已够大，搜索范围已达3600m².为了计算差分相位的统计特性，需要定义一个像元的统计邻域窗口，如图1.中阴影所示，这里取7x7个像元，以保证统计特性有较多样本单元。

步骤五：计算差分干涉图中每个像元的强度累积分布函数

像元位置在强度图和差分干涉图中是相同的，故差分干涉图中像元的累积分别函数是在序列强度图中进行的。首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈，如图2.所示。对每一个像元构造一个强度向量:

U_s=(u_s1,u_s2,……u_sN)^T   (6)上式中，U_s表示像元s的强度向量，u_s1,u_s2,……u_sN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序。然后计算每一个像元的强度累积分布函数：

$F_s\left(x\right)=\left\{\begin{array}{ccc}0& & x<u_{s1}\\ \frac{k}{N}& \mathrm{if}& u_{\mathrm{sk}}\&le;x<u_{s(k+1)}\\ 1& & x\&GreaterEqual;u_{\mathrm{sN}}\end{array}\right.---\left(7\right)$

上式中，F_s(x)表示像元s在N张雷达强度图中的累积分布函数，x是累积分布函数的自变量，如图3.所示。理论上，雷达强度服从Gamma分布，但在大数据应用中，这种参数分布，即条件似然函数，计算复杂，耗时长，故本发明中，用无参数累积分布函数以近似的方式替代其理论分布函数。这一步的输出，是差分干涉图中的每一个像元的累积分布函数F_s(x)。

步骤六：确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元

这一步，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元，换句话说，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元具有相同，至少是相似的强度累积分布函数。如果F_s(x)和F_t(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数，它们之间的最大距离定义为这两个分布的相似性的量度，见下式，并如图3.所示：

$D=\underset{s,t\&Element;\mathrm{\&Omega;},\mathrm{allx}}{\max }|F_s\left(x\right)-F_t\left(x\right)|---\left(8\right)$

在差分干涉图中，对每一个像元定义了一个以该像元为中心的大小相同的搜索窗口，其大小在步骤四中选定为21x21,共有441个像元，那么，除了该像元本身外，在其余440个像元中，哪些像元与该像元强度同质呢？这就要比较这两个像元的累积分布函数的距离。

由统计检验理论得知，两个累积分布的距离D服从如下式所示的KS分布：

$\begin{array}{c}{Q}_{\mathrm{ks}}\left(\mathrm{\&lambda;}\right)=2\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{j-1}\&CenterDot;e^{-2j^2{\mathrm{\&lambda;}}^2}\\ \mathrm{\&lambda;}=(\sqrt{N/2}+0.12+0.11/\sqrt{N/2})\&CenterDot;D\end{array}---\left(9\right)$

KS分布Q_ks(λ)具有如下的性质：

Q_ks(0)=1及Q_ks(∞)=0   （10）可见，Q_ks(λ)是一个单调函数，如果两个累积分布函数的距离D对应的Q_ks(λ)大于某一个阈值，即D是足够的小，则认为被比较的累积分布是相似的，用到两个像元的强度的比较上，则认为这两个像元在强度上是同质的。图4.是搜索强度同质像元后其结果的一个示意图，在以像元i为中心的搜索窗口内有j，k，l，m，n，o，p等6个像元与它是强度意义下的同质像元，这些像元构成搜索窗口内以像元i相关联的强度同质区域，它是搜索窗口Ω的一个子集Ω’，这个子集可以是不规则排列的。每一个像元都具有它自己的强度同质子集，并且对堆栈中所有差分干涉图都有效。

根据以上描述，对于搜索窗口内某一像元t，搜索出它的累积分布函数F_t(x)与中心像元s的累积分布函数F_s(x)的距离D，如图3.所示，将其带入（9）式，计算出对应的Q_ks(λ)，如果其值大于或等于阈值，就确定像元t是中性像元s在强度概率分布意义下的同质像元。阈值通常取0.95或更高。

这一步的输出，是确定了差分干涉图中每一个像元具有的在搜索窗口Ω内的强度同质像元的集合Ω'。

步骤七：相位同质的度量

通过前面的步骤，差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集，斑点噪声在这个子集内被去除了。接下来，是要确定，在这个子集内，在差分相位概率分布的意义下，哪些像元与搜索窗口中心像元是同质的呢？在强度同质子集内，差分干涉相位服从高斯分布，它的噪声是一个均值为零的高斯噪声。居于此，选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量，用W来表示：

$W_{s,t}=\exp [-\frac{\underset{k\&Element;\mathrm{\&omega;}}{\mathrm{\&Sigma;}}{({\mathrm{\&phi;}}_{s,k}-{\mathrm{\&phi;}}_{t,k})}^2}{{2\mathrm{\&sigma;}}^2}]---\left(11\right)$

在式（11）中，W_s,t可以看作是像元s和t之间差分相位分布的相似性，ω是每一个像元的统计邻域，如图4.中同质像元周边的阴影所示，通常是一个5x5或7x7的小窗口。φ_s,k和φ_t,k分别是像元s和t在其各自邻域ω内对应像元k的差分相位。需要强调的是，式（11）中求和运算是在邻域窗口内进行的。另外，与式（20）不同，因邻域内噪声的方差σ未知，在实际应用时，用滤波参数h来替代，它起到一个控制滤波程度的作用。由（11）式的指数函数性质可知：

0<W_s,t≤1   （12）

根据W_s,t在式（11）和（12）的性质，选用它作为对像元s进行相位滤波时像元t的权。这一步的输出，即是利用（11）式，计算在强度同质子集内，所有强度同质像元与搜索窗口中心像元间的相位分布的相似性。

步骤八：差分干涉图的空间域滤波

上一步骤中，像元间差分相位分布相似的程度可以用（11）式中的W_s,t来度量，亦可以理解为它们在差分相位意义上同质性的度量，W_s,t还具有（12）式的性质，因此，利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权（weight）：

${\hat{u}}_s=\frac{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)\&CenterDot;u_t}{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)}---\left(13\right)$

上式中，是差分干涉图中像元s的滤波后的复数值，Ω'_s是属于像元s的强度同质子集，u_t是像元s的强度同质子集Ω'_s内的像元t的复数值，Weight(s,t)是像元t对于像元s所具有的权，它等于（11）式中的W_s,t。

对每一张待滤波的差分干涉图，对每一个像元s，用它的搜索窗口内所有强度同质像元t

的像元值u_t的加权平均值

来代替像元s原来的值u_s。像元t对像元s的权按式(11)计算，从而完成差分干涉图的时空同质滤波。

图6.a是以范家坪滑坡TerraSAR-X的2012年1月7日图像和2012年1月29日图像干涉后的差分干涉图，差分干涉相位几乎被各种噪声淹没。图6.b是图6.a经过本发明的时空同质滤波以后的结果，可以见到十分清晰的高分辨率的雷达差分干涉相位。从黑到白的灰度变化显示了一个以2π为周期的差分相位变化。

可以看出，本发明提出的一个序列的差分干涉图的滤波同其它图像滤波方法的基础是一样的，都是利用了被滤波的像元的周边像元的加权平均。但本发明的特点在于，它没有将这个周边的形状固定，而是自适应地由周边那些在强度和差分相位概率分布意义下的同质像元所确定。因此，本发明提出的滤波方法，不仅可以同时消除差分干涉图中的乘性和加性噪声，同时还在很大程度上保留了雷达图像的原始分辨率。

Claims (1)

1.一种合成孔径雷达差分干涉图序列的时空同质滤波方法，其特征在于：该方法具体步骤如下：

步骤一：主副图像的自适应精密配准及副图像的重采样

在主副图像进行配准时，不仅要考虑雷达干涉系统的成像几何、干涉几何、雷达运行轨道引起的地面同一目标在主副图像上的位置的偏移，还要考虑地面目标因较大的高程或较大的形变位移引起的主副图像上的位置偏移；在此，副图像利用Sinc函数进行重采样，Sinc函数的位移量由下式决定：

Offset(i,j)=(△i,△j)=((△i,△j)_i,j)_System+((△i,△j)_i,j)_Object   (10)上式中，offset(i,j)表示副图像中位于第i行第j列的像元相对于主图像上同一像元位置的偏移量，它是系统偏移和目标本身因高程或形变引起的偏移量的和；offset(i,j)由副图像s中位于(i,j)的像元的邻域与相应的主图像p中像元的邻域的协方差函数取最大值时的位移量决定，见下面的(11)和(12)式

$c_{\mathrm{ps}}{(n,m)}_{i.j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\&Sigma;}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\&Sigma;}}}(p(k,l)-\tilde{p})\&CenterDot;(s(k+n,l+m)-\tilde{s})---\left(11\right)$

$\left\{\begin{array}{c}\mathrm{\&Delta;i}=n_{\max }\\ \mathrm{\&Delta;j}=m_{\max }\end{array}\right.---\left(12\right)$

上两式中，p和s分别表示主副图像，c表示在主副图像中以第i行第j列为中心的N行M列的两个窗口的协方差函数，n和m分别表示行和列偏移量，k和l表示窗口中参与相关计算像元的行列位置，(△i，△j)表示当c取得最大值时的行列号，即（nmax，mmax）；

步骤二：生成每张雷达图像的强度图

单视复数雷达图像的每一个像元都是复数，它包含了两个信息，即雷达目标的后向散射强度，简称强度信息，和表征雷达与目标距离的相位信息；为了计算像元强度的统计特征，要先计算出每一个像元的强度：

$I_{i,j}={\mathrm{re}}_{i,j}^2+{\mathrm{im}}_{i,j}^2---\left(13\right)$ 上式中，I_i,j表示雷达图像中位于第i行第j列像元的强度值，re表示像元的实部值，im表示像元的虚部值；

步骤三：生成序列差分干涉图

$\begin{array}{c}{f}_{\mathrm{ms}}=(f_m\&CenterDot;f_s^{*})\&CenterDot;\exp (-{\mathrm{\&phi;}}_{\mathrm{sim}})\\ =\sqrt{I_m{I}_s}\exp \left(\mathrm{j\&phi;}\right)\end{array}---\left(14\right)$

上式中，f_m，f_s和f_ms分别表示主图像，副图像和它们进行差分干涉后生成的差分干涉图，φ_sim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位，I_m和I_s分别表示主副图像的强度，φ表示差分干涉相位；当雷达图像较多，多于30张时，采用一张主图像的模式，否则采用多主图像模式，也就是所谓的短基线模式；

步骤四：确定滤波搜索窗口和统计邻域窗口的大小

所谓图像滤波，就是对图像中每一个像元要用一个新的值去替代它原来的值，这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果，这个窗口叫滤波搜索窗口，通常取21行21列，或31行31列；为了计算差分相位的统计特性，需要定义一个像元的统计邻域窗口，通常取5x5,或7x7个像元；

步骤五：计算差分干涉图中每个像元的强度累积分布函数

首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈，对每一个像元构造一个强度向量:

U_s=(u_s1,u_s2,……u_sN)^T   (15)上式中，U_s表示像元s的强度向量，u_s1,u_s2,……u_sN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序，然后计算每一个像元的强度累积分布函数：

$F_s\left(x\right)=\left\{\begin{array}{ccc}0& & x<u_{s1}\\ \frac{k}{N}& \mathrm{if}& u_{\mathrm{sk}}\&le;x<u_{s(k+1)}\\ 1& & x\&GreaterEqual;u_{\mathrm{sN}}\end{array}\right.---\left(16\right)$

上式中，F_s(x)表示像元s在N张雷达强度图中的累积分布函数，理论上，雷达强度服从Gamma分布，但在大数据应用中，这种参数分布，即条件似然函数，计算复杂，耗时长，故用无参数累积分布函数来近似替代其理论分布函数；

步骤六：确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元

这一步，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元，换句话说，是要确定，在一个搜索窗口内，哪些像元与该搜索窗口中心像元具有相同，至少是相似的强度累积分布函数；如果F_s(x)和F_t(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数，它们之间的最大距离定义为这两个分布的相似性的量度，见下式，

$D=\underset{s,t\&Element;\mathrm{\&Omega;},\mathrm{allx}}{\max }|F_s\left(x\right)-F_t\left(x\right)|---\left(17\right)$ 由统计检验理论得知，两个累积分布的距离D服从如下式所示的KS分布：

$\begin{array}{c}{Q}_{\mathrm{ks}}\left(\mathrm{\&lambda;}\right)=2\underset{j=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{j-1}\&CenterDot;e^{-2j^2{\mathrm{\&lambda;}}^2}\\ \mathrm{\&lambda;}=(\sqrt{N/2}+0.12+0.11/\sqrt{N/2})\&CenterDot;D\end{array}---\left(18\right)$ KS分布Q_ks(λ)具有如下的性质：

Q_ks(0)=1及Q_ks(∞)=0   （19）可见，Q_ks(λ)是一个单调函数，如果两个累积分布函数的距离D对应的Q_ks(λ)大于某一个阈值，即D是足够的小，则认为被比较的累积分布是相似的，用到两个像元的强度的比较上，则认为这两个像元在强度上是同质的；在以像元i为中心的搜索窗口内有j，k，l，m，n，o，p6个像元与它是强度意义下的同质像元，这些像元构成搜索窗口内以像元i相关联的强度同质区域，它是搜索窗口Ω的一个子集Ω’，这个子集是不规则排列的，每一个像元都具有它自己的强度同质子集，并且对堆栈中所有差分干涉图都有效；

步骤七：相位同质的度量

通过前面的步骤，差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集，斑点噪声在这个子集内被去除了；步骤七是要确定，在这个子集内，在差分相位概率分布的意义下，哪些像元与搜索窗口中心像元是同质的，在强度同质子集内，差分干涉相位服从高斯分布，它的噪声是一个均值为零的高斯噪声；居于此，选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量，用W来表示：

$W_{s,t}=\exp [-\frac{\underset{k\&Element;\mathrm{\&omega;}}{\mathrm{\&Sigma;}}{({\mathrm{\&phi;}}_{s,k}-{\mathrm{\&phi;}}_{t,k})}^2}{{2\mathrm{\&sigma;}}^2}]---\left(20\right)$ 上式中，W_s,t看作是像元s和t之间差分相位分布的相似性，ω是每一个像元的统计邻域，通常是一个5x5或7x7的小窗口；φ_s,k和φ_t,k分别是像元s和t在其各自邻域ω内对应像元k的差分相位，σ是邻域内噪声的方差；可见：

0<W_s,t≤1   （21）

步骤八：差分干涉图的空间域滤波

上一步骤中，像元间差分相位分布相似的程度用（20）式中的W_s,t来度量，亦理解为它们在差分相位意义上同质性的度量，W_s,t还具有（21）式的性质，因此，利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权（weight）：

$u_s=\frac{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)\&CenterDot;u_t}{\underset{t\&Subset;{{\mathrm{\&Omega;}}^{\&prime;}}_s}{\mathrm{\&Sigma;}}\mathrm{Weight}(s,t)}---\left(22\right)$ 上式中，u_s是差分干涉图中像元s的复数值，Ω'_s是属于像元s的强度同质子集，u_t是像元s的强度同质子集Ω'_s内的像元t的复数值，Weight(s,t)是像元t对于像元s所具有的权，它等于（20）式中的W_s,t。

Images