CN106097292B - 一种用于遥感影像处理的时空邻域高斯加权中值滤波方法 - Google Patents

Abstract

本发明公开了一种用于遥感影像处理的时空邻域高斯加权中值滤波方法，根据不同时间序列SAR影像的观测值同时与其时间邻域和空间邻域的观测值高度相关，且时间间隔越大、空间距离越大相关性越小，提出根据其时间和空间距离进行高斯加权的中值滤波算法，并根据SAR影像的统计分布特征，预估中值所在的区间，缩小中值搜索范围，加快中值求解速度，本发明能较好保持影像上地物细节特征，抑制斑点噪声，同时较大程度去除时序影像中车、船等临时地物的影响，解决了时序影像斑点噪声抑制问题。

Description

一种用于遥感影像处理的时空邻域高斯加权中值滤波方法

【技术领域】

本发明属于遥感影像处理技术领域,是一种合成孔径雷达(Synthetic ApertureRadar,SAR)影像相干斑噪声抑制方法，涉及一种对重轨时间序列SAR影像中的相干斑噪声进行滤除的新方法，具体涉及一种应用于时序SAR影像斑点噪声抑制的时空邻域高斯加权中值滤波快速算法。

【背景技术】

合成孔径雷达不受天气地理和时间等因素的限制，能够对地面进行高分辨率成像，并且具有一定的穿透力，因而被广泛应用于军事侦察、资源探测、环境监测、测绘制图、地理国情监测等对地遥感应用中，随着多颗雷达卫星的成功发射，时序雷达影像成为当前遥感应用的一大重要方面。

由于SAR基于相干成像机理进行工作，导致SAR影像中存在严重的斑点噪声，对SAR图像的解译判读和信息提取带来非常大的影响。早期的SAR影像斑点噪声抑制是通过空间多视平均处理来实现的，该方法的主要缺点是牺牲了影像的空间分辨率，不适用于当前的高分辨率SAR系统。目前主要采用空域滤波方法来实现斑点噪声抑制，即一般是利用一个滑动窗口，对窗口内的像素进行邻域处理得到窗口中心点的像素值。而中值滤波是一种非线性滤波算法，能较好地抑制斑点噪声、椒盐噪声等非高斯噪声，保持影像细节信息，在SAR影像中具有广泛的应用。

现有技术中，中值滤波仅运用在空间邻域，部分算法引入了不同形状和大小的邻域窗口，并采用高斯等方法进行了加权，但此类方法的不足是：没有充分利用时间序列观测值之间的相关性，且计算效率较低。其它一些考虑时间域相关性的滤波方法，多对各时相同等对待，未充分顾及时间间隔越大其相关性越低的问题，且不能克服临时性车、船等着某一时相偶然出现地物的影响，因此滤波结果对后续分类解译效果的提升依然有限。这些改进算法均未能彻底解决时序SAR影像中斑点噪声滤除的问题。

本发明即针对现有技术的不足而研究提出。

【发明内容】

本发明一种用于遥感影像处理的时空邻域高斯加权中值滤波方法，充分利用SAR时间序列影像中时间邻域和空间邻域观测值的相关性，发挥中值滤波算法对斑点噪声等非高斯噪声抑制和细节信息保持的良好性能，在抑制斑点噪声的同时，抑制时序影像中临时地物的影响，有效解决现有方法未能很好解决时序SAR影像中斑点噪声快速滤波的问题。

为解决上述技术问题，本发明一种用于遥感影像处理的时空邻域高斯加权中值滤波方法,时序影像的表示形式为强度影像，宽和高分别为n_w和n_h个像素，共有n_t个时相，其中第k个时相的观测时刻记为T_k，通过滑动窗口遍历每个像素(i,j,k)，其中i、j表示影像的行列号，记滤波前像素(i,j,k)的强度影像为P(i,j,k)，对其时空邻域进行处理即得到像素(i,j,k)的滤波结果其特征在于包括以下步骤：

步骤1：设置时空邻域窗口参数：即将空间邻域窗口的大小设置为W_s，时间邻域窗口的大小为W_t，由此得到通过滑动窗口遍历每个像素(i,j,k)时的时空邻域大小为W_s×W_s×W_t,满足：51≥W_s≥3，且W_s为奇数；且W_t至少为最大相邻时相观测时间间隔的3倍；

