CN112561901A - 基于有噪谱平稳的sar图像显著区域检测方法 - Google Patents

Abstract

本发明公开了一种基于有噪谱平稳的SAR图像显著区域检测方法，所述方法如下：步骤一：输入标准图像格式的高分辨率SAR图像；步骤二：计算有噪估计条件下局部图像谱平稳特征，最终生成与图像等大小的有噪估计谱平稳特征图；步骤三：使用误差为e的数据格式无关p‑tile分割法处理有噪估计谱平稳特征图，取分割前景部分为疑似显著区域；步骤四：获取疑似显著区域的连通域结果，利用尺度滤波滤除像素数少于c的连通域图块，余下部分作为高分辨率SAR图像显著区域输出。本发明具有较好的性能与效率，能够适用于任意格式数据，具有较高的实用价值。

Description

基于有噪谱平稳的SAR图像显著区域检测方法

技术领域

本发明属于遥感图像处理技术领域，涉及一种SAR图像显著区域检测方法，具体涉及一种基于有噪谱平稳的SAR图像显著区域检测方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar，SAR)是一种具有全天时、全天候特性的主动式微波雷达，是遥感领域中一种不可或缺的重要手段，其在城镇遥感、灾害监测、海洋探索、人类活动痕迹观测等涉及到国防民生的应用中起着关键作用。与传统光学遥感图像成像机制不同，SAR图像主要反映了地面物体的后向散射信息，这一特点使得SAR图像对于一些人造目标或者其他异常、地面突变突起等目标或区域比较敏感，因此相对于光学图像，针对这一类地物，即显著性地物，使用SAR图像进行检测具有明显优势。同时，显著性地物虽然不能直接代表特定类别人造地物，但是往往包含了较多的有价值目标或者信息。因此显著性区域的检测提取在遥感图像处理中，对于高价值目标快速检测、冗余信息去除等具有重要意义，同时也为进一步的精细处理以及其他应用提供高质量数据。但是，SAR图像受到相干斑噪声干扰，使得大多经典图像处理方法、经典细节描述特征在SAR图像中不再适用，导致SAR图像的处理质量下降，这种影响在高分辨率SAR图像中更加明显。

目前针对高分辨率SAR图像中显著目标、显著区域的检测问题，主流解决方案是根据SAR图像性质直接构建显著性特征，并利用分割方法从特征图中提取显著性区域。这一方法中的难点是显著性特征的构建。现有显著性特征一般从光学显著性特征发展而来，例如中心-邻域对比度方法。但是SAR图像局部细节及纹理与地表粗糙度、地面物体分布相关，且由于相干斑噪声影响，使得部分非显著区域，例如草地或沙地背景，在图像中呈现出了高散射、复杂纹理的形式，使其经常被误认为是显著性区域。现有针对这一问题的解决方法，主要通过引入复杂的统计量来滤除由于背景散射能量较高以及复杂纹理的影响，但是当背景对比度过高或者纹理模式过于复杂时，统计量的滤除能力有限。

另外，针对特征图的分割方法中，大量经典分割方法主要针对标准图像格式设计，例如要求输入图像数据格式为8bit量化或固定比特数整形量化。尤其是直方图分割法，需要统计图像像素分布，而这一步骤通过扫描所有图像量化值并进行统计实现。但是大多特征图在未进行量化调整时难以满足特定量化要求，这将直接导致两个问题，一是进行格式转换、整形量化时可能导致信息与细节丢失、分辨率降低；而不经格式转化直接进行统计，则排序步骤会消耗大量时间，降低了显著区域检测作为中间处理过程的实用性，增加了SAR图像处理链路负担。此外，不同特征图数据分布区间差异大，难以通过设置经验性硬阈值进行数据筛选。

发明内容

本发明的目的是提供一种基于有噪谱平稳的SAR图像显著区域检测方法，该方法利用噪声拟合构建了频域稳定性特征，并设计了数据格式无关的e精度p-tile分割方法，以较高的处理效率实现了高分辨率SAR图像中显著区域的检测，有效地滤除了图像中具有高散射复杂纹理的背景区域，避免了非显著性复杂背景的干扰，为进一步的精细处理以及其他应用提供高质量的参考数据。

本发明的目的是通过以下技术方案实现的：

