一种逆高斯分布主动轮廓模型海岸线检测的计算方法与流程

本发明涉及一种图像处理技术，特别是一种主动轮廓模型海岸线检测的计算方法。

背景技术：

目前在海岸线检测问题逐渐得到了大家的关注，在海岸线的检测计算方法上，通常采用gamma分布的sar图像模型，由于基于gamma分布的sar图像模型是假设雷达散射面近似均匀，斑点噪声服从单位均值的gamma分布，同理，假设雷达散射面为均匀，斑点噪声分量服从高斯分布、负指数分布和rayleigh分布分别组合形成不同的sar图像统计模型且假设乘性模型为

i(x,y)＝f(x,y)n(x,y)

其中，i(x,y)表示最终雷达实测得到sar图像在点(x,y)的像素值，f(x,y)为点(x,y)后向散射系数，n(x,y)表示独立于f(x,y)均值为1、方差为的随机噪声，设ui(x,y)表示i(x,y)的均值，则由乘性模型得

ui(x,y)＝uf(x,y)

当单位空间散射系数为常量时则

对于海陆分布比较均匀且海陆分界明显的radarsat卫星合成的sar图像，用基于gamma分布的区域主动轮廓模型检测海岸线可以得到理想的结果，但是对于envisat卫星合成的sar图像，该方法不仅处理速度慢，需要较多的迭代次数；而且对于存在弱边缘的sar图像，利用基于gamma分布的区域主动轮廓模型检测海岸线不能得到理想的结果。

由上分析gamma分布中的参数中图像视数无法影响水平集演化结果，而影响演化的关键只是区域均值，在envisat图像中区域均值差异较小时，以gamma分布为基础的水平集演化方程无法得到较好的结果。如果概率模型单一参数无法得到有效检测结果，是否存在多参数，并且具有近似gamma分布非对称特性的概率模型可以解决这一问题？

技术实现要素：

为解决现有技术存在的上述问题，本发明要设计一种不仅处理速度快、迭代次数较少的逆高斯分布主动轮廓模型海岸线检测的计算方法，以解决存在弱边缘的sar图像传统方法识别效果不好的问题。

为了实现上述目的，本发明的技术方案如下：一种逆高斯分布主动轮廓模型海岸线检测的计算方法，包括以下步骤：

a、初始化水平集函数φ(x,y)；

基于符号距离函数初始化水平集函数φ(x,y)并记为φ⁰(x,y)；

其中，x,y表示图像中的像素位置坐标；

b、计算目标区域的均值u1(x,y)和背景区域的均值u2(x,y)；

读取sar图像像素i(x,y)，按以下公式计算u1(x,y)和u2(x,y)：

设逆高斯分布概率密度函数如下：

其中，λ0(x,y)表示均值的三次方与方差的比值，u0(x,y)表示目标区域和背景区域的均值，当目标区域和背景区域均值相近时通过参数λ0(x,y)区别目标和背景像素点。

假设斑点噪声分布近似为均值为1的逆高斯分布，雷达后向散射系数为常数，依据联合概率密度理论，sar图像像素i(x,y)指的是整幅图像的所有像素的集合服从逆高斯分布。假设sar图像指的是整幅图像区域并用ω表示，且区域ω被一闭合可变曲线c分为区域ω1和区域ω2两类区域，假设任意两个不同像素位置(x1,y1)≠(x2,y2)，且假设像素i(x1,y1)指的是整幅图像中给定位置的像素值与像素i(x2,y2)相互独立，区域ω1和区域ω2没有重合部分，sar图像服从逆高斯分布p(i(x,y))，基于最大似然估计分割区域ω1和区域ω2，则得到：

对ψig求最小时等价于求区域ω1最大似然估计，从而实现从ω区域中分割ω1和ω2，引入逆高斯概率密度函数pig(i(x,y))得到如下表达式：

从整幅sar图像考虑表达式由于其不随分割区域ω1和ω2的变化影响ψig的结果，因此，式(3)简化为：

用e(ui(x,y),λi(x,y))表示区域ω的能量泛函，则ψig1表达式改写为：