步骤2：设置时空邻域高斯加权函数的参数：即空间邻域加权参数σ_s＞0，且2σ_s<W_s/2；时间邻域加权参数σ_t＞0，且2σ_t<W_t/2；

步骤3：初始化当前处理影像的时相，k＝1；

步骤4：计算当前时相滑动窗口的高斯权矩阵；

步骤4.1：计算当前时相滑动窗口的高斯权矩阵，先计算参与当前时相滤波的相邻时相，按照时间先后顺序遍历所有时相的观测时间T_t，当时相t与当前时相k的时间间隔T_t-T_k≤W_t/2时，将时相t加入时相k时间邻域集合，并记时相k的时间邻域集合为S_nt,k，S_nt,k中元素的个数，即为时相k时间邻域的时相个数，记为W_tk；

步骤:4.2：计算当前处理时相k的时空邻域滑动窗口内的三维权矩阵各元素的值w(d_x,d_y,d_t),满足：

其中d_x,d_y分别是方形模板上(x,y)处到模板中心的行方向和列方向的距离，其中x,y＝1,2,...,W_s；d_t为时相k的邻近时相t的观测时刻到时相k的时间间隔，记为d_t＝T_t-T_k，且T_t∈S_nt,k，t＝1,2,...,W_tk，由此得到大小为W_s×W_s×W_tk的三维权矩阵w₀；

步骤4.3：计算三维权矩阵元素的总和，记总权值为w_all，则有：

其中w(i,j,k)是三维权矩阵w₀第i行第j列第k个时相的权重；

步骤4.4：计算高斯权矩阵，记为w,则有

步骤5：初始化当前处理像素(i_c,j_c,k)为第一个像素，即(i_c,j_c)＝(0,0)，并将滑动窗口的中心与(i_c,j_c,k)对齐，开始处理第一行第一列像素；

步骤6：根据时空邻域内最大值和最小值确定中值的搜索范围为

步骤6.1：记当前处理像素时空邻域W_s×W_s×W_tk个影像强度值组成的集合为S_N，计算S_N的最大值P_max,c和最小值P_min,c，W_tk为时相个数；

步骤6.2：计算当前处理像素时空邻域内的影像幅度值的最大值和最小值

步骤6.3：SAR幅度影像的视数为n_L，由于视数n_L越大，越接近正态分布，其中值与均值越接近，当n_L＞100，幅度影像近似服从正态分布，其中值在均值附近的概率极大；随着视数n_L降低，其中值逐渐向均值的左边移动，缩小后的中值搜索范围可根据不同的视数n_L进行计算：

步骤7：判断中值是否位于缩小的搜索范围

步骤7.1：判断中值位于三者中的何区间。

步骤7.2：计算满足和时空邻域元素的总权重和

步骤7.3：当时，中值位于区间并记w₁＝0，记位于区间的元素组成的集合为S；否则当时，中值位于区间并记记位于区间的元素组成的集合为S；否则中值位于区间并记w₂＝0，记位于区间的元素组成的集合为S；

步骤8：在步骤7确定的搜索区间J＝(J_min,J_max)内进行二分法递归快速搜索高斯时空加权的中值,其中J_min、J_max分别是步骤7.3中区间J的下限值和上限值；

步骤8.1：终止条件是集合S的元素个数为1，该元素即为求得的中值满足终止条件时，进入步骤9；不满足终止条件时，继续执行步骤8.2；

步骤8.2：根据中值所在搜索区间J＝(J_min,J_max)，将该区间内的元素分为两部分，并分别记为S₁和S₂，满足：S₁<(J_min+J_max)/2，S₂≥(J_min+J_max)/2，则S₁和S₂的权值分别为△₁和△₂；

步骤8.3：如果w₁+△₁<0.5,则更新w₁＝w₁+△₁，更新S＝S₂，更新J＝((J_min+J_max)/2,J_max)；否则更新w₂＝w₂+△₂，更新S＝S₁，更新J＝(J_min,(J_min+J_max)/2)；

步骤8.4:继续执行步骤8.1；

步骤9：保存结果