一种基于有噪谱平稳的SAR图像显著区域检测方法，包括如下步骤：

步骤一：输入标准图像格式的高分辨率SAR图像；

步骤二：利用尺寸为(2d+1)×(2d+1)的滑动窗口对高分辨率SAR图像进行局部数据采集，计算有噪估计条件下局部图像谱平稳特征，最终生成与图像等大小的有噪估计谱平稳特征图，其中：d为正整数，用来控制滑动窗口尺寸；

步骤三：使用误差为e的数据格式无关p-tile分割法(Format-free e-accuracyp-tile segmentation,FF ep-tile)处理有噪估计谱平稳特征图，取分割前景部分为疑似显著区域；

步骤四：获取疑似显著区域的连通域结果，利用尺度滤波滤除像素数少于c的连通域图块，余下部分作为高分辨率SAR图像显著区域输出。

相比于现有技术，本发明具有如下优点：

本发明的方法作为高分辨率SAR图像粗处理中关键步骤，具有较好的性能与效率，能够适用于任意格式数据，具有较高的实用价值。

附图说明

图1是本发明方法的整体流程图；

图2是MiniSAR实验图像，分辨率0.1m×0.1m，机场环境，8bit量化Bitmap图像；

图3是d＝5、R＝3.2时有噪估计谱平稳特征数据分布图；

图4是d＝5、R＝3.2时有噪估计谱平稳特征数据线性动态压缩8bit量化图像；

图5是d＝5、R＝3.2时有噪估计谱平稳特征数据对数动态压缩8bit量化图像；

图6是d＝5、R＝3.2、e＝0.1％、p＝5％的特征图分割结果；

图7是d＝5、R＝3.2、e＝0.1％、p＝5％时SAR图像显著区域检测结果；

图8是d＝5、R＝10、e＝0.1％、p＝10％的特征图分割结果；

图9是d＝5、R＝10、e＝0.1％、p＝10％时SAR图像显著区域检测结果。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

本发明提供了一种基于有噪估计谱平稳特征的高分辨率SAR图像显著区域检测方法，如图1所示，所述方法具体实施步骤如下：

步骤一：输入标准图像格式的待处理高分辨率SAR图像。

步骤二：构建尺寸为(2d+1)×(2d+1)的滑动窗口对高分辨率SAR图像进行局部数据采集，计算有噪估计条件下的局部图像谱平稳特征，生成与图像等大小的有噪估计谱平稳特征图，其中：d为正整数，用来控制滑动窗口的尺寸。

特征图中数值表示局部图像的频谱平稳特征，该特征由频谱分布均匀性以及一致性共同表征，特征数值越大表明局部图像频域平稳性越差，越有可能是显著区域，数值越小则平稳性越高、越接近纯噪声，越可能是背景区域。

本步骤中，有噪估计条件下的局部图像谱平稳特征的计算方法为：

步骤二A：输入待处理高分辨率SAR图像，图像尺寸为M×N，其中：M与N表示图像的长宽，单位为像素数，使用(2d+1)×(2d+1)滑动窗口采集图像(x，y)处及其d邻域内像素构成局部像素矩阵I，其维度为(2d+1)×(2d+1)，矩阵元素为I_ij，i、j分别为矩阵行列坐标；

步骤二B：利用矩阵标准差估计全噪声假设下的噪声功率σ_N，并以σ_N为标准差生成一组维度为(2d+1)×(2d+1)高斯随机数，用来拟合功率为σ_N的噪声矩阵N，σ_N估计方法如下所示：

其中

式中，N_ij为噪声矩阵N中元素，i、j分别为矩阵行列坐标；

步骤二C：利用二维傅里叶变换(2D-FT)分别处理图像矩阵I与噪声矩阵N，获得对应的图像频谱矩阵F_I与噪声频谱矩阵F_N，计算公式如下所示：

式中，u、v为频域矩阵行列坐标；

步骤二D：计算图像频谱矩阵的分布方差σ² _F，计算方法如下所示：

其中

步骤二E：分别计算图像频谱矩阵F_I与噪声频谱矩阵F_N的分布熵ENT_F与ENT_N，计算方法如下所示：

式中，ε为限位参数，避免取对数过程中出现零值导致计算崩溃，ε取较小值以避免对结果造成干扰，本实施案例中令ε＝10^-5；

