基于递归OTSU算法的SAR图像水域自动检测方法与流程

本发明涉及合成孔径雷达探测技术领域，特别涉及一种基于递归otsu算法的sar图像水域自动检测方法，即基于递归最大类间方差(otsu)算法的sar图像水域自动检测方法，适用于自适应提取sar图像中的水域部分，进而完成对水体动态变化的实时监测。

背景技术：

合成孔径雷达(syntheticapertureradar,sar)因全天候、全天时的监测能力，得到了广泛的应用。至今，世界各国已成功发射多颗先进的星载sra卫星，如alos-2、radarsat-2,、cosmo-skymed、terrasar-x、sentinel-1a/b以及高分三号等，它们为合成孔径雷达各项应用提供了有效的数据支持。

sar图像中含有水域、植被区、城区等区域，在对sar图像分类的过程中，可以对水域进行单独检测，而且水域检测结果作为数字高程模型(dem)编辑过程中的补充信息，也可用于监测自然灾害(洪水等)以及水资源(如湖泊，河流等)的分布变化。

t.hahmann等人在文章“automaticextractionofwaterbodiesfromterrasar-xdata”中使用基于阈值的检测方法进行水域提取，即将sar图像中的后向散射系数值小于或等于给定阈值的像素标记为水域；基于阈值的检测方法具有较快的检测速度，因而非常适合于大规模数据的批量处理；但是，基于阈值的检测方法仍然存在以下问题：在一般情况下，由于镜面散射的原因，水域的后向散射系数较低；但是，当水体受环境(如风速、降水等)影响而表面有起伏时，水域的后向散射系数也会有所提高，这意味着，当利用不同时刻录取的sar图像数据进行水域检测时，不同的sar图像的水域检测阈值需要单独确定；当批量处理不同的sar图像数据时，水域检测阈值的确定将耗费大量的时间，严重降低批量处理多幅sar图像数据的效率。

针对以上问题，a.wendleder等人在文章“tandem-xwaterindicationmask:generationandfirstevaluationresults”中提出通过设置高低两个阈值的方式解决不同环境情况下的水域检测问题，最后将高低两个阈值各自得到的水域检测结果进行有效融合以获得最终的结果。但是，文中高低两个阈值的获取，是在分析全球大约1700幅sar图像的基础上，通过统计方法得到的，缺乏客观性。

目前，针对sar图像水域检测中批量高效的处理需求，基于阈值的检测方法不失为一种好的选择，但却无法实现sar图像水域检测的高效和自动处理。

技术实现要素：

针对上述现有技术存在的不足，本发明提出了一种基于递归otsu算法的sar图像水域自动检测方法，即基于递归最大类间方差(otsu)算法的sar图像水域自动检测方法，该种基于递归otsu算法的sar图像水域自动检测方法对sar图像的功率信息递归使用otsu算法，能够自适应得到sar图像对应的阈值。

为达到上述目的，本发明的实施例采用如下技术方案予以实现。

一种基于递归otsu算法的sar图像水域自动检测方法，包括以下步骤：

步骤1，获取sar图像功率数据，并对所述sar图像功率数据进行相干斑滤波，得到相干斑滤波后的sar图像数据；所述sar图像功率数据和相干斑滤波后的sar图像数据分别为m×p维矩阵，m、p分别为大于0的正整数；

步骤2，对相干斑滤波后的sar图像数据进行归一化，得到归一化后的sar图像数据；

步骤3，对归一化后的sar图像数据进行灰度化，得到sar灰度图像；

步骤4，对sar灰度图像进行阈值估计，进而确定sar灰度图像中的最优阈值t；

步骤5，构造与sar图像功率数据同样大小的矩阵，记为标记矩阵；然后根据sar灰度图像中的最优阈值t，得到sar图像功率数据对应的二值化标记矩阵；最后对sar图像功率数据对应的二值化标记矩阵进行形态学处理，进而得到sar图像功率数据的水域检测标记图。

本发明与现有的技术相比具有以下优点：本发明对于大多数sar图像进行水域检测时，都能自动的选取阈值，而且检测结果的正确性较高，具有很好的自动性和鲁棒性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的一种基于递归otsu算法的sar图像水域自动检测方法流程图；

