一种高光谱溢油图像的有效分割方法与流程

技术特征：

1.一种高光谱溢油图像的有效分割方法，其特征在于:

一、定义初始水平集函数及其他相关函数

初始水平集函数一般定义为：

φ(x,t＝0)＝±c (1)

其中，t＝0表示初始轮廓C，x为自变量，φ为水平集函数，若点x在轮廓内函数取值为-c，反之则取值为+c，

定义Heaviside函数H及其一阶导数δ：

$\begin{array}{c}H\left(z\right)=\left\{\begin{array}{c}1,z\≥0\\ 0,z\<0\end{array}\right.\\ \mathrm{\δ}\left(z\right)=\frac{d}{dz}H\left(z\right)\end{array}---\left(2\right)$

其中，H代表Heaviside函数，δ代表H函数的一阶导数，z为实数自变量，定义函数H_ε,δ_ε如下：

$\begin{array}{c}{H}_{\mathrm{\ϵ}}\left(z\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan (\frac{z}{\mathrm{\ϵ}}\left)\right)\\ {\mathrm{\δ}}_{\mathrm{\ϵ}}\left(z\right)\frac{d}{dz}{H}_{\mathrm{\ϵ}}\left(z\right)\end{array}---\left(3\right)$

其中，z为实数自变量，ε定义为一个很小的实数，H代表Heaviside函数，δ代表H函数的一阶导数在模型能量函数中，用函数H_ε,δ_ε分别近似表示H,δ；

二、结合Fisher准则得到新的拟合项

Fisher准则定义一个好的分类方法应该达到最大的类间间距和最小的类内差异，用公式可以表示为：

$min\left\{\frac{{\mathrm{\σ}}_1^2+{\mathrm{\σ}}_2^2}{{(\stackrel{\‾}{c_1}-\stackrel{\‾}{c_2})}^2}\right\}---\left(4\right)$

其中，分别为两类样本的平均向量，分别为两类样本的统计方差，当J_fisher取最大值时达到最佳分类效果；

原CV模型中的拟合项是类内方差和，目标函数为类内方差最小，将Fisher准则的基本思想引入到CV模型，将原模型代替类内方差，在此基础上加入类间方差)中，则改进后的拟合项可以表示为：

$\begin{array}{c}{J}_{fisher}=\frac{1}{{({\stackrel{\‾}{c}}_1-{\stackrel{\‾}{c}}_2)}^2}\[{\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{|I(x,y)-{\stackrel{\‾}{c}}_1|}^2{H}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x,y\left)\right)dxdy\\ +{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|I(x,y)-{\stackrel{\‾}{c}}_2|}^2(1-H_{\mathrm{\ϵ}}(\mathrm{\φ}\left(x,y\right)\left)\right)dxdy\]\end{array}---\left(5\right)$

其中，J_fisher结合Fisher原则的拟合项函数，(x,y)为像素点坐标值，I表示高光谱图像，分别为轮廓曲线C内外的均值向量，λ₁,λ₂为内外拟合项的常量系数，Ω表示图像区域，H_ε,δ_ε分别近似表示H,δ，φ表示水平集函数；

三、构造边缘停止函数得到新的长度项

光谱角是指影像像元光谱矢量与样本参考光谱矢量之间的夹角，用来度量两向量之间的相似程度，与两矢量的大小无关，其计算公式如下：

$SAM={\cos }^{-1}\frac{A\·B}{\left|A\right|*\left|B\right|}={\cos }^{-1}\frac{\underset{i=1}{\overset{N}{\Σ}}{A}_i*B_i}{\sqrt{\underset{i=1}{\overset{N}{\Σ}}{A}_i^2}*\sqrt{\underset{i=1}{\overset{N}{\Σ}}{B}_i^2}}---\left(6\right)$

式中，SAM为光谱角，N为波段数，A＝(A₁,A₂,…,A_N)和B＝(B₁,B₂,…,B_N)分别代表样本空间中的两个影像像元的光谱值，当光谱角SAM越小，则表明两光谱矢量越接近；

边缘停止函数只需满足是一个正的关于光谱角梯度递减的函数，所以将其定义为：