步骤二F：给定限位参数R，计算有噪估计条件下的局部图像谱平稳特征值SS(Spectral Stationary)，并将SS值赋给特征图中与图像坐标(x，y)相对应的位置，当滑动窗口取遍所有像素后生成与原图像尺寸等大小的有噪估计谱平稳特征图S，计算方法如下：

S＝{SS_xy|x＝1,2,…M,y＝1,2,…N}。

本发明给定一种基于K的限位参数R估计方式，但实际R的选择方法不必局限于此，估计方式如下所示：

R＝2×K×ln[2×(2d+1)+1]；

式中，K是调节系数，取值范围为正实数，无物理意义，用来调节限位参数R，(2d+1)是滑动窗口尺寸，exp()表示自然常数e为底的指数函数，ln()表示取自然对数。

R参数控制了图像频域分布熵对特征的决定性，即只有当图像频域熵小于该值时才进行平稳性计算，否则直接认为频谱为均匀分布，呈现白噪声特性，特征值赋为零；有噪估计谱平稳特征图S的一个数据分布示例如图3所示，该例子中特征图数据在10⁰至10⁶数量级间分布，且数据格式为浮点数，若进行8bit量化则数据对比度大幅压缩、细节丢失，若通过传统直方图统计需进行排序，消耗大量时间，且不同特征图数据分布区间较大，难以设置硬阈值进行数据筛选。图4与图5分别为对原始特征图进行线性与对数的动态范围压缩、8bit量化显示，可以看出，图像能粗略看出特征图反映出的信息，显著区域相对于背景呈现高强度值，但是大量细节损失、对比度被压缩。

步骤三：使用误差为e的数据格式无关p-tile分割法(Format-free e-accuracyp-tile segmentation,FF ep-tile)处理有噪估计谱平稳特征图，取分割结果前景部分为疑似显著区域。

本步骤中，误差为e的数据格式无关p-tile分割(Format-free e-accuracy p-tile segmentation,FF ep-tile)的处理方法为：

步骤三A：初始化迭代参数t＝1，输入待处理特征图S，其尺寸为M×N，将其按照元素下标序列化，得到一组等价的数值序列S，序列中元素表示为S_t，其中下标t表示该数值在序列中的序号，取值范围为1至M×N，按精度需求设置误差率e，给定经验性估计前景稀疏率p，令最小分割阈值T_min为特征图或等价序列S中最小值，最大分割阈值T_max为特征图或等价序列S中最大值，结果分割阈值T^*＝(T_max+T_min)/2，公式如下所示：

T_max＝max(S),T_min＝min(S),T^*＝(T_max+T_min)/2,t＝1；

式中，max与min分别为取最大与取最小函数；

本步骤中误差率e与前景稀疏率p为人工设置，其中参数e控制分割算法的精度，e越小则误差越小，分割精度越高，p为前景稀疏率，即图像中前景像素数量与全部像素数量的比值为p，p越大分割结果中前景比例越大，候选结果越多。

步骤三B：将序列值S_t赋值给临时分割阈值T，并进行判断：若T大于最小分割阈值T_min且小于最大分割阈值T_max，则进入步骤三C，否则令t＝t+1，重新执行本步骤(返回步骤三B)。

步骤三C：利用临时分割阈值T分割特征图或等价序列S，生成二值化的分割结果图或等价序列R，令分割前景为1、分割背景为0，并计算前景比例FR，公式化表述如下所示：

步骤三D：判断：|p-FR|是否小于等于给定误差率e，若是，则令结果分割阈值T^*等于当前临时分割阈值T，并进入步骤三G，否则进入步骤三E。

步骤三E：判断：前景稀疏率p是否大于前景比例FR，若是，则令最大分割阈值T_max等于当前临时分割阈值T，结果分割阈值T^*＝(T_max+T_min)/2，否则令最小分割阈值T_min等于当前临时分割阈值T，T^*＝(T_max+T_min)/2，公式化表述如下所示：

其中

步骤三F：判断：迭代参数t是否等于M×N，若是则进入步骤三G，否则令t＝t+1，返回步骤三B。

步骤三G：利用结果阈值T^*分割输入图像S，得到分割结果二值图，输出前景部分作为高分辨率SAR图像疑似显著区域。