图2为terrasar-x录取的西安某区域(覆盖沣河区域)的未经相干斑滤波的sar图像功率图；

图3为对图2依次进行相干斑滤波和灰度化后得到的结果示意图；

图4为本发明实施例进行初步水域检测的结果图；

图5为形态学处理后的水域检测结果图；

图6为西安地区sar图像中对应的水域区域参考图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1，为本发明的一种基于递归otsu算法的sar图像水域自动检测方法流程图；其中所述基于递归otsu算法的sar图像水域自动检测方法，包括以下步骤：

步骤1，获取sar图像功率数据，并对所述sar图像功率数据进行相干斑滤波，得到相干斑滤波后的sar图像数据；本实施例采用lee滤波方法，选取n×n的滤波窗口对所述sar图像数据进行相干斑滤波，n为大于0的奇数，本实施例中n取值为7，所述sar图像功率数据和相干斑滤波后的sar图像数据分别为m×p维矩阵，m、p分别为大于0的正整数。

步骤2，sar图像数据归一化。若是直接利用相干斑滤波后的sar图像数据中所有像素点值中的最大值对该相干斑滤波后的sar图像数据进行归一化，则可能会造成相干斑滤波后的sar图像数据中的水域区域与背景之间的对比度严重下降，从而影响后续的算法性能；为了避免这种情况，提高相干斑滤波后的sar图像数据中水域区域与背景之间的对比度，本发明对相干斑滤波后的sar图像数据进行归一化，其子步骤如下：

2.1计算相干斑滤波后的sar图像数据的第90百分位数，即将相干斑滤波后的sar图像数据按照像素点值进行从小到大排序，得到从小到大排序后的m×p个数据，并从小到大排序后的m×p个数据中选取位置处的像素点值，记为参考像素点值z90；表示向上取整。

2.2将相干斑滤波后的sar图像数据中所有大于参考像素点值z90的像素点值分别用z90代替，其余小于或等于参考像素点值z90的像素点值保持不变，进而得到修改后的sar图像数据。

2.3对修改后的sar图像数据进行归一化，即将修改后的sar图像数据中所有像素点的值分别除以z90，进而得到归一化后的sar图像数据，所述归一化后的sar图像数据中每一个像素点的值为0到1之间的数据。

通过上述步骤，本发明可完成相干斑滤波后的sar图像数据的归一化。

步骤3，相干斑滤波后的sar图像数据的灰度化，所述相干斑滤波后的sar图像数据的灰度化是将归一化后的sar图像数据中每一个像素点的值分别从0到1之间扩展为0到l之间，进而得到sar灰度图像，所述sar灰度图像的最大灰度级为l，l为大于0的正整数，其子步骤为：

3.1将归一化后的sar图像数据中所有像素点的值分别乘以l，得到扩展后的sar图像数据，所述扩展后的sar图像数据中每个像素点的值域范围分别在0和l之间，l表示sar灰度图像中的最大灰度级，本实施例中l取值为255。

3.2对扩展后的sar图像数据中所有像素点的值分别进行取整操作，此处取整操作为四舍五入操作，使得扩展后的sar图像数据中所有像素点的值均为整数；进而得到sar灰度图像，所述sar灰度图像的最大灰度级为l。

通过上述步骤，本方案可完成归一化后的sar图像数据的灰度化，从而得到sar灰度图像，所述sar灰度图像的最大灰度级为255。

下面结合实测数据实验对本发明的效果做进一步的说明：图2为terrasar-x录取的西安某区域(覆盖沣河区域)的未经相干斑滤波的sar图像功率图，图3为对图2依次进行相干斑滤波和灰度化后得到的结果示意图；将图2和图3进行对比可以发现：相干斑滤波后的sar图像数据中的噪声得到了较好的去除；在不影响水域区域与背景的对比度的情况下，对sar图像数据依次进行了归一化和灰度化。