假设区域ω1为目标区域，区域ω2为背景区域，在sar图像海岸线检测中，海洋和陆地区域均作为目标区域，不必特意指定。为了很好地分割两区域，引用mumford-shah模型中轮廓曲线c长度项length(c)作为能量泛函的约束项，在能量泛函最小化过程中尽可能使轮廓曲线c贴近真实目标区域边界，同时保证曲线c的长度尽可能短，能量泛函e(ui(x,y),λi(x,y))改写为如下表达式：

其中，θ表示与sar图像像素无关的区域能量泛函的权重系数，μ表示与sar图像像素无关的长度约束项的权重系数。为了使得该能量泛函能够适应拓扑结构变化，将水平集方法引入区域能量泛函e(ui(x,y),λi(x,y))，则轮廓曲线c隐式地表现为水平集函数中的任意一个等值曲线，假设水平集函数用一个三维连续曲面函数φ(x,y)表示，且将水平集函数等于零对应的曲线位置作为轮廓曲线c随轮廓内外的均值演化后的形状，直至轮廓曲线演化到真正目标边缘处使得能量泛函e(ui(x,y),λi(x,y))处于局部极小，在sar图像海岸线检测中引入水平集函数约束项，选取符号距离函数作为水平集演化方程即水平集函数的更新方程的约束项，并且水平集演化方程初始化时规则曲线作为初始轮廓曲线；所述的规则曲线为圆形或长方形；

由于符号距离函数具有特性，为了使水平集函数在更新过程中保持符号距离函数特性，引入在迭代过程中无需重复初始化的水平集函数作为能量泛函的约束项，通常称为能量泛函惩罚项，其具体表达式如下：

在式(6)中引入水平集函数，并将式(7)引入式(8)得能量泛函如下：

其中，m1(φ(x,y))＝h(φ(x,y))，m2(φ(x,y))＝1-h(φ(x,y))，为梯度，h(φ(x,y))表示以水平集函数为参数的heaviside函数，当φ(x,y)为正值时m1(φ(x,y))有效，λ1(x,y)、u1(x,y)和φ(x,y)为正值时对应sar图像像素i(x,y)参与水平集函数演化，当φ(x,y)为负值时m2(φ(x,y))有效，λ2(x,y)、u2(x,y)和φ(x,y)为负值对应像素点的sar图像像素i(x,y)参与水平集函数演化。v表示水平集函数惩罚项的权重系数，至此，将sar图像区域ω1和区域ω2的分割问题转化为求区域能量泛函e(λi(x,y),ui(x,y),φ(x,y))的最小值问题。

由式(8)可知，能量泛函e(λi(x,y),ui(x,y),φ(x,y))取极小值不仅和sar图像像素i(x,y)有关，还取决于λi(x,y)、ui(x,y)和φ(x,y)的大小，在水平集函数演化过程中假设λi(x,y)和水平集函数φ(x,y)均为已知常数，ui(x,y)未知，且能量泛函e(λi(x,y),ui(x,y),φ(x,y))随ui(x,y)的变化而变化，由式(8)可知，能量泛函e(λi(x,y),ui(x,y),φ(x,y))中第二项和第三项表达式均不涉及ui(x,y)，因此将式(8)简化为：

将看作参数为1/ui(x,y)、λi(x,y)和φ(x,y)的能量泛函，假设λi(x,y)为大于零的已知常数，当能量泛函存在极小值时，能量泛函与1/ui(x,y)满足以下关系：

通过求式(10)得u1(x,y)和u2(x,y)的表达式如下：

同理，设φ(x,y)和ui(x,y)为已知常数，能量泛函e(λi(x,y),ui(x,y),φ(x,y))中λi(x,y)未知，为求解λi(x,y)，对能量泛函关于λi(x,y)求偏导表示为：

c、计算λ1(x,y)和λ2(x,y)；

将步骤b中计算得到的u1(x,y)、u2(x,y)和sar图像像素i(x,y)代入式(14)和式(15)分别计算λ1(x,y)和λ2(x,y)；

由式(11)-(13)得λ1(x,y)和λ2(x,y)的表达式如下：

最后，假设u1(x,y)和u2(x,y)、λ1(x,y)和λ2(x,y)均已知，求式(8)能量泛函对应的水平集演化方程。基于经典的求水平集演化函数方法，假设能量泛函满足如下表达式：