步骤四：获取疑似显著区域的连通域结果，利用尺度滤波滤除像素数少于c的连通域图块，余下部分作为高分辨率SAR图像显著区域输出。本实施案例中c通过人工选择设置，但实际c的选择方法不必局限于此。

1、实验数据

本实验采用美国桑迪亚国家实验室(Sandia National Laboratories)提供的MiniSAR图像为实验数据，如图2所示，该图像为Ku波段、机场背景，图像分辨率为0.1米，数据格式为8bit量化Bitmap格式；显著性区域主要包含地面飞机目标、地面涂覆、地表设施、不同地物交界等；主要待滤除部分为跑道背景，如图2中a框所示，以及其他高散射复杂纹理背景，如b、c、d框所示。

2、实验内容和分析

本实验以图2所示MiniSAR图像为实验数据，以d＝5，即窗口尺寸为11构建滑动窗口进行有噪估计谱平稳特征提取。其中，取R＝3.2时，特征提取结果原始数据如图3所示，对应的线性及对数8bit量化图像如图4、图5所示。利用该特征图，取误差率e＝0.1％、前景稀疏率p＝5％得到的分割结果如图6所示。该分割步骤中，迭代次数t＝10，与图像总像素数比例为2.43×10^-6，分割阈值T^*＝696.534，实际稀疏率p^*＝2.59％，实际误差率e＝2.41％，该组结果说明本特征图中，存在大量相等数值，不存在使的误差率在2.41％以内的分割阈值，任何偏离696.534的分割阈值都将获得相等的误差率或者使误差率更大。基于该分割结果，取尺寸滤波阈值c＝50，获得显著区域检测结果如图7所示。调整R＝10，即降低图像分布熵的限位值，并取误差率e＝0.1％、前景稀疏率p＝10％，获得特征分割结果，如图8所示。该分割步骤，迭代次数12次，与总像素数比例为2.91×10^-6，分割阈值T^*＝0.344，实际稀疏率p^*＝10.04％，实际误差率e＝0.04113％。对比图8与图6分割结果，发现R与p的调整影响了分割结果中显著区域的候选结果数量，R越大、p越大则候选区域越多，噪点也越多。基于图8所示的分割结果，取尺寸滤波阈值c＝100，获得的显著区域检测结果如图9所示。图7、图9中结果中，包含地面目标、地面突起突变、地面涂覆等高价值显著区域均被检测出来，其中图7结果噪声点较少，图9噪声点较多但是区域一致性更好，二组结果均不同程度地为更进一步的处理提供了高质量参考数据，并滤除了背景冗余信息。

此外，针对本发明的方法步骤二与步骤三，步骤二中获得的特征图能够一定程度上看出显著区域，表明特征能够有效滤除高能背景、保留显著区域，实现设计需求。步骤三中FF ep-tile分割步骤，能够针对任意格式数据进行基于直方图统计的分割，且效率较高。两次实验中，分割算法迭代次数约为10次左右，相比于经典自适应阈值方法，例如大津阈值分割法(OTSU算法)，该方法要求输入数据为固定bit数量化，其中8bit量化时分割迭代次数为255，FF ep-tile分割方法在本实验中，在迭代次数上仅为OTSU算法的3.92％，即速度为OTSU算法25.5倍，大幅节约了分割时间，使得本基于有噪估计谱平稳特征的高分辨率SAR图像显著区域检测方法具有较高的执行效率。步骤二与步骤三使得本发明的方法作为高分辨率SAR图像粗处理中的显著区域检测方法，具有较好的性能与效率，具有较高的实用性。

Claims (4)

1.一种基于有噪谱平稳的SAR图像显著区域检测方法，其特征在于所述方法包括如下步骤：

步骤一：输入标准图像格式的高分辨率SAR图像；

步骤二：利用尺寸为(2d+1)×(2d+1)的滑动窗口对高分辨率SAR图像进行局部数据采集，计算有噪估计条件下局部图像谱平稳特征，最终生成与图像等大小的有噪估计谱平稳特征图，其中：d为正整数，用来控制滑动窗口尺寸；

步骤三：使用误差为e的数据格式无关p-tile分割法处理有噪估计谱平稳特征图，取分割前景部分为疑似显著区域；

步骤四：获取疑似显著区域的连通域结果，利用尺度滤波滤除像素数少于c的连通域图块，余下部分作为高分辨率SAR图像显著区域输出。