步骤4，利用递归的otsu阈值估计算法对sar灰度图像进行阈值估计；当sar灰度图像中的水域区域面积大于sar灰度图像面积的30％时，标准的otsu阈值估计算法具有良好的阈值分割性能，但是当sar灰度图像中的水域区域面积小于sar灰度图像面积的10％时，标准的otsu阈值估计算法性能迅速下降。一般来说，sar图像功率数据中的低功率区域面积总是远远小于较高功率区域面积，因而在这种情况下，利用标准的otsu阈值估计算法很难得到最优阈值。

为了克服上述情况，本发明实施例引入递归的otsu阈值估计算法来确定最优阈值；所述递归的otsu阈值估计算法的关键思想为：根据标准的otsu阈值估计算法得到阈值，移除sar灰度图像中灰度值高于阈值的所有像素点，然后利用sar灰度图像中灰度值小于或等于阈值的剩余像素点，以及标准的otsu阈值估计算法重新估计阈值；不断重复进行上述过程直至达到终止条件，进而得到sar灰度图像中的最优阈值。

4.1初始化：令目标区域表示sar灰度图像中的水域区域，背景区域表示sar灰度图像中的非水域区域；令otsu(a,b)表示利用sar灰度图像中灰度值范围在a和b之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值，本实施例中a取值为0，b的初值为l；分别令k表示第k次递归，k的初始值为1，令tk表示第k次递归后利用sar灰度图像中灰度值范围在0和tk-1之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值，且t0＝l，l表示sar灰度图像中的最大灰度级，本实施例中l取值为255；令t1表示第1次递归后利用sar灰度图像中灰度值范围在0和l之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值，t1和tk的定义表达分别为：

t1＝otsu(0,255)

tk＝otsu(0,tk-1)

其中，tk-1表示第k-1次递归后利用sar灰度图像中灰度值范围在0和tk-2之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值。

4.2计算第k次递归后利用sar灰度图像中灰度值在0和tk-1之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值tk，其过程如下：

4.2.1计算灰度值i在sar灰度图像中灰度值范围为{0,...,tk-1}的像素点中出现的概率pi，i∈{0,...,tk-1}。

4.2.2当令tk表示第k次递归后利用sar灰度图像中灰度值范围在0和tk-1之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值时，阈值的选取范围为{0,...,tk-1}；令待估计的阈值为t，并分别将sar灰度图像中灰度值小于或等于t的所有像素点作为目标区域中的像素点，将sar灰度图像中灰度值大于t的所有像素点作为背景区域中的像素点；t的初始值为0。

4.2.3分别计算待估计的阈值为t时目标区域在sar灰度图像中灰度范围为{0,...,tk-1}的像素点中出现的概率ω0(t)、待估计的阈值为t时背景区域在sar灰度图像中灰度范围为{0,...,tk-1}的像素点中出现的概率ω1(t)、待估计的阈值为t时目标区域的平均灰度值μ0(t)和待估计的阈值为t时背景区域的平均灰度值μ1(t)，其计算公式分别为：

其中，i'∈{0,...,t}，pi'表示灰度值i'在sar灰度图像中灰度值范围为{0,...,tk-1}的像素点中出现的概率，i″∈{t+1,...,tk-1}，pi″表示灰度值i″在sar灰度图像中灰度值范围为{0,...,tk-1}的像素点中出现的概率。

进而计算待估计的阈值为t时sar灰度图像中目标区域与背景区域间的类间方差

4.2.4令t分别取0至tk-1，重复执行4.2.3，进而分别得到待估计的阈值为0时sar灰度图像中目标与背景间的类间方差至待估计的阈值为tk-1时sar灰度图像中目标区域与背景区域间的类间方差记为sar灰度图像中的tk-1+1个类间方差

将sar灰度图像中的tk-1+1个类间方差取得最大值时对应的待估计的阈值，作为第k次递归后利用sar灰度图像中灰度值范围在0和tk-1之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值tk，其计算表达式为：

其中，arg表示取最大值时对应的参数值。

4.3随着递归次数的增加，与sar灰度图像中剩余像素里的背景区域面积相比，目标区域的面积不断的增加，因而sar灰度图像中的最优阈值从集合中选取，下标为递归终止时对应递归次数；递归的otsu阈值估计算法主要面临两个问题：一个问题是何时停止递归过程；一般而言，当两次连续的递归过程之间的阈值变化小于设定的阈值d时，即如果|tk-tk-1|≥d，则令k加1，返回4.2。

