时序sar时空邻域高斯加权中值滤波斑噪抑制快速算法

时序sar时空邻域高斯加权中值滤波斑噪抑制快速算法
【专利摘要】本发明公开了一种时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法，根据不同时间序列SAR影像的观测值同时与其时间邻域和空间邻域的观测值高度相关，且时间间隔越大、空间距离越大相关性越小，提出根据其时间和空间距离进行高斯加权的中值滤波算法，并根据SAR影像的统计分布特征，预估中值所在的区间，缩小中值搜索范围，加快中值求解速度，本发明能较好保持影像上地物细节特征，抑制斑点噪声，同时较大程度去除时序影像中车、船等临时地物的影响，解决了时序影像斑点噪声抑制问题。
【专利说明】时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法 【技术领域】
[0001 ] 本发明属于遥感影像处理技术领域，是一种合成孔径雷达（Synthetic Aperture Radar，SAR)影像相干斑噪声抑制方法，涉及一种对重轨时间序列SAR影像中的相干斑噪声 进行滤除的新方法，具体涉及一种应用于时序SAR影像斑点噪声抑制的时空邻域高斯加权 中值滤波快速算法。 【【背景技术】】
[0002] 合成孔径雷达不受天气地理和时间等因素的限制，能够对地面进行高分辨率成 像，并且具有一定的穿透力，因而被广泛应用于军事侦察、资源探测、环境监测、测绘制图、 地理国情监测等对地遥感应用中，随着多颗雷达卫星的成功发射，时序雷达影像成为当前 遥感应用的一大重要方面。
[0003] 由于SAR基于相干成像机理进行工作，导致SAR影像中存在严重的斑点噪声，对SAR 图像的解译判读和信息提取带来非常大的影响。早期的SAR影像斑点噪声抑制是通过空间 多视平均处理来实现的，该方法的主要缺点是牺牲了影像的空间分辨率，不适用于当前的 高分辨率SAR系统。目前主要采用空域滤波方法来实现斑点噪声抑制，即一般是利用一个滑 动窗口，对窗口内的像素进行邻域处理得到窗口中心点的像素值。而中值滤波是一种非线 性滤波算法，能较好地抑制斑点噪声、椒盐噪声等非高斯噪声，保持影像细节信息，在SAR影 像中具有广泛的应用。
[0004] 现有技术中，中值滤波仅运用在空间邻域，部分算法引入了不同形状和大小的邻 域窗口，并采用高斯等方法进行了加权，但此类方法的不足是:没有充分利用时间序列观测 值之间的相关性，且计算效率较低。其它一些考虑时间域相关性的滤波方法，多对各时相同 等对待，未充分顾及时间间隔越大其相关性越低的问题，且不能克服临时性车、船等着某一 时相偶然出现地物的影响，因此滤波结果对后续分类解译效果的提升依然有限。这些改进 算法均未能彻底解决时序SAR影像中斑点噪声滤除的问题。
[0005] 本发明即针对现有技术的不足而研究提出。 【
【发明内容】
】
[0006] 本发明一种时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法，充分利用SAR 时间序列影像中时间邻域和空间邻域观测值的相关性，发挥中值滤波算法对斑点噪声等非 高斯噪声抑制和细节信息保持的良好性能，在抑制斑点噪声的同时，抑制时序影像中临时 地物的影响，有效解决现有方法未能很好解决时序SAR影像中斑点噪声快速滤波的问题。
[0007] 为解决上述技术问题，本发明一种时序SAR时空邻域高斯加权中值滤波斑噪抑制 快速算法，时序影像的表示形式为强度影像，宽和高分别为n4Pn h个像素，共有nt个时相，其 中第k个时相的观测时刻记为Tk，通过滑动窗口遍历每个像素（i，j，k)，其中（i，j)表示影像 的行列号，记滤波前像素（i，j，k)的强度影像为P(i，j，k)，对其时空邻域进行处理即得到像 素(i，j，k)的滤波结果,包括以下具体步骤：
[0008] 步骤1:设置时空邻域窗口参数：即将空间邻域窗口的大小设置为Ws，时间邻域窗 口的大小为Wt，由此得到通过滑动窗口遍历每个像素（i，j，k)时的时空邻域大小为W s XWS X Wt，满足：51多Ws彡3，且1为奇数；两n； - ?；，且wt至少为最大相邻时相观测时间间隔的3 倍。
[0009] 步骤2:设置时空邻域高斯加权函数的参数：即空间邻域加权参数〇s>〇，且2〇s〈W s/ 2;时间邻域加权参数〇t>0，且2〇t〈Wt/2。
[0010 ] 步骤3:初始化当前处理影像的时相，k = 1。
[0011] 步骤4:计算当前时相滑动窗口的高斯权矩阵。
[0012] 步骤4.1:计算参与当前时相滤波的相邻时相，按照时间先后顺序遍历所有时相的 观测时间Tt，当时相t与当前时相k的时间间隔T t-Tk<Wt/2时，将时相t加入时相k时间邻域 集合，并记时相k的时间邻域集合为S nt,k，Snt,k的大小即为时间邻域，时相个数为W tk。
[0013] 步骤:4.2:计算当前处理时相k的时空邻域滑动窗口内的三维权矩阵各元素的值w (dx，dy，dt)，满足：
[0015] 其中dx，dy分别是方形模板上(x，y)处到模板中心I
的行方向和列方 向的距离，其中x，y=l，2,. . .，Ws;dt为时相k的邻近时相t的观测时刻到时相k的时间间隔dt =Tt-Tk，且Tt e Snt, k，t = 1，2，…，Wtk。由此得到大小为Ws X Ws X Wtk的权矩阵Wo。
[0016] 步骤4.3:计算三维权矩阵元素的总和，即总权值Wall满足：
[0017] Wa"= J； η·(/..Μ)，其中w(i，j，k)是三维权矩阵Wo第i行第j列第k个时相的权重。 ij\k=l
[0018] 步骤4 · 4:计算归一化权矩阵W，满足：
[0019] 步骤5:初始化当前处理像素（i。」。)为第一个像素，即（i。」。）= (0,0)，并将滑动 窗口的中心与（i。，」。)对齐，开始处理第一行第一列像素。
[0020] 步骤6:根据时空邻域内最大值和最小值确定，并缩小的中值搜索范围。
[0021] 步骤6.1:记当前处理像素时空邻域Ws X Ws X Wtk个影像强度值组成的集合为SN，计 算SN的最大值P max, c和最小值Pmin, c。
[0022] 步骤6.2:计算当前处理像素时空邻域内的影像幅度值的最大值4?^ 和 最小值= °
[0023]步骤6.3:由于SAR幅度影像的视数越大，越接近正态分布，其中值与均值越接近， 当nL>100，幅度影像近似服从正态分布，其中值在均值附近的概率极大；随着视数降低，其 中值逐渐向均值的左边移动。且通过实验数据统计表明，缩小后的中值搜索范围可根据不 同的视数nL进行计算：
[0030] 步骤7:判断中值是否位于缩小的搜索范围
[0031] 步骤7.1:判断中值位于(/^,力_)、dyU、三者中的何区 间，理想情况是位于区间(€η,,尤_)，且该区间的元素个数最少。
[0032] 步骤7.2:计算Ρ < 和/^时空邻域元素的权重，以及计算胃尸< 和
[0033] 步骤7 · 3 :当W(P?；之_) > 0.5时，中值位于区间/ = dd J，并记w! = 0， w2 = 1 - W(K之m.()，记位于区间d,J的元素组成的集合为3;否则当贾/^之_)>0.5 时，中值位于区间·/=(>：_，_)，并记％二W(P #_)，>v= 1 - W(K $nax>e)，记位于区 间dr ）的元素组成的集合为S ;否则中值位于区间《/=(之_, )，并记 in = w(p<4a、j，W2=0,记位于区间的元素组成的集合为So
[0034] 步骤8:在步骤7确定的搜索区间J= (Jmin，Jmax)内进行二分法递归快速搜索高斯时 空加权的中值。
[0035] 步骤8.1:终止条件是集合S的元素个数为1，该元素即为求得的中值身。满足终止 条件时，进入步骤9;不满足终止条件时，继续执行步骤8.2。
[0036] 步骤8.2:根据中值所在搜索区间1=(1_，1_)，将该区间内的元素分为两部分， 并分别记为S4PS2，满足:SK( Jmin+Jmax)/2，S2彡（Jmin+Jmax)/2，则S4PS2的权值分别为八4口 Δ2〇
[0037] 步骤8 · 3 :如果Wl+Al〈0 · 5，则更新Wl = Wl+Al，更新S = S2，更新 J = ( ( Jmin+Jmax)/2， Jmax);否则更新W2 = W2+A2，更新S = Si，更新 J = ( Jmin，（ Jmin+JmaX)/2) 〇
[0038] 步骤8.4:继续执行步骤8.1。
[0039] 步骤9:保存结果汽，:,./,4) = Α。
[0040] 步骤10:当前时相像素全部处理完成，进入步骤11;否则移动滑动窗口中心到下一 个像素，更新当前处理像素 (i。，jc)，进入步骤6。
[0041] 步骤11:当时相k影像的全部像素都处理完时，更新当前处理时相k = k+l，如果 m进入步骤4，开始对下一时相进行处理，否则处理结束。
[0042] 与现有技术相比较，本发明一种时序SAR时空邻域高斯加权中值滤波斑噪抑制快 速算法，采用上述步骤，充分利用SAR时间序列影像中时间邻域和空间邻域观测值的相关 性，发挥中值滤波算法对斑点噪声等非高斯噪声抑制和细节信息保持的良好性能，在抑制 斑点噪声的同时，抑制时序影像中临时地物的影响，有效解决现有方法未能很好解决时序 SAR影像中斑点噪声快速滤波的问题。 【【附图说明】】
[0043]下面结合附图对本发明的【具体实施方式】作进一步详细说明，其中：
[0044]图1为本发明的流程图。 【【具体实施方式】】
[0045]下面结合附图对本发明的实施方式作详细说明。
[0046]如图1所示，本发明一种时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法， 要求时序影像为经过精配准的强度影像，配准精度不低于〇. 1像素，而且数据形式为经过定 标的浮点型强度数据。各相邻时相的间隔一般比较短，可以在几小时到数日甚至数月，但考 虑到间隔越久，影像的相关性越差，地物发生了较大变化，不满足观测同一地物的假设条 件，斑点噪声抑制和细节保持性能则会下降，本发明优选时间间隔最大不超过半年。
[0047]下面以华南某区域的岗哨的6月到11月共12个时相的影像进行处理，经过了 4视处 理，像元大小为12m X 12m，图像大小为1400像素X 5000像素。除8月14日至1」9月7日两时相时 间间隔为24天，其余影像相邻间隔均为12天。采用步骤如下：
[0048] 步骤1:设置时空邻域窗口参数，即将空间邻域窗口的大小设置为Ws = 9像素，时间 邻域窗口的大小设置为Wt= 12*6 = 72天。
[0049] 步骤2:设置时空邻域高斯加权函数的参数:空间邻域加权参数为〇s = 2像素，时间 邻域加权参数为〇t= 16天。
[0050 ] 步骤3:初始化当前处理影像的时相，k = 1。
[0051 ]步骤4:计算当前时相滑动窗口的高斯权矩阵。
[0052]步骤4.1:计算参与当前时相滤波的相邻时相，即按照时间先后顺序遍历所有时相 的观测时间Tt，当时相t与当前时相k的时间间隔Tt-Tk<36天时，将时相t加入时相k时间邻 域集合，并记时相k的时间邻域集合为S nt,k，Snt,k的大小即时间邻域，时相个数为W tk，当k值 不同时，Wtk = 4,5,6,7。
[0053]步骤:4.2:计算当前处理时相k的时空邻域滑动窗口内的三维权矩阵各元素的值w (dx，dy，dt)，满足：
[0055] 其中dx，dy分别是方形模板上(X，y)处到模板中心(5，5)的行方向和列方向的距离， 其中x，y = l，2,. . .，9;dt为时相k的邻近时相t的观测时刻到时相k的时间间隔dt = Tt-Tk，且 Ttesnt,k，t = l，2，...，Wtk，由此可得大小为9X9XWtk的权矩阵Wo。
[0056] 步骤4.3:计算三维权矩阵元素的总和，即总权值wall: 9.9馬
[0057] Wk = Σ吨从），其中w(i J，k)是三维权矩阵WQ第i行第j列第k个时相的权重。
[0058] 步骤4.4:计算归一化权矩阵：
[0060]步骤5:初始化当前处理像素（i。，」。)为第一个像素，即（i。，」。）= (0,0)，并将滑动 窗口的中心与（i。，」。)对齐，开始处理第一行第一列像素。
[0061 ]步骤6:根据时空邻域内最大值和最小值确定，并缩小的中值搜索范围。
[0062]步骤6.1:记当前处理像素时空邻域9X9 XWtk个影像强度值组成的集合为SN，计算 Sn的取大值Pmax，c和取小值Pmin，c。
[0063]步骤6.2:计算当前处理像素时空邻域内的影像幅度值的最大值和 最小值= 。
[0064]步骤6.3:由于SAR幅度影像的视数越大，越接近正态分布，其中值与均值越接近， 此实施例中α = 4，故缩小后的中值搜索范围为： 毛in.t=(毛邮 + 〇·1 5 · (4m'r ))
[0065] 、2 Pmm,c + 0：·32 ' (4naXjc _ ν Ο
[0066] 步骤7:判断中值是否位于缩小的搜索范围
[0067] 步骤7.1:判断中值位于三者中的何区 间，理想情况是位于区间()Lu，，且该区间的元素个数最少。
[0068] 步骤7.2:计算iL,和尸< 时空邻域元素的权重，以及计算贾尸< 之￥)和 W("之).
[0069] 步骤 7 · 3 :当 W(i^ ) > 0.5 时，中值位于区间 / = (fnm,,之π,￡)，并记w: = 0， % = 1 - W(P2患_):，记位于区间dd,)的元素组成的集合为S;否则当之^) > 0J 时，中值位于区间J=(Uma、J，并记Μ = W(/><) * h2 = 1 - 之~)，记位于区间 (見ιη?，Α~)的元素组成的集合为S ;否则中值位于区间/=(之，并记 ^ = ，:W2=〇,记位于区间d,.，Cx,)的元素组成的集合为S。
[0070] 步骤8:在步骤7确定的搜索区间J= (Jmin，Jmax)内进行二分法递归快速搜索高斯时 空加权的中值。
[0071] 步骤8.1:终止条件是集合S的元素个数为1，该元素即为求得的中值#，满足终止 条件时，进入步骤9;不满足终止条件时，继续执行步骤8.2。
[0072] 步骤8.2:根据中值所在搜索区间1=(11"，1^)，将该区间内的元素分分别记为5 1 和 S2，满足:SK (Jmin+Jmax)/2，S2 彡（Jmin+Jmax)/2，则 SjPS2 的权值分别为 Λ4ΡΛ2。
[0073] 步骤8 · 3 :如果Wl + Al〈0 · 5则更新Wl=Wl + Al，更新S= S2，更新J = ( ( Jmin+Jmax)/2， Jmax);否则更新W2 = W2+A2，更新S = Si，更新 J = ( Jmin，（ Jmin+JmaX)/2) 〇
[0074] 步骤8.4:继续执行步骤8.1。
[0075] 步骤9:保存结果Α，:,./·，，幻=#。
[0076] 步骤10:当前时相像素全部处理完成，进入步骤11;否则移动滑动窗口中心到下一 个像素，更新当前处理像素(i。，jc)，进入步骤6。
[0077] 步骤11:当时相k影像的全部像素都处理完时，更新当前处理时相k = k+l;如果 m进入步骤4，开始对下一时相进行处理;否则处理结束。
[0078] 本发明采用上述步骤，除了应用于SAR强度影像外，对于SAR幅度影像或者光学遥 感影像均具有一定的噪声抑制效果，与现有技术相比本发明方法以下优势：
[0079] (1)由于是非线性滤波器，且引入了多时相观测值，因此具有更好的斑点噪声抑制 性能；
[0080] (2)由于并非所有像素的值都能影响滤波输出，只有像素的权重和大小位置能影 响输出，因此能较大程度上抑制车、船等临时地物的影响，滤波后影像尽可能保持了稳定的 地物信息；
[0081 ] (3)由于利用了不同时相相关性不同的特点，时间间隔越近，相关性越高，同时做 获得相同等效视数时，空间邻域窗口更小，且权重与空间距离成反比，因此具有更好的空间 信息。
[0082] (4)由于基于SAR强度分布先验知识，缩小了均值搜索范围，大幅减小了比较与交 换操作，因此与已有算法相比，具有更好的计算效率，做实施例中，计算时间约为非快速算 法的1/4。
【主权项】
1. 时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法，时序影像的表示形式为强 度影像，宽和高分别为nw和nh个像素，共有nt个时相，其中第k个时相的观测时刻记为化，通 过滑动窗口遍历每个像素（i，j，k)，其中（i，j)表示影像的行列号，记滤波前像素(i，j，k)的 强度影像为p(i，j，k)，对其时空邻域进行处理即得到像素 α，j，k)的滤波结果戶(/，./，/;)，其 特征在于包括W下步骤： 步骤1:设置时空邻域窗口参数:即将空间邻域窗口的大小设置为Ws，时间邻域窗口的大 小为Wt，由此得到通过滑动窗口遍历每个像素（i，j，k)时的时空邻域大小为WsXWsXWt，满 足:51 >Ws>3，且Ws为奇数；?? 7： - 7]，且Wt至少为最大相邻时相观测时间间隔的3倍； 步骤2:设置时空邻域高斯加权函数的参数：即空间邻域加权参数〇s〉0,且2〇s<Ws/2;时 间邻域加权参数〇t〉0，且2〇t<Wt/2; 步骤3:初始化当前处理影像的时相，k = 1; 步骤4:计算当前时相滑动窗口的高斯权矩阵； 步骤5:初始化当前处理像素(ic，jc)为第一个像素，即。。^'。）= (0,0)，并将滑动窗口的 中屯、与（ic，jc)对齐，开始处理第一行第一列像素； 步骤6:根据时空邻域内最大值和最小值确定，并缩小的中值捜索范围； 步骤7:判断中值是否位于缩小的捜索范围(iL,,); 步骤8:在步骤7确定的捜索区间J=(Jmin，J"ax)内进行二分法递归快速捜索高斯时空加 权的中值； 步骤9:保存结果约4,克，巧=#. 步骤10:当前时相像素全部处理完成，进入步骤11;否则移动滑动窗口中屯、到下一个像 素，更新当前处理像素（iC，jc)，进入步骤6; 步骤11:当时相k影像的全部像素都处理完时，更新当前处理时相k=k+l;如果k《nt进 入步骤4，开始对下一时相进行处理，否则处理结束。2. 根据权利要求1所述时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法，其特征 在于所述步骤4还包括： 步骤4.1:计算参与当前时相滤波的相邻时相，按照时间先后顺序遍历所有时相的观测 时间Tt，当时相t与当前时相k的时间间隔Tt-Tk《Wt/^2时，将时相t加入时相k时间邻域集合， 并记时相k的时间邻域集合为Snt,k，Snt,k的大小即为时间邻域，时相个数为Wtk; 步骤:4.2:计算当前处理时相k的时空邻域滑动窗口内的Ξ维权矩阵各元素的值w(dx, dy,山），满足：其中山，dy分别是方形模板上(X，y)处到模板中屯的行方向和列方向的距 离，其中x，y=l，2，...，Ws;dt为时相k的邻近时相t的观测时刻到时相k的时间间隔dt = Tt- Tk，且TteSnt,k，t = l，2,...，Wtk，由此得到大小为WsXWsXWtk的权矩阵Wo; 步骤4.3 :计算Ξ维权矩阵元素的总和，即总权值Wall满足：，其中w(i，j，k)是Ξ维权矩阵Wo第i行第j列第k个时相的权重； 步骤4.4:计算归一化权矩阵W，满足:W = ^3. 根据权利要求1所述时序SA則寸空邻域高斯加权中值滤波斑噪抑制快速算法，其特征 在于所述步骤6包括： 步骤6.1:记当前处理像素时空邻域WsXWsXWtk个影像强度值组成的集合为Sn，计算Sn 的最大值Pmax, C和最小值Pmin, C ; 步骤6.2:计算当前处理像素时空邻域内的影像幅度值的最大值和最小值步骤6.3:由于SAR幅度影像的视数越大，越接近正态分布，其中值与均值越接近，当nL〉 100，幅度影像近似服从正态分布，其中值在均值附近的概率极大;随着视数降低，其中值逐 渐向均值的左边移动，且通过实验数据统计表明，缩小后的中值捜索范围可根据不同的视 数nL进行计算：4. 根据权利要求1所述时序SA則寸空邻域高斯加权中值滤波斑噪抑制快速算法，其特征 在于所述步骤7包括： 步骤7.1 :判断中值位于Wnin.点η ,)、（么in.:，iU )、（么a、,。，f胃.。)亏者中的何区间，理 想情况是位于区间(么1。,，/Lw )，且该区间的元素个数最少； 步骤7 . 2 :计算F<么,。.^和F <么。、。时空邻域元素的权重，W及计算W(F空iU。)和 W(尸先,)； 步骤7 . 3 :当W(P含么。.1.)>化5时，中值位于区间/ =巧。心息。,> 并记W1 = 0， W2=l-W(/j<i么山记位于区间巧。的元素组成的集合为S;否则当W(/^么。、,.)>0.5 时，中值位于区间/=(么。，并记= W(戶 < 么。X.:.)，w'2 = 1 - W(戶?Ξ么。，记位于区间 的元素组成的集合为S ;否则中值位于区间/=(么，并记 M'l = W(P<么。、>W2 = 0,记位于区间(iLu.，气g、,。)的元素组成的集合为S。5.根据权利要求1所述时序SAR时空邻域高斯加权中值滤波斑噪抑制快速算法，其特征 在于所述步骤8包括： 步骤8.1:终止条件是集合S的元素个数为1，该元素即为求得的中值食,满足终止条件 时，进入步骤9;不满足终止条件时，继续执行步骤8.2; 步骤8.2:根据中值所在捜索区间1=^1。^。3〇，将该区间内的元素分为两部分，并分 别记为Si和S2，满足:Sl<( Jmin+Jmax)/2，S2> ( Jmin+Jmax)/2，则Si和S2的权值分别为Δ?和Δ2 ; 步骤8.3 :如果￥1+Δ：1<0.5 ,则更新W1 = Wl+Δ?，更新S = S2，更新 J = ( ( Jmin+Jmax)/2 , Jmax); 否则更新化二'WS+As ,更新S = Si,更新 J = ( Jmin , ( Jmin+Jmax)/2); 步骤8.4:继续执行步骤8.1。
【文档编号】G06T5/50GK106097292SQ201610640069
【公开日】2016年11月9日
【申请日】2016年8月4日 公开号201610640069.X, CN 106097292 A, CN 106097292A, CN 201610640069, CN-A-106097292, CN106097292 A, CN106097292A, CN201610640069, CN201610640069.X
【发明人】邓少平, 甘宗平, 孙盛, 李胜, 吴泽洪, 方志民, 刘学林
【申请人】甘宗平