e(λi(x,y),ui(x,y),φ(x,y))＝∫∫ωf(φ(x,y),φx(x,y),φy(x,y))dxdy(16)

假设式(16)中f(φ(x,y),φx(x,y),φy(x,y))由f1(φ(x,y),φx(x,y),φy(x,y))、f2(φ(x,y),φx(x,y),φy(x,y))与f3(φ(x,y),φx(x,y),φy(x,y))三个泛函之和组成，且fi(φ(x,y),φx(x,y),φy(x,y))(i＝1,2,3)的表达式分别如下：

为了简化fi(φ(x,y),φx(x,y),φy(x,y))书写，将其用fi(·)表示，对式(17)求euler-lagrange方程得：

同理，对式(18)、式(19)求euler-lagrange方程分别得：

其中，div(·)表示散度，δ表示拉普拉斯算子，δ(φ(x,y))为h(φ(x,y))关于φ(x,y)的导数，在求能量泛函e(λi(x,y),ui(x,y),φ(x,y))对应的水平集演化方程之前，先对式(22)进行以下说明，假设式(22)对应的水平集函数为φp(x,y)则φp(x,y)满足如下关系：

由于δφ(x,y)写为对求梯度，即式(23)改写如下：

令

则式(23)看作扩散系数为α关于φ(x,y)的扩散方程

当时，得到α＞0，属于正向扩散，水平集函数趋于平滑，减小；反之，当1时，得到α＜0，属于反向扩散，增大，因此惩罚项对水平集函数演化起到调节作用，使得保持为符号距离函数特性。能量泛函表达式最小时对应的水平集函数演化方程满足如下等式：

将式(20)-(22)代入(25)得到基于逆高斯分布主动轮廓模型的水平集演化方程：

其中，

d、基于逆高斯分布主动轮廓模型的数值求解

由于式(14)、式(15)和式(26)中图像像素i(x,y)均作为分母参与运算，将图像所有像素都增加单位像素，避免当原图像像素为零时作为分母没有意义，同理，由于u1(x,y)和u2(x,y)也作为分母参与运算，所以在计算u1(x,y)和u2(x,y)时分别增加小的正数，以避免分母为零的运算出现。这里λ1(x,y)和λ2(x,y)因取对数运算参与式(26)水平集演化，为避免为零也需要增加小正数处理。沿用水平集演化方程数值求解的定义和符号说明，水平集函数的时间求偏导采用前向差分形式，即

水平集函数对空间距离求偏导采用中心差分的形式，用和分别表示水平方向和垂直方向的偏导数，则

式(26)水平集演化方程对应的数值求解如下：

其中θ、μ和ν为权值参数，表示长度规则项的数值求解项，表示水平集函数惩罚项的数值求解项；和的表达式如下：

其中η为小正数，使得在水平集函数梯度模为零时仍然有意义。

将步骤c和d得到的u1(x,y)、u2(x,y)和λ1(x,y)、λ2(x,y)代入式(30)更新水平集函数，若记当前水平集函数为φⁿ(x,y)，则更新后的水平集函数记为φⁿ⁺¹(x,y)

e、检查水平集函数是否达到收敛条件，若收敛，停止水平集方程演化，并将得到的轮廓曲线c作为最终轮廓边界，若不收敛，返回步骤b，直到满足水平集函数收敛条件。

与现有技术相比，本发明具有以下有益效果：

针对海洋区域灰度图像像素分布均匀、陆地区域灰度图像像素分布不均匀、且海陆分界区域存在弱边缘、整体图像像素偏向于低灰度值区域的envisat卫星图像，本发明提出了一种基于逆高斯分布的区域主动轮廓模型。首先分析了基于gamma分布的区域主动轮廓模型原理，然后基于该模型的原理和能量泛函构建方法，假设噪声服从单位均值的逆高斯分布，雷达后向散射系数为常数，通过最大似然估计方法构建基于逆高斯分布的区域能量泛函，然后引入水平集函数、长度规则项和惩罚项得到基于逆高斯分布的区域主动轮廓模型，进一步理论推导该模型对应的水平集函数演化方程。通过实验验证了本发明不仅可以处理海洋区域图像像素分布均匀、陆地区域图像像素分布不均匀、整体图像像素偏向于低灰度区域的envisat卫星图像，还适用于海陆边界区域存在弱边缘的envisat卫星图像，并且本发明的运行时间比基于gamma分布的区域主动轮廓模型运行时间要少。因为用传统的gamma分布去模拟最后求得的能量泛函对于水平集函数的偏导数要小于逆高斯分布的能量泛函对于水平集函数的偏导数，而更新函数是与能量泛函对水平集的偏导数成正比例的所以后者要比前者迭代的次数更少。

