一种超像素的概率因子TMF的SAR图像海岸线检测算法的制作方法

本发明涉及一种超像素的概率因子tmf的sar图像海岸线检测算法，属于海岸线检测领域。

背景技术：

海岸线检测是海域使用和管理中的一个重要的环节，也是国家管理海洋和陆地资源的一个重要的手段，因此研究海岸线的变化有着重要的意义，导致海岸线变化的原因主要由以下几方面：(1)地壳的运动，这是海岸线变化的最主要原因。(2)冰川的融化和扩展。合成孔径雷达(syntheticapertureradar,sar)是一种主动式微波探测器，利用合成孔径原理，信号处理方法和脉冲压缩技术，通过较小尺寸的真实天线孔径合成较大的等效天线孔径来成像。相对于光学成像，合成孔径雷达具有全天候成像等优势，因而使得sar图像在战略目标识别与探测，灾害控制，国土资源监测，海域使用管理，地图测绘，船舰目标识别等领域起到了重要的作用，虽然sar图像虽然sar图像有全天候的特性，但是sar图像的解析是较为困难的，从海岸线检测角度来看，其原因是由于sar图像有很强的相干斑，使得较精确提取海岸线成为较困难的问题，因此研究sar图像的海岸线检测具有重要的意义。

目前为止，常见的用于海岸线检测的技术有：基于邻域相似性局部窗的超像素算法和三重马尔科夫随机场(tmf)算法。但是传统的超像素算法在计算像素特征时采用的是恒定大小的矩形窗作为邻域窗，将邻域窗内所有像素特征的均值作为中心点的特征，当中心点位于边缘处时，此时矩形窗内含有边缘，就会导致计算出的特征很模糊，进而使边缘处超像素的边缘贴合度较低。由于在已有基于超像素算法文献中，未发现基于超像素的区域合并海岸线检测算法。而已有马尔科夫随机场算法具有贴合度不佳的特点，在海岸线检测时特别是低对比度边界容易脱离海岸边缘。但是现有的算法尚不能有效的解决这一问题。

技术实现要素：

本发明针对以上问题的提出，提出一种超像素的概率因子tmf的sar图像海岸线检测算法，其特征在于包括如下步骤：

s1：读取合成孔径雷达采集到的sar图像，生成对应的纹理图像itexture；

s2：根据图像的大小输入种子点数k，其取值范围为所述纹理图像itexture大小的10％-20％，并在所述纹理图像中均匀选取种子点；

s3：遍历所述的纹理图像itexture，每一个点的7像素×7像素的局部窗内计算局部纹理图像ipart，判断当前局部窗内是否含有边缘并根据所述局部窗内是否含有边缘确定用于计算特征的邻域点的集合c，将无边缘的局部窗作为邻域点的集合c中的特征点；

s4：根据所述步骤s3所得的邻域点的集合c计算所述步骤s2中所选取的种子点所有点的特征均值和方差；

s5：计算每个邻域点和选取的种子点之间的特征距离dci,j，邻域点和种子点之间位置的欧氏距离dsi,j及种子点与邻域点空间距离和特征距离的加权距离di,j，通过比较di,j的大小将当前所述di,j，将di,j较小的邻域点j归入较小的种子点中形成超像素；根据均值法更新每个选取的种子点特征；

s6：重复步骤s3-s5直至前后两次超像素的边缘不再发生变化，超像素恒定后进行步骤s7；

s7：初始化超像素标签，同时初始化辅助场u；遍历所述合成孔径雷达采集到的sar图像计算出超像素s属于第i类的概率ps,i初始化二维向量的标记场x并计算第i类超像素的均值μi和第i类超像素方差σi；

s8：根据步骤s7计算出的第i类超像素的均值μi计算gamma分布

s9：根据计算的二维向量的标记场x，辅助场u计算势能函数，并更新超像素和辅助场u；

s10：根据均值法更新每一类标签的特征和每一个超像素的特征，并计算出ps,i，进而更新二维向量的标记场x；

s11：重复步骤s7-s10，直至标记场不变时停止更新；更新后使用icm算法和sg算法对参数集θ＝{αx,αua,αub}更新，最后通过icm算法进行条件迭代，得到海岸的检测结果。

进一步的，步骤s3中：当局部窗不含有边缘时，则使用全部的像素值计算特征；当局部窗含有边缘时，则使用和中心点属于同一类的像素集合计算特征。

进一步的，每一个像素点的统计量st形成的图像即所述纹理图像itexture的计算过程包括：

其中，μi表示中心点i局部窗内的均值，ij表示中心点i局部窗内的第j像素的像素值，ni表示局部窗内的像素的集合，m表示ni中像素的个数。

进一步的，所述邻域点和种子点之间的距离dc1,2：

其中，μ1表示种子点的统计量st的均值，μ2表示邻域点的统计量st的均值，var1表示种子点的统计量st的方差，var2表示邻域点的统计量st的方差，wvar表示均值项的权值，w2表示方差项的权值；