$\begin{array}{c}\mathrm{\α}=\frac{1}{2}({\cos }^{-1}\frac{I(x,y)\·I(x+\mathrm{\Δ}x,y)}{||I(x,y)||*||I(x+\mathrm{\Δ}x,y)||}+{\cos }^{-1}\frac{I(x,y)\·I(x,y+\mathrm{\Δ}y)}{||I(x,y)||*||I(x,y+\mathrm{\Δ}y)||})\\ g_{sam}\left(\mathrm{\α}\right)=\frac{1}{1+\left|\▿\mathrm{\α}\right|}\end{array}---\left(7\right)$

式中，I表示高光谱图像，为光谱角梯度，(x，y)代表像素点的坐标值，g_sam(α)是光谱角梯度函数，α表示图像任意一点的像元向量I(x,y)相邻点的像元向量I(x+Δx,y)以及I(x,y+Δy)之间的光谱角距离，当处于图像边界时，光谱角梯度的大小取值较大，取值较小，

引入上述基于光谱角梯度信息的边缘函数g_sam(α)，则改进后的长度项可以表示为：

$L\left(\mathrm{\φ}\right)=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}{g}_{sam}\left(\mathrm{\α}\right)\mathrm{\δ}\left(\mathrm{\φ}\right(x,y\left)\right)|\▿\mathrm{\φ}(x,y)|dxdy---\left(8\right)$

其中，L代表长度项，Ω表示图像区域，g_sam(α)是光谱角梯度函数，(x，y)代表像素点的坐标值，δ表示Heaviside函数的一阶导数，φ为水平集函数，表示φ的增量，μ为长度项系数，

四、结合端元提取算法的改进

CV模型一般只能用于二类划分，即只能将图像分割为目标和背景两个区域，而现实情况是图像中往往含有多类地物，当图像中含有多类地物时，主动轮廓模型无法实现对感兴趣的特定地物的划分，因此，将CV模型与端元提取算法相结合，从而实现从含多种地物的图像中划分出特定地物，

相比于传统的端元提取算法，ATGP算法能够在无需先验信息的情况下，从目标图像中提取出端元向量，ATGP算法是一种基于非监督正交子空间投影理论的端元提取算法，该算法首先按照凸面几何学理论，将正交子空间中亮度最大的像元作为候选端元，从而得到初始端元向量，以此构造子空间，求取相应的正交子空间，并将所有像元投影到该正交子空间，求取正交子空间中投影最大的像元作为下一个端元向量；以此类推，直到找到指定个数的端元；

采用ATGP算法自动提取目标端元，首先用ATGP算法从图像中提取出特定类的端元向量，然后以该端元向量的位置为中心点设置初始轮廓，并且用该端元向量替换改进主动轮廓模型能量泛函中的轮廓曲线内部的均值向量，若设t为ATGP算法从高光谱图像中提取出的特定类端元向量，则此时改进模型的能量泛函为：

$\begin{array}{c}{E}_{\mathrm{\ϵ}}(t,{\stackrel{\‾}{c}}_2,\mathrm{\φ})=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}{g}_{sam}\left(\mathrm{\α}\right){\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x,y\left)\right)|\▿\mathrm{\φ}(x,y)|dxdy\\ +\mathrm{\ν}{\∫}_{\mathrm{\Ω}}{E}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x,y\left)\right)dxdy\\ +\frac{1}{{(t-{\stackrel{\‾}{c}}_2)}^2}({\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{\left|I\left(x,y\right)-t\right|}^2{E}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\left(x,y\right)\right)dxdy\\ +{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|I(x,y)-{\stackrel{\‾}{c}}_2|}^2(1-E_{\mathrm{\ϵ}}(\mathrm{\φ}\left(x,y\right)\left)\right)dxdy)\\ +\mathrm{\η}{\∫}_{\mathrm{\Ω}}{\left(\right|\▿\mathrm{\φ}|-1)}^{2}dxdy\end{array}---\left(9\right)$

其中，E_ε表示目标能量函数，α表示由(7)式决定的光谱角，表示表示水平集函数φ的增量，(x，y)代表像素点的坐标值，g_sam(α)是光谱角梯度函数，μ,ν,λ₁,λ₂均为常系数且满足μ≥0,ν≥0,λ₁,λ₂＞0，I为矢量图像中(x,y)处的像元向量，φ为水平集函数，t,分别为轮廓曲线C内部区域和外部区域的灰度均值向量，H_ε,δ_ε分别表示近似Heaviside函数及其一阶导数,Ω表示图像空间；

五、引入水平集正则项避免水平集函数重新初始化