步骤10：当前时相像素全部处理完成，进入步骤11；否则移动滑动窗口中心到下一个像素，更新当前处理像素(i_c,j_c,k)，进入步骤6；

步骤11：当时相k影像的全部像素都处理完时，更新当前处理时相k＝k+1；如果k≤n_t进入步骤4，开始对下一时相进行处理，否则处理结束。

与现有技术相比较，本发明一种用于遥感影像处理的时空邻域高斯加权中值滤波方法，采用上述步骤，充分利用SAR时间序列影像中时间邻域和空间邻域观测值的相关性，发挥中值滤波算法对斑点噪声等非高斯噪声抑制和细节信息保持的良好性能，在抑制斑点噪声的同时，抑制时序影像中临时地物的影响，有效解决现有方法未能很好解决时序SAR影像中斑点噪声快速滤波的问题。

【附图说明】

下面结合附图对本发明的具体实施方式作进一步详细说明，其中：

图1为本发明的流程图。

【具体实施方式】

下面结合附图对本发明的实施方式作详细说明。

如图1所示，本发明一种用于遥感影像处理的时空邻域高斯加权中值滤波方法，要求时序影像为经过精配准的强度影像，配准精度不低于0.1像素，而且数据形式为经过定标的浮点型强度数据。各相邻时相的间隔一般比较短，可以在几小时到数日甚至数月，但考虑到间隔越久，影像的相关性越差，地物发生了较大变化，不满足观测同一地物的假设条件，斑点噪声抑制和细节保持性能则会下降，本发明优选时间间隔最大不超过半年。

下面以华南某区域的岗哨的6月到11月共12个时相的影像进行处理，经过了4视处理，像元大小为12m×12m，图像大小为1400像素×5000像素。除8月14日到9月7日两时相时间间隔为24天，其余影像相邻间隔均为12天。采用步骤如下：

步骤1：设置时空邻域窗口参数，即将空间邻域窗口的大小设置为W_s＝9像素，时间邻域窗口的大小设置为W_t＝12*6＝72天。

步骤2：设置时空邻域高斯加权函数的参数：空间邻域加权参数为σ_s＝2像素，时间邻域加权参数为σ_t＝16天。

步骤3：初始化当前处理影像的时相，k＝1。

步骤4：计算当前时相滑动窗口的高斯权矩阵；

步骤4.1：计算当前时相滑动窗口的高斯权矩阵。先计算参与当前时相滤波的相邻时相，即按照时间先后顺序遍历所有时相的观测时间T_t，当时相t与当前时相k的时间间隔T_t-T_k≤36天时，将时相t加入时相k时间邻域集合，并记时相k的时间邻域集合为S_nt,k，S_nt,k中元素的个数即为时相k时间邻域的时相个数，记为W_tk，当k值不同时，W_tk＝4,5,6,7。

步骤:4.2：计算当前处理时相k的时空邻域滑动窗口内的三维权矩阵各元素的值w(d_x,d_y,d_t),满足：

其中d_x,d_y分别是方形模板上(x,y)处到模板中心(5,5)的行方向和列方向的距离，其中x,y＝1,2,...,9；d_t为时相k的邻近时相t的观测时刻到时相k的时间间隔d_t＝T_t-T_k，且T_t∈S_nt,k，t＝1,2,...,W_tk，由此可得大小为9×9×W_tk的权矩阵w₀。

步骤4.3：计算三维权矩阵元素的总和，记总权值为w_all：

其中w(i,j,k)是三维权矩阵w₀第i行第j列第k个时相的权重。

步骤4.4：计算高斯权矩阵：

步骤5：初始化当前处理像素(i_c,j_c,k)为第一个像素，即(i_c,j_c)＝(0,0)，并将滑动窗口的中心与(i_c,j_c,k)对齐，开始处理第一行第一列像素。

步骤6：根据时空邻域内最大值和最小值确定中值的搜索范围为

步骤6.1：记当前处理像素时空邻域W_s×W_s×W_tk个影像强度值组成的集合为S_N，计算S_N的最大值P_max,c和最小值P_min,c，W_tk为时相个数；

步骤6.2：计算当前处理像素时空邻域内的影像幅度值的最大值和最小值