所述均值项的权值wvar：

其中，邻域点和中心点统计量方差的相似程度varmin，μmin表示中心点局部窗内统计量均值和邻域点局部窗内统计量均值的比值的较小值；

所述方差项的权值w2：

进一步的，所述邻域点和种子点之间位置的欧氏距离dsi,j：

所述邻域点和种子点空间距离和特征距离的加权距离di,j：

其中，m表示dci,j和dsi,j之间的权重系数，s表示第一次选取种子点时，两个种子点之间的距离。

进一步的，二维向量的标记场x每一维表示当前超像素属于两类之一的概率为：

其中，i表示超像素的标签，s表示超像素的位置；

所述第s位置的超像素和第i类的相似程度第s位置的超像素的特征和第i类特征之间的距离ds,i为：

其中，i表示超像素的标签，isp表示当前超像素的均值，σsp表示当前超像素的方差，μi表示第i类超像素的均值，σi表示第i类超像素方差。

进一步的，将ps,ii∈{1,2}形成一个向量作为标记场x的实现，计算μi和σi：

标记场的后验概率服从均值为μi，视数为l的gamma分布为：

其中，i表示标签，i的取值范围为i∈{1,2}，l表示sar图像的视数，isp表示超像素内像素均值，s表示的是整个图像中所有超像素位置的集合。

进一步的，基于标记场x的势能函数模型为：

其中，ps,1表示第s位置超像素属于第一类的概率，ps,2表示第s位置超像素属于第二类的概率，t是s的邻域超像素，pt,1表示超像素t属于第一类的概率，pt,2表示超像素t属于第二类的概率，{αx,αua,αub}表示参数集；

wedge(xs,xt)计算如下：

将初始的超像素图像通过初始聚类将图像中的超像素分成两类即标签值是1和2两类当超像素s和t都属于第一类即1时，

进一步的，当超像素s和超像素t的特征均介于第一类和第二类中间且偏向于第一类略大于第二类时，根据辅助场u的值判断其种类；当超像素s和超像素t中有一个与第一/二类相似，所述与第一类相似就是指ps,1/pt,1远大于ps,2/pt,2，与第二类非常相似就是指ps,2/pt,2远大于ps,1/pt,1，另一个与该类比较相似时，根据标签值和辅助场u的值联合判断其种类。

本发明的优点在于：本专利不再采用传统的增加特征维数的方法，而是由液体扩散系数的启发提出了用统计量表示图像的方法，提出这种纹理，使用该统计量可以在一定程度上降低sar图像中海洋和陆地的纹理复杂性，同时在一定程度上增强海洋和陆地的对比度。同时，在本专利模型中，标记场不再仅仅是一个标签值而是一个概率因子本专利模型中将标记场变成一个二维向量，能够有效的解决贴合度不佳的特点，在海岸线检测时特别是低对比度边界容易脱离海岸边缘的问题。

附图说明

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

图1为本发明的整体流程图。

图2(a)-(d)为本发明得到的海岸线检测图。

具体实施方式

为使本发明的实施例的目的、技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1所示，一种超像素的概率因子tmf的sar图像海岸线检测算法，作为优选的实施方式，包括如下步骤：

s1：读取合成孔径雷达采集到的sar图像，生成对应的纹理图像itexture；

s2：根据图像的大小输入种子点数k，其取值范围为纹理图像itexture大小的10％-20％，并在纹理图像itexture中均匀选取种子点；

s3：遍历的纹理图像itexture，在选取的每一个点的7像素×7像素的局部窗内计算局部纹理图像ipart，判断当前局部窗内是否含有边缘并根据局部窗内是否含有边缘确定用于计算特征的邻域点的集合c；

s4：根据步骤s3所得的邻域点的集合c计算步骤s2中所选取的种子点所有点的特征均值和方差；

s5：计算每个邻域点和选取的种子点之间的距离dci,j，邻域点和种子点之间位置的欧氏距离dsi,j及di,j，通过比较di,j的大小将当前di,j，将di,j较小的邻域点j归入较小的种子点中形成超像素；根据均值法更新每个选取的种子点特征；

s6：重复步骤s3-s5直至前后两次超像素的边缘不再发生变化，超像素恒定后进行步骤s7；

s7：初始化超像素标签，同时初始化辅助场u；遍历原图像计算出ps,i，初始化二维向量的标记场x并计算第i类超像素的均值μi和第i类超像素方差σi；

s8：根据步骤s7计算出的第i类超像素的均值μi计算gamma分布

s9：根据计算的二维向量的标记场x，辅助场u计算势能函数，并对超像素和辅助场u进行更新；

s10：根据均值法更新每一类标签的特征和每一个超像素的特征，并计算出ps,i，进而更新二维向量的标记场x；

s11：重复步骤s7-s10，直至标记场不变时停止更新；更新后使用icm算法和sg算法对参数集θ＝{αx,αua,αub}更新，最后通过icm算法进行条件迭代，得到海岸的检测结果。