当两次连续的递归过程之间的阈值变化小于设定的阈值d时，即如果|tk-tk-1|<d，则终止递归，并将递归停止时对应的第k次递归后利用sar灰度图像中灰度值范围在0和tk-1之间的所有像素点以及标准的otsu阈值估计算法进行阈值估计得到的阈值，记为递归终止阈值下标为递归终止时对应递归次数，此时意味着递归过程已十分的稳定，且递归终止时对应得到个阈值，记为分别为递归终止前对应得到的个阈值，tn表示递归停止前得到的第n个阈值；本实施例中d＝3；第二个问题是如何从集合中确定sar灰度图像中的最优阈值。

4.4为了解决第二个问题，引入归一化的类间方差，使得每次递归过程中的阈值分割性能具有可比性；令η(tn)表示递归停止前得到第n个阈值tn时对应的归一化类间方差，其表达式为：

其中，当阈值为tn时，在sar灰度图像中，只考虑像素点的灰度值在0和tn-1之间的所有像素点，将0和tn-1之间的像素点灰度值小于或等于tn的所有像素点设为目标区域，其余为背景区域；表示阈值为tn时sar灰度图像中目标区域与背景区域间的类间方差，表示sar灰度图像中灰度值范围在0和tn-1之间的所有像素点的总方差，其有如下定义

其中，表示灰度值在sar灰度图像中灰度值范围为{0,...,tn}的像素点中出现的概率，μt表示sar图像灰度图中灰度值在0和tn-1之间的所有像素点的平均灰度值。

4.5令n分别取1至重复执行4.4，进而分别得到递归停止前得到第1个阈值t1时对应的归一化类间方差η(t1)至递归终止时得到第个阈值时对应的归一化类间方差记为个归一化类间方差η，

需要注意的是，上式中灰度值在sar灰度图像中灰度值范围为{0,...,tn-1}的像素点中出现的概率在每次递归过程中已经被更新。一般而言，当目标区域(水域)和背景区域(除水域外的其他区域)面积相当时，标准的otsu阈值估计算法具有良好的性能，且归一化类间方差可达到最大值。

4.6由于在递归的过程中sar灰度图像中目标区域与背景区域间的差异越来越小，甚至到最后仅剩下目标区域，因此归一化类间方差持续下降；因此根据所述个归一化类间方差η，利用下述准则确定sar灰度图像中的最优阈值t：

其中，argmax{η}表示η取最大值时对应的参数值。

通过上述操作发现，利用递归的otsu阈值估计算法完成阈值的精确估计，其在处理小面积目标图像时仍然具有良好的估计性能。

步骤5，构造与sar图像功率数据同样大小的矩阵，记为标记矩阵；然后根据sar灰度图像中的最优阈值t对sar图像功率数据进行水域区域检测，将小于或等于最优阈值t的所有像素点分别标记为水域区域，并分别在对应标记矩阵中的相应位置处标记为1；将大于最优阈值t的剩余像素点分别标记为背景区域，并分别在对应标记矩阵中的相应位置处标记为0；进而得到sar图像功率数据对应的二值化标记矩阵；最后对sar图像功率数据对应的二值化标记矩阵进行形态学处理，得到sar图像功率数据的水域检测标记图。

具体地，得到sar图像功率数据的水域检测标记图即水域检测完成，然后计算水域检测完成后其检测结果的完整性和正确性。

参照图4和图5，图4为本发明实施例进行初步水域检测的结果图；图5为形态学处理后的水域检测结果图；分别令a1表示使用本发明方法检测到的水域区域面积与事先设定的参考水域区域面积的交集,令a2表示事先设定的参考水域区域面积,令a3表示使用本发明方法检测到的水域区域面积，本实施例中事先设定的参考水域区域面积如图6所示，图6为西安地区sar图像中对应的水域区域参考图；令pcom为水域区域检测结果的完整性参数，pcor为水域区域检测结果的正确性参数，其表达式分别为：

计算本发明在实施例中水域区域检测结果的正确性参数为85.5％，水域区域检测结果的完整性参数为77.1％，实验结果显示本发明方法在水域检测时有较好的检测效果。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。