步骤6.3：SAR幅度影像的视数为n_L，由于视数n_L越大，越接近正态分布，其中值与均值越接近，当n_L＞100，幅度影像近似服从正态分布，其中值在均值附近的概率极大；随着视数n_L降低，其中值逐渐向均值的左边移动，缩小后的中值搜索范围可根据不同的视数n_L进行计算：

步骤7：判断中值是否位于缩小的搜索范围

步骤7.1：判断中值位于三者中的何区间；

步骤7.2：计算和时空邻域元素的权重，且分别记为和

步骤7.3：当时，中值位于区间并记w₁＝0，记位于区间的元素组成的集合为S；否则当时，中值位于区间并记记位于区间的元素组成的集合为S；否则中值位于区间并记w₂＝0，记位于区间的元素组成的集合为S；

步骤8：在步骤7确定的搜索区间J＝(J_min,J_max)内进行二分法递归快速搜索高斯时空加权的中值,其中J_min、J_max分别是步骤7.3中区间J的下限值和上限值；

步骤8.1：终止条件是集合S的元素个数为1，该元素即为求得的中值满足终止条件时，进入步骤9；不满足终止条件时，继续执行步骤8.2；

步骤8.2：根据中值所在搜索区间J＝(J_min,J_max)，将该区间内的元素分为两部分，并分别记为S₁和S₂，满足：S₁<(J_min+J_max)/2，S₂≥(J_min+J_max)/2，则S₁和S₂的权值分别为△₁和△₂；

步骤8.3：如果w₁+△₁<0.5,则更新w₁＝w₁+△₁，更新S＝S₂，更新J＝((J_min+J_max)/2,J_max)；否则更新w₂＝w₂+△₂，更新S＝S₁，更新J＝(J_min,(J_min+J_max)/2)；

步骤8.4:继续执行步骤8.1；

步骤9：保存结果

步骤10：当前时相像素全部处理完成，进入步骤11；否则移动滑动窗口中心到下一个像素，更新当前处理像素(i_c,j_c,k)，进入步骤6；

步骤11：当时相k影像的全部像素都处理完时，更新当前处理时相k＝k+1；如果k≤n_t进入步骤4，开始对下一时相进行处理，否则处理结束。

本发明采用上述步骤，除了应用于SAR影像外，对于光学遥感影像均具有一定的噪声抑制效果，与现有技术相比本发明方法以下优势：

(1)由于是非线性滤波器，且引入了多时相观测值，因此具有更好的斑点噪声抑制性能；

(2)由于并非所有像素的值都能影响滤波输出，只有像素的权重和大小位置能影响输出，因此能较大程度上抑制车、船等临时地物的影响，滤波后影像尽可能保持了稳定的地物信息；

(3)由于利用了不同时相相关性不同的特点，时间间隔越近，相关性越高，同时做获得相同等效视数时，空间邻域窗口更小，且权重与空间距离成反比，因此具有更好的空间信息。

(4)由于基于SAR分布先验知识，缩小了均值搜索范围，大幅减小了比较与交换操作，因此与已有算法相比，具有更好的计算效率，做实施例中，计算时间约为非快速算法的1/4。

Claims (1)

1.一种用于遥感影像处理的时空邻域高斯加权中值滤波方法,时序影像的表示形式为强度影像，宽和高分别为n_w和n_h个像素，共有n_t个时相，其中第k个时相的观测时刻记为T_k，通过滑动窗口遍历每个像素(i,j,k)，其中i、j表示影像的行列号，记滤波前像素(i,j,k)的强度影像为P(i,j,k)，对其时空邻域进行处理即得到像素(i,j,k)的滤波结果其特征在于包括以下步骤：

步骤1：设置时空邻域窗口参数：即将空间邻域窗口的大小设置为W_s，时间邻域窗口的大小为W_t，由此得到通过滑动窗口遍历每个像素(i,j,k)时的时空邻域大小为W_s×W_s×W_t,满足：51≥W_s≥3，且W_s为奇数；且W_t至少为最大相邻时相观测时间间隔的3倍；

步骤2：设置时空邻域高斯加权函数的参数：即空间邻域加权参数σ_s＞0，且2σ_s<W_s/2；时间邻域加权参数σ_t＞0，且2σ_t<W_t/2；