在本实施方式中，步骤s3中：当局部窗不含有边缘则使用全部的像素值计算特征；当局部窗含有边缘时，则使用和中心点属于同一类的像素集合计算特征。

在本实施方式中，每一个像素点的统计量st形成的图像即纹理图像itexture的计算过程包括：

其中，μi表示中心点i局部窗内的均值，ij表示中心点i局部窗内的第j像素的像素值，ni表示局部窗内的像素的集合，m表示ni中像素的个数。可以理解为在其它实施方式中，纹理图像itexture的可以用其它的物理量代替或表示，只要能够满足能够清楚的区别开所需的纹理图像并且能够有效的表示图像即可。

在本实施方式中，邻域点和种子点之间的距离dc1,2：

其中，μ1表示种子点的统计量st的均值，μ2表示邻域点的统计量st的均值，var1表示种子点的统计量st的方差，var2表示邻域点的统计量st的方差，wvar表示均值项的权值，w2表示方差项的权值；

均值项的权值wvar：

其中，邻域点和中心点统计量方差的相似程度varmin，μmin表示种子点的统计量均值和邻域点的统计量均值的比值的较小值；方差项的权值w2：

在本实施方式中对于dsi,j的计算公式如下：

di,j的计算方式为：

二维向量的标记场x每一维表示当前超像素属于两类之一的概率为：

其中，i表示超像素的标签，s表示超像素的位置；第s位置的超像素和第i类的相似程度。

ds,i的计算方式为：

其中，i表示超像素的标签，isp表示当前超像素的均值，σsp表示当前超像素的方差，μi表示第i类超像素的均值，σi表示第i类超像素方差。

作为优选的实施方式，将ps,ii∈{1,2}形成一个向量作为标记场x的实现，计算μi和σi：

假设标记场的后验概率服从均值为μi，视数为l的gamma分布为：

其中，i表示标签，i的取值范围为i∈{1,2}，l表示sar图像的视数，isp表示超像素内像素均值。

在本实施方式中，基于标记场x的势能函数模型为：

其中，ps,1表示第s位置超像素属于第一类的概率，ps,2表示第s位置超像素属于第二类的概率，t是s的邻域超像素，pt,1表示超像素t属于第一类的概率，pt,2表示超像素t属于第二类的概率，{αx,αua,αub}表示参数集；

在本实施方式中，wedge(xs,xt)计算如下：

将初始的超像素图像通过初始聚类将图像中的超像素分成两类即标签值是1和2两类当超像素s和t都属于第一类即1时，

在本实施方式中，当超像素s和超像素t的特征均介于第一类和第二类中间且偏向于第一类略大于第二类时，根据辅助场u的值判断其种类；当超像素s和超像素t中有一个与第一/二类相似，与第一类相似就是指ps,1/pt,1远大于ps,2/pt,2，与第二类非常相似就是指ps,2/pt,2远大于ps,1/pt,1，另一个与该类比较相似时，根据标签值和辅助场u的值联合判断其种类。

实施例一：envisat图像检测结果对比

取种子数k分别为：150，150，250，250。每个邻域点和选取的种子点之间的距离dc与邻域点和种子点之间位置的欧氏距离ds之间的权重系数m为0.5，局部窗大小是5×5，最大迭代次数为5次。对于tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。

对比试验gmrf的参数设置如下：条件迭代模式的迭代次数为5，滤波窗的大小是3×3。

对比试验基于超像素的tmf的参数设置如下：种子数k取值为150，150，250，250，dc与ds之间的权重系数m为0.5，局部窗大小是3×3，最大迭代次数为5次，tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。得到图2(a)-(d)所示为本具体实施方式得到的海岸线检测图。

表1针对envisat图像的三种算法rmse对比

实施例二：terra图像检测结果对比与分析

种子数k分别是150，150，150，150。dc与ds之间的权重系数m为0.5，局部窗大小是5×5，最大迭代次数为5次。对于tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。对比算法gmrf的参数设置如下：条件迭代模式的迭代次数为5，滤波窗的大小是3×3。对比算法基于超像素的tmf的参数设置如下：种子数k取值为150，150，150，150，dc与ds之间权重系数m为0.5，局部窗大小是3×3，最大迭代次数为5次。对于tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。

表2针对terra图像的三种算法rmse对比

实施例三：radasat图像检测结果对比与分析

种子数k分别是250，250，250，250。dc与ds之间的权重系数m为0.5，局部窗大小是5×5，最大迭代次数为5次。对于tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。对比算法gmrf的参数设置如下：条件迭代模式的迭代次数为5，滤波窗的大小是3×3。对比算法基于超像素的tmf的参数设置如下：种子数k取值为250，250，250，250，dc与ds之间的权重系数m为0.5，局部窗大小是3×3，最大迭代次数为5次。对于tmf部分αx的初始值为1，αua的值为0.6，αub的值为0.6，外迭代的迭代次数是10，内迭代的迭代次数是10，用于判断不再发生变化的门限值是0.0001。

以上，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。