附图说明

图1为本发明的模型和gamma模型检测海岸线的实验结果。

图2为本发明的模型和gamma模型检测海岸线的实验结果对比。

图3为rmse曲线图。

图4为本发明的模型和gamma模型运行时间比较图。

图5为本发明的方法流程图。

具体实施方式

下面结合附图和实验结果对本发明进行进一步地描述。

本发明基于windows7系统、matlab2012a作为实验平台，实验使用envisat-1和envisat-2卫星合成的sar图像作为研究数据，图像大小选取255×255像素，采用imp_v/v模式和imp_h/h模式。

1、envisat图像对比和分析

为了验证本发明的有效性，按照图5所示流程，选取基于gamma分布的区域主动轮廓模型作为对比实验对象，由于本发明和基于gamma分布的区域主动轮廓模型都是基于全局角度考虑问题，理论上检测结果不受初始轮廓的位置影响。

基于gamma分布的区域能量模型中θ＝1，μ取值范围为0.05～0.5，图1-2中μ＝0.2，由于在能量泛函中增加惩罚项不仅使得水平集演化方程演化速度较为缓慢，且对实验结果没有明显改善，所以对比实验计算方法中不加入水平集函数的惩罚项。本发明中参数λ1＝λ2＝1，μ取值范围为0.1～1.0，图1-2中μ＝0.5，ν的取值范围为0.1～1.0，图1-2中ν＝0.2，本发明选取5000次作为终止条件，基于gamma分布的区域主动轮廓模型实验选取10000次作为终止条件。

由合成孔径雷达成像原理可知，对于实际均匀目标经雷达散射成像后像素低，对于海面区域分布均匀陆地区域分布不均匀且海陆分界区域不存在明显弱边缘的sar图像，通过gamma分布构建能量泛函和本发明均可取得比较满意的检测结果，当海洋区域存在小岛或船只时，本发明和对比实验gamma分布的模型均将其检测出来。影响检测结果的准确性，如图1所示，图1中的图像编号从上至下依次为1-3号。在图2中看到对于海洋和陆地分界区域存在弱边缘、且海洋区域和陆地区域灰度对比度不明显，且海洋区域均匀，陆地区域不均匀的sar图像，利用gamma分布构建能量模型检测的结果没有本发明检测结果理想。图2中的图像编号从上至下依次为4-7号。

2、实验误差分析和运行时间比较

本发明同样参考shengf和wangl等通过计算均方根误差分析精度的方法分析本发明检测海岸线的准确性，选取人工标注海岸线结果为基准，使用计算机检测海岸线与人工标注海岸线比较计算均方根误差，并用rmse表示均方根误差。本发明选取同样大小的7幅envisat卫星图像计算rmse，结算结果如图3所示。

图3中的图像1-7分别对应图1-2中的1-7号图像，当海洋和陆地区域图像像素对比度低，且图像像素主要集中在低灰度值时，本发明可以很好地拟合边缘，通过图3表明本发明和人工标注方法得到的海岸线均方误差平均在2个像素左右。

同时，基于本发明与基于gamma分布的区域主动轮廓模型均能检测出理想海岸线的前提下，选取7幅envisat图像比较两种计算方法中轮廓曲线到达真实海岸线基本稳定时所消耗时间，两种计算方法均基于本发明参数假设条件粗略统计运行时间，运行时间结果如图4所示。

通过分析图4的本发明与基于gamma分布的区域主动轮廓模型计算方法检测海岸线的系统运行时间，本发明检测海岸线所需要的时间较少，而基于gamma分布的区域主动轮廓模型对于envisat卫星图像检测海岸线所需要的时间较多。