步骤3：初始化当前处理影像的时相，k＝1；

步骤4：计算当前时相滑动窗口的高斯权矩阵；

步骤4.1：计算当前时相滑动窗口的高斯权矩阵，先计算参与当前时相滤波的相邻时相，按照时间先后顺序遍历所有时相的观测时间T_t，当时相t与当前时相k的时间间隔T_t-T_k≤W_t/2时，将时相t加入时相k时间邻域集合，并记时相k的时间邻域集合为S_nt,k，S_nt,k中元素的个数，即为时相k时间邻域的时相个数，记为W_tk；

步骤:4.2：计算当前处理时相k的时空邻域滑动窗口内的三维权矩阵各元素的值w(d_x,d_y,d_t),满足：

其中d_x,d_y分别是方形模板上(x,y)处到模板中心的行方向和列方向的距离，其中x,y＝1,2,...,W_s；d_t为时相k的邻近时相t的观测时刻到时相k的时间间隔，记为d_t＝T_t-T_k，且T_t∈S_nt,k，t＝1,2,...,W_tk，由此得到大小为W_s×W_s×W_tk的三维权矩阵w₀；

步骤4.3：计算三维权矩阵元素的总和，记总权值为w_all，则有：

其中w(i,j,k)是三维权矩阵w₀第i行第j列第k个时相的权重；

步骤4.4：计算高斯权矩阵，记为w,则有

步骤5：初始化当前处理像素(i_c,j_c,k)为第一个像素，即(i_c,j_c)＝(0,0)，并将滑动窗口的中心与(i_c,j_c,k)对齐，开始处理第一行第一列像素；

步骤6：根据时空邻域内最大值和最小值确定中值的搜索范围为

步骤6.1：记当前处理像素时空邻域W_s×W_s×W_tk个影像强度值组成的集合为S_N，计算S_N的最大值P_max,c和最小值P_min,c，W_tk为时相个数；

步骤6.2：计算当前处理像素时空邻域内的影像幅度值的最大值和最小值

步骤6.3：SAR幅度影像的视数为n_L，由于视数n_L越大，越接近正态分布，其中值与均值越接近，当n_L＞100，幅度影像近似服从正态分布，其中值在均值附近的概率极大；随着视数n_L降低，其中值逐渐向均值的左边移动，缩小后的中值搜索范围可根据不同的视数n_L进行计算：

步骤7：判断中值是否位于缩小的搜索范围

步骤7.1：判断中值位于三者中的何区间；

步骤7.2：计算满足和时空邻域元素的总权重和

步骤7.3：当时，中值位于区间并记w₁＝0，记位于区间的元素组成的集合为S；否则当时，中值位于区间并记记位于区间的元素组成的集合为S；否则中值位于区间并记w₂＝0，记位于区间的元素组成的集合为S；

步骤8：确定的搜索区间J＝(J_min,J_max)内进行二分法递归快速搜索高斯时空加权的中值,其中J_min、J_max分别是步骤7.3中区间J的下限值和上限值；

步骤8.1：终止条件是集合S的元素个数为1，该元素即为求得的中值满足终止条件时，进入步骤9；不满足终止条件时，继续执行步骤8.2；

步骤8.2：根据中值所在搜索区间J＝(J_min,J_max)，将该区间内的元素分为两部分，并分别记为S₁和S₂，满足：S₁<(J_min+J_max)/2，S₂≥(J_min+J_max)/2，则S₁和S₂的权值分别为△₁和△₂；

步骤8.3：如果w₁+△₁<0.5,则更新w₁＝w₁+△₁，更新S＝S₂，更新J＝((J_min+J_max)/2,J_max)；否则更新w₂＝w₂+△₂，更新S＝S₁，更新J＝(J_min,(J_min+J_max)/2)；

步骤8.4:继续执行步骤8.1；

步骤9：保存结果

步骤10：当前时相像素全部处理完成，进入步骤11；否则移动滑动窗口中心到下一个像素，更新当前处理像素(i_c,j_c,k)，进入步骤6；

步骤11：当时相k影像的全部像素都处理完时，更新当前处理时相k＝k+1；如果k≤n_t进入步骤4，开始对下一时相进行处理，否则处理结束。