2.根据权利要求1所述的基于有噪谱平稳的SAR图像显著区域检测方法，其特征在于所述步骤二中，有噪估计条件下的局部图像谱平稳特征的计算方法为：

步骤二A：输入待处理高分辨率SAR图像，图像尺寸为M×N，其中：M与N表示图像的长宽，单位为像素数，使用(2d+1)×(2d+1)滑动窗口采集图像(x，y)处及其d邻域内像素构成局部像素矩阵I，其维度为(2d+1)×(2d+1)，矩阵元素为I_ij，i、j分别为矩阵行列坐标；

步骤二B：利用矩阵标准差估计全噪声假设下的噪声功率σ_N，并以σ_N为标准差生成一组维度为(2d+1)×(2d+1)高斯随机数，用来拟合功率为σ_N的噪声矩阵N；σ_N估计方法如下所示：

其中

式中，N_ij为噪声矩阵N中元素，i、j分别为矩阵行列坐标；

步骤二C：利用二维傅里叶变换分别处理图像矩阵I与噪声矩阵N，获得对应的图像频谱矩阵F_I与噪声频谱矩阵F_N，计算公式如下所示：

式中，u、v为频域矩阵行列坐标；

步骤二D：计算图像频谱矩阵的分布方差σ² _F，计算方法如下所示：

其中

步骤二E：分别计算图像频谱矩阵F_I与噪声频谱矩阵F_N的分布熵ENT_F与ENT_N，计算方法如下所示：

式中，ε为限位参数；

步骤二F：给定限位参数R，计算有噪估计条件下的局部图像谱平稳特征值SS，并将SS值赋给特征图中与图像坐标(x，y)相对应的位置，当滑动窗口取遍所有像素后生成与原图像尺寸等大小的有噪估计谱平稳特征图S，计算方法如下：

S＝{SS_xy|x＝1,2,…M,y＝1,2,…N}。

3.根据权利要求2所述的基于有噪谱平稳的SAR图像显著区域检测方法，其特征在于所述R＝2×K×ln[2×(2d+1)+1]，式中，K是调节系数，(2d+1)是滑动窗口尺寸。

4.根据权利要求1所述的基于有噪谱平稳的SAR图像显著区域检测方法，其特征在于所述步骤三中，误差为e的数据格式无关p-tile分割的处理方法为：

步骤三A：初始化迭代参数t＝1，输入待处理特征图S，其尺寸为M×N，将其按照元素下标序列化，得到一组等价的数值序列S，序列中元素表示为S_t，其中下标t表示该数值在序列中的序号，取值范围为1至M×N，按精度需求设置误差率e，给定经验性估计前景稀疏率p，令最小分割阈值T_min为特征图或等价序列S中最小值，最大分割阈值T_max为特征图或等价序列S中最大值，结果分割阈值T^*＝(T_max+T_min)/2，公式如下所示：

T_max＝max(S),T_min＝min(S),T^*＝(T_max+T_min)/2,t＝1；

式中，max与min分别为取最大与取最小函数；

步骤三B：将序列值S_t赋值给临时分割阈值T，并进行判断：若T大于最小分割阈值T_min且小于最大分割阈值T_max，则进入步骤三C，否则令t＝t+1，返回步骤三B；

步骤三C：利用临时分割阈值T分割特征图或等价序列S，生成二值化的分割结果图或等价序列R，令分割前景为1、分割背景为0，并计算前景比例FR，公式化表述如下所示：

步骤三D：判断：|p-FR|是否小于等于给定误差率e，若是，则令结果分割阈值T^*等于当前临时分割阈值T，并进入步骤三G，否则进入步骤三E；

步骤三E：判断：前景稀疏率p是否大于前景比例FR，若是，则令最大分割阈值T_max等于当前临时分割阈值T，结果分割阈值T^*＝(T_max+T_min)/2，否则令最小分割阈值T_min等于当前临时分割阈值T，T^*＝(T_max+T_min)/2，公式化表述如下所示：

其中

步骤三F：判断：迭代参数t是否等于M×N，若是则进入步骤三G，否则令t＝t+1，返回步骤三B；

步骤三G：利用结果阈值T^*分割输入图像S，得到分割结果二值图，输出前景部分作为高分辨率SAR图像疑似显著区域。

Images