本发明涉及遥感影像处理技术领域，是一种合成孔径雷达(Synthetic Aperture Radar,简称SAR)影像相干斑噪声抑制方法，具体为一种SAR影像的中值滤波快速近似算法。
背景技术：
合成孔径雷达不受天气地理和时间等因素的限制，能够对地面进行高分辨率成像，并且具有一定的穿透力，因而被广泛应用于军事侦察、资源探测、环境监测、测绘制图、地理国情监测等对地遥感应用中。随着多颗雷达卫星的成功发射，时序雷达影像成为当前遥感应用的一大重要方面。
由于SAR基于相干成像机理进行工作，导致SAR强度影像中存在严重的斑点噪声，对SAR强度影像的解译判读和信息提取带来非常大的影响。早期的SAR强度影像斑点噪声抑制是通过空间多视平均处理来实现的，该方法的主要缺点是牺牲了影像的空间分辨率，不适用于当前的高分辨率SAR系统。
目前，主要采用空域滤波方法来实现斑点噪声抑制，即一般是利用一个滑动窗口，对窗口内的象素进行邻域处理得到窗口中心点的象素值。其中中值滤波是一种非线性滤波算法，能较好地抑制斑点噪声、椒盐噪声等非高斯噪声，保持影像细节信息，在SAR强度影像中具有广泛的应用。
由于中值滤波需要进行排序，其中比较和交换操作的次数与元素个数的平方成正比，虽然提出了不少快速算法，但当窗口增大或参与运算的波段或时相增加时，效率改善仍然有限。
本发明即针对现有技术的不足而研究提出。
技术实现要素：
本发明要解决的技术问题是提供一种SAR影像的中值滤波快速近似算法，即解决SAR强度影像中进行多时相或大滑动窗中值滤波计算效率不高的问题，本发明根据多视SAR强度影像近似服从Gamma分布的特点，先设置邻域窗口参数，将强度影像转换为幅度影像，并统计幅度影像的最大值和最小值，然后估算幅度影像中值所处的区间，并基于该范围对幅度影像进行线性量化到较少的灰度级，以便通过计算滑动窗口分布累积函数实现幅度影像中值所在灰度级的准确定位，有效提高中值滤波的计算效率。
为解决上述技术问题，本发明一种SAR影像的中值滤波快速近似算法,影像的表示形式为强度影像，数据类型为浮点型，宽和高分别为nw和nh个像素。通过滑动窗口遍历每个像素(i,j)，其中(i,j)表示影像的行列号，记滤波前像素(i,j)的强度影像为P(i,j)，对其时空邻域进行处理即得到(i,j)的滤波结果包括以下步骤：
步骤1：设置邻域窗口参数,即设置邻域窗口为方形，边长为Ws，满足：51≥Ws≥3，且Ws为奇数，则滑动窗口中用于求中值的元素个数为nn＝Ws·Ws，对于多时相数据参与中值滤波，nn还应乘上时相个数。
步骤2：对整个强度影像进行开方处理，得到幅度影像。
步骤3：统计幅度影像的最大值Amax和最小值Amin。
步骤4：根据SAR幅度影像的分布，计算中值所在区间，并对幅度影像进行量化。
步骤4.1：根据SAR幅度影像的分布，计算中值所在区间：
步骤4.2：将幅度影像从0至n进行了量化，且500≥n≥100，优选n≈200：
步骤5：初始化当前处理像素(ic,jc)为第一个像素，即(ic,jc)＝(0,0)，并将滑动窗口的中心与(ic,jc)对齐，开始处理第一行第一列像素。
步骤6：计算滑动窗口内局部幅度值的最大值Amax,c和最小值Amin,c。
步骤7：如果Amin,c>Amin+0.35·(Amax-Amin)或Amax,c<Amin+0.10·(Amax-Amin)时，则认为滑动窗口内像素强度值位于最大值或最小值附近，需对幅度值根据滑动窗口内的最大值和最小值将浮点型的幅度数据重新量化为0到100的整数：
步骤8：遍历滑动窗口内的各像素，统计量化后的累积分布函数及各级强度均值。
步骤8.1：设置一个长度为n+1的计数器数组C，并初始化计数器数组C各元素值为0；设置一个长度为n+1的强度均值数组M，并初始化强度均值数组M各元素值为0。
步骤8.2：遍历滑动窗口内幅度量化值g(i,j),重新量化后的则基于g'(i,j)分别更新计数器数组C和强度均值数组M的元素:C[g(i,j)]＝C[g(i,j)]+1，M[g(i,j)]＝M[g(i,j)]+P(i,j)。
步骤8.3：将计数器数组C和强度均值数组M内的值转换成累积分布函数值：M[0]＝M[0]/C[0]，C[0]＝C[0]/nn，M[k]＝M[k]/C[k]，C[k]＝C[k]/nn+C[k-1],其中k＝1,2,...,n。
步骤9：遍历累积分布函数数组，即计数器数组C，当C[k]<0.5且C[k+1]≥0.5时，中值位于k+1灰度级对应范围。
步骤10：如果(C[k+1]-C[k])nn≤4且C[k+1]-C[k]≤0.1，取M[k]作为中值滤波的近似输出值；否则进入步骤12。
步骤11：对滑动窗口内g(i,j)＝＝k+1的强度值进行排序，取第(0.5-C[k])nn个作为中值滤波的输出。
步骤12：当前时相像素全部处理完成，处理结束；否则移动滑动窗口中心到下一个像素，更新当前处理像素(ic,jc)，进入步骤6。
与现有技术相比较，本发明一种SAR影像的中值滤波快速近似算法，通过上述步骤，缩小了幅度影像中值所在区间范围，并进行了量化，有助于准确对中值进行准确定位，大幅减小搜索计算量，同时取得了与经典中值滤波几乎相同的滤波结果，解决了大窗口和多时相SAR强度影像中值滤波效率低下的问题。
【附图说明】
下面结合附图对本发明的具体实施方式作进一步详细说明，其中：
图1为本发明流程图。
【具体实施方式】
下面结合附图对本发明的实施方式作详细说明。
以下将结合附图和具体实施例对本发明做进一步的阐述。
如图1所示，本发明一种SAR影像的中值滤波快速近似算法，实施例以其中某一区域的5个时相精配准后的强度影像中值滤波为例。该强度影像经过了4视处理，图像大小为5000像素×5000像素，处理在7×7的邻域窗口进行，求取7×7×5个时相(三维邻域)共245个元素的中值，具体包括以下处理步骤：
步骤1：设置邻域窗口参数，将邻域窗口设置为方形，边长Ws＝7。由于要计算5个时相的中值，故滑动窗口中用于求中值的元素个数为nn＝7·7·5。
步骤2：对整个强度影像进行开方处理，得到幅度影像。
步骤3：统计幅度影像的最大值Amax和最小值Amin。
步骤4：根据SAR幅度影像的分布，计算中值所在区间，并对幅度影像进行量化。
步骤4.1：计算中值所在区间，即：
步骤4.2：将幅度影像从0至n进行了量化，且500≥n≥100，优选n≈200：
步骤5：窗口中心像素为(ic,jc)，初始化窗口中心像素(ic,jc)，即(ic,jc)＝(0,0)，并将滑动窗口的中心与(ic,jc)对齐，开始处理第一行第一列像素。
步骤6：计算滑动窗口内局部幅度值的最大值Amax,c和最小值Amin,c。
步骤7：如果Amin,c>Amin+0.35·(Amax-Amin)或
Amax,c<Amin+0.10·(Amax-Amin)时，则认为滑动窗口内像素强度值位于最大值或最小值附近，需对幅度值根据滑动窗口内的最大值和最小值从0至100进行重新量化：
步骤8：遍历滑动窗口内的各像素，统计量化后的累积分布函数及各级强度均值。
步骤8.1：设置一个长度为n+1的计数器数组C，并初始化计数器数组C各元素值为0；设置一个长度为n+1的强度均值数组M，并初始化强度均值数组M各元素值为0。
步骤8.2：遍历滑动窗口内幅度量化值g(i,j),重新量化后的则基于g'(i,j)分别更新计数器数组C和强度均值数组M的元素:C[g(i,j)]＝C[g(i,j)]+1，M[g(i,j)]＝M[g(i,j)]+P(i,j)。
步骤8.3：将计数器数组C和强度均值数组M内的值转换成累积分布函数值：M[0]＝M[0]/C[0]，C[0]＝C[0]/nn，M[k]＝M[k]/C[k]，C[k]＝C[k]/nn+C[k-1],其中k＝1,2,...,n。
步骤9：遍历累积分布函数数组，即计数器数组C，当C[k]<0.5且C[k+1]≥0.5时，中值位于k+1灰度级对应范围。
步骤10：如果(C[k+1]-C[k])nn≤4且C[k+1]-C[k]≤0.1，取M[k]作为中值滤波的近似输出值；否则进入步骤12。
步骤11：对滑动窗口内g(i,j)＝＝k+1的强度值进行排序，取第(0.5-C[k])nn个作为中值滤波的输出。
步骤12：当前时相像素全部处理完成，处理结束；否则移动滑动窗口中心到下一个像素，更新当前处理像素(ic,jc)，进入步骤6。
本发明除了应用于SAR强度影像外，对于SAR幅度影像或者光学遥感影像均具有一定的噪声抑制效果。通过上述实施步骤，本实施例的结果如下：
(1)本实施例所得中值与精确中值存在误差，并对该误差与本中值的比值进行了评估，该比值的最大值为0.174，均值为0.031，则表明本发明中绝大部分像素的结果与精确均值一致。
(2)需要进行二次量化的滑窗个数占比为9％,需要进行步骤11精确排序的比例为1.6％。
(3)本实施例中耗时为其全排序方法的0.07倍，为快速排序算法的0.31倍。即由于本发明方法的耗时与需要排序的元素个数无显著关系，若待排序元素个数增加，耗时增速则缓慢，但全排序算法和快速算法耗时均呈窗口边长的几何级数增加，当参与中值滤波的滑窗内元素较大时，特别处理多时相数据或大滑窗时，故本发明效率提升更加明显。