为了避免水平集函数重新初始化，节约算法时间复杂度，进而提高算法效率，引入如下水平集正则项：

${\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\φ}|-1)}^{2}dxdy---\left(10\right)$

Ω表示图像空间,表示表示水平集函数φ的增量；

六、能量泛函最小化得到欧拉-拉格朗日方程

综上所述，新建模型的能量泛函为

$\begin{array}{c}{E}_{\mathrm{\ϵ}}(\stackrel{\‾}{c},{\stackrel{\‾}{c}}_2,\mathrm{\φ})=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}{g}_{sam}\left(\mathrm{\α}\right){\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x,y\left)\right)|\▿\mathrm{\φ}(x,y)|dxdy\\ +\mathrm{\ν}{\∫}_{\mathrm{\Ω}}{E}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right(x,y\left)\right)dxdy\\ +\frac{1}{{(\stackrel{\‾}{c}-{\stackrel{\‾}{c}}_2)}^2}({\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{\left|I\left(x,y\right)-\stackrel{\‾}{c}\right|}^2{H}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\left(x,y\right)\right)dxdy\\ +{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|I(x,y)-{\stackrel{\‾}{c}}_2|}^2(1-H_{\mathrm{\ϵ}}(\mathrm{\φ}\left(x,y\right)\left)\right)dxdy)\\ +\mathrm{\η}{\∫}_{\mathrm{\Ω}}{\left(\right|\▿\mathrm{\φ}|-1)}^{2}dxdy\end{array}---\left(11\right)$

其中，E_ε表示目标能量函数，Ω表示图像空间，I为矢量图像，μ,v,λ₁,λ₂均为常系数且满足μ≥0,v≥0,λ₁,λ₂＞0，I(x,y)为矢量图像，φ为水平集函数，表示表示水平集函数φ的增量，分别为轮廓曲线C内部区域和外部区域的均值向量，t为对应特定类的端元向量，或t，g_sam(α)是光谱角梯度函数，H_ε,δ_ε分别表示近似Heaviside函数及其一阶导数；

将能量泛函最小化得到欧拉-拉格朗日方程，将其动态化，得到相应的梯度下降流为

$\frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\[{\mathrm{\μg}}_{sam}\left(\mathrm{\α}\right)div\left(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}\right)-\mathrm{\ν}-\frac{1}{{(\stackrel{\‾}{c}-{\stackrel{\‾}{c}}_2)}^2}({\mathrm{\λ}}_1{\left(I-\stackrel{\‾}{c}\right)}^2-{\mathrm{\λ}}_2{\left(I-{\stackrel{\‾}{c}}_2\right)}^2)\]+\mathrm{\η}({\▿}^2\mathrm{\φ}-div(\frac{\▿\mathrm{\φ}}{\left|\▿\mathrm{\φ}\right|}\left)\right)---\left(12\right)$

其中，或t，表示表示水平集函数φ的增量，α表示由(7)式决定的光谱角，μ代表水平集正则项系数，μ,v,λ₁,λ₂均为常系数且满足μ≥0,v≥0,λ₁,λ₂＞0，I为矢量图像，φ为水平集函数，分别为轮廓曲线C内部区域和外部区域的均值向量，t为对应特定类的端元向量，或t，g_sam(α)是光谱角梯度函数，δ_ε表示近似Heaviside函数的一阶导数，div(x)表示散度；

七、设置各个参数

迭代时间步长Δt，函数H_ε的参数ε，长度项系数μ，面积项系数v，拟合项系数λ₁,λ₂，水平集正则项系数η；

八、选择显示波段以及初始轮廓

从几个波段中选择出一个对比度相对较高的波段作为显示最终轮廓的波段，然后设置初始轮廓，一般选择以特定圆心的圆或均匀分布半径较小的多个圆；

九、显示分割结果图

将方程离散化，并按设定的次数进行迭代，显示轮廓最终在图像中的位置；

十、计算各种分割精度评价指标

先对所得的分割结果图像进行二值化，然后结合真实的二值化图像，计算出错分率、漏分率、错分漏分率之和、Kappa系数以及总体精度等指标；

十一、对分割结果精度进行比较、评价

而对于模拟高光谱图像和已知真实地貌的高光谱图像，我们可以通过计算并比较相应的评价指标来评价其分割效果；而对于未知真实地貌的高光谱图像，只能通过人眼直观评价。