CN102609709A - 一种基于极化sar数据融合的海面溢油分割方法 - Google Patents

Abstract

本发明公开了合成孔径雷达海洋遥感应用技术领域中的一种极化SAR数据融合的海面溢油图像分割方法。该方法首先构造基于分割区域最大后验概率准则的活动轮廓能量泛函，并将分割区域的分布表示成Gibbs先验概率模型，然后将活动轮廓模型嵌入到高维的水平集函数中，利用Euler-Lagrange公式得到发展方程，模型中包含CFAR边缘检测加权的边界长度项和融合数据统计距离项。同时，本发明提出了确定演化参数的方法，结合油膜衰减特性初始化的水平集函数，通过方程演化得到分割结果。发明可以有效的融合各种SAR、极化SAR数据，实现了对海面暗区域的自动分割。

Description

一种基于极化SAR数据融合的海面溢油分割方法

技术领域

本发明属于合成孔径雷达(SAR)海洋遥感应用技术领域，尤其涉及一种基于极化SAR数据融合的海面溢油分割方法。

背景技术

随着人们在海上的活动越来越频繁，海面浮油已经成为了人们日益关注的重要环境污染问题。目前，海上溢油污染主要是由两种因素造成的，一是由重大事故造成的污染，如：油船事故和油气田平台的漏油事故；二是由油船在运输过程中的非法排放引起。由于卫星系统可以提供连续的海洋表面监测数据并且可以进行大范围的测绘，因此很适宜对海洋进行监测。对海况进行监测有很多遥感手段：光学、红外、以及多波段雷达等。其中，合成孔径雷达(Synthetic Aperture Radar，缩写为SAR)由于可以全天时、全球范围内、可穿透、高分辨成像，且成像不受天气海洋等因素的影响，已经成为海洋监测的主要手段。在SAR的海面溢油监测任务中，油膜由于具有粘弹性，对海洋表面的重力毛细波具有衰减作用，在SAR图像上表现为暗区域。海洋油膜检测一般首先需要对暗区域进行检测和分割，之后提取典型特征，将油膜与其他具有低后向散射特性的海洋现象(油膜类似区域)加以区分。因此，在海面溢油检测中，海面暗区域的自动检测通常是溢油检测的第一步，分割出来的暗区域可以作为后续油膜识别的待选区域。本发明结合油膜的衰减特性以及统计特性，主要关注于极化SAR融合数据的海面溢油分割问题，并提供一种适用于各种SAR数据的统一海面暗区域自动分割的新方法。

目前，关于海面暗区域的分割与分类方法大都是基于单极化SAR强度图像的方法，在海面油膜的检测中主要依据的是海面浮油的衰减特征。近年来，随着全极化SAR系统的广泛应用，获取的海洋表面油膜监测数据越来越多，发展了很多利用极化SAR数据对海面油膜进行检测的方法。从数据信息的获取中看，单极化SAR系统仅能提供目标的一个通道的测量数据，全极化SAR系统可以提供目标完整的后向散射信息，双极化SAR系统对目标观测所获得的数据量处于两者之间。获取的目标观测信息越多越有利于对目标的检测，因此利用多波段、多极化、多时相的SAR数据以及极化SAR融合数据可以提高海面油膜的检测概率。但是，此项技术的关键问题就是如何融合不同数据源的SAR数据，目前关于极化SAR数据融合技术方面的研究还比较少。

在图像的自动分割研究中，活动轮廓法由于融入了人类的感知理解机制，对图像的分割具有较好的效果。活动轮廓模型是基于变分法以及偏微分方程的一种分割模型，通过在可以活动的轮廓曲线上定义一种能量泛函，将图像分割问题转化为能量最小化问题，通常将此能量函数隐式地嵌入水平集函数中，这就是变分水平集分割方法。对于光学图像的分割，绝大多数活动轮廓方法采用了图像分片光滑假设以及高斯噪声模型，因此发展成了统计活动轮廓模型。很多人用此加性噪声模型对SAR图像进行分割，但此模型不符合SAR数据的统计分布且对SAR图像以及极化SAR图像没有普适性。SAR数据的统计模型是由雷达回波信号衰落引起的乘性斑点噪声(Speckle)模型。为了实现较好的图像分割效果，必须使分割模型与实际数据的统计分布一致。Ayed等人分别于2005年和2006年发展了基于Gamma分布以及Wishart分布的单通道SAR强度图像以及全极化SAR散射相干矩阵的水平集分割方法。这些分割模型中包含了区域的统计信息以及分割边界的长度信息，但是没有利用到目标的边缘信息。图像的边缘是目标的基本纹理特征之一，活动轮廓模型最初就是利用目标的边缘信息驱动演化方程，从而实现对目标的分割。边缘信息是海面油膜的基本特征之一，若在分割模型中加入海面浮油的边缘信息可以提高油膜的分割精度。统计活动模型通常包含几个部分，各部分对最终图像分割结果的影响不同，只有使各组成部分成比例才能有较好的分割效果。但在目前的基于活动轮廓模型的图像分割方法中，还没有给出过明确的参数选择解决方案。另外，极化SAR数据融合可以改善SAR系统的应用性能，如地物分类、目标检测与识别等，但目前的方法还没有关于对极化融合数据的统计活动轮廓分割模型方面的讨论。

发明内容

本发明的目的在于，提供一种基于活动轮廓模型的极化SAR数据融合的海面溢油图像分割方法，用以解决基于加性噪声的区域统计模型对SAR图像以及极化SAR图像没有普适性的问题，同时本发明可以有效地融合不同的SAR数据源，实现了海面溢油的自动分割。

为实现上述目的，本发明提供的技术方案是，一种基于极化SAR数据融合的海面溢油图像方法，其特征是所述方法包括：

步骤1：选取不同合成孔径雷达数据源对海面溢油图像的监测数据；

步骤2：根据选取的监测数据确定所述图像的散射相干矩阵；

步骤3：依据活动轮廓法对所述图像进行分割。

所述散射相干矩阵为分块对角矩阵，具体形式为

其中，C_i为1×1、2×2或3×3的半正定Hermitian矩阵，i＝1，2，...，n，n为半正定Hermitian矩阵的个数，并且C_i满足P(T|V)＝P₁(C₁|V₁)P₂(C₂|V₂)...P_n(C_n|V_n)，P(·)是矩阵T的概率密度分布函数，V＝E(T)是散射相干矩阵T的均值，P_i(·)是矩阵C_i的概率密度分布函数，V_i＝E(C_i)是矩阵C_i的均值。

所述步骤3具体是：

步骤101：利用活动轮廓曲线c将所述图像分割成2个区域；

步骤102：利用公式g＝(1+k)/(1+k(R/T_f)²)计算每个分割区域的边缘截止函数；其中，g是边缘截止函数，R是由恒虚警率边缘检测算子得到的所述图像边缘密度，T_f由虚警率得到的边缘检测阈值，k是设定常数；

步骤103：根据公式

$E_{\mathrm{POL}}(V_1,V_2,c)=\mathrm{ulength}\left(c\right)\·g(x,y)+{\mathrm{\λ}}_1{\∫}_{\mathrm{inside}\left(c\right)}\left(\log \right|V_1|+\mathrm{Tr}(V_1^{-1}\·T(x,y)))\mathrm{dxdy}+$

${\mathrm{\λ}}_2{\∫}_{\mathrm{outside}\left(c\right)}\left(\log \right|V_2|+\mathrm{Tr}(V_2^{-1}\·T(x,y)))\mathrm{dxdy}$

计算所述图像分割后的分割能量函数；其中，c为活动轮廓曲线，u、λ₁和λ₂为设定常数，V₁和V₂分别是分割区域的均值，length(c)表示分割边界的长度，∫_inside(c)(·)dxdy和∫_outside(c)(·)dxdy分别表示对活动轮廓曲线c内和活动轮廓曲线c外所有的像素点进行积分，Tr(·)表示求矩阵的迹，(x，y)是图像中的像素点的坐标，T(x，y)表示像素点的散射相干矩阵；

步骤104：将所述分割能量函数嵌入到水平集函数φ(x，y)中，令λ＝λ₁＝λ₂，得到水平集分割能量泛函为

$E_{\mathrm{POL}}(V_1,V_2,\mathrm{\φ})=U{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)g(x,y)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}$

$+\mathrm{\λ}\left(\begin{array}{c}{\∫}_{\mathrm{\Ω}}\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+\\ {\∫}_{\mathrm{\Ω}}\left(\log \right|V_2|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}\end{array}\right)$

其中，Ω是待分割的图像区域，δ(·)表示Dirac函数，

表示求φ(x，y)的梯度，H(·)是Heaviside函数；

步骤105：初始化水平集函数并设定水平集函数的发展方程的参数；

所述初始化水平集函数具体是设定 ${\mathrm{\φ}}_0(x,y)=\left\{\begin{array}{cc}d& (x,y)\∈{\mathrm{\Ω}}_0\\ 0& (x,y)\∈{\∂\mathrm{\Ω}}_0\\ -d& (x,y)\∈\mathrm{\Ω}-{\mathrm{\Ω}}_0\end{array}\right.;$

所述设定水平集函数的发展方程的参数具体是设定t＝0、λ＝10、k＝9并且利用公式 $u\left({\mathrm{\Σ}}_{{\mathrm{\Ω}}_h}g\right)/\left(\mathrm{\λ}{\mathrm{\Σ}}_{{\mathrm{\Ω}}_h}\mathrm{Tr}\left(V^{-1}T\right)\right)=0.25$ 计算u的值，Ω_h为选定的均匀区域；其中，t为迭代初值；

步骤106：计算每个分割区域的均值；其计算公式分别为

$V_1=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}$ 和 $V_2=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}};$

步骤107：利用公式

${\mathrm{\φ}}_{t+1}={\mathrm{\φ}}_t+\mathrm{dt}\×\frac{\∂\mathrm{\φ}}{\∂t}$

更新水平集函数；其中，t＝0，1，2...为迭代次数，dt为步长并设定dt＝0.5；

为水平集函数的发展方程

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(u\mathrm{div}(g\·\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|})-\mathrm{\λ}(\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T\right))-\left(\log \right|V_2|+\mathrm{Tr}\left(V_2^{-1}T\right))))$

其中，div(·)表示求向量的散度；

步骤108：判断|φ_t+1-φ_t|＜ε是否成立，如果|φ_t+1-φ_t|＜ε成立，则执行步骤109；否则，令t＝t+1，返回步骤106；其中，ε为设定阈值；

步骤109：分割结束。

本发明通过融合不同SAR数据源的数据建立散射相干矩阵，并依据最大后验概率准则建立了基于区域统计的活动轮廓能量泛函，将该活动轮廓模型用于海面溢油图像的分割中，分割的结果可作为后续油膜判定的依据，降低了油膜检测的虚警率；同时，本发明适用于任何形式的SAR数据、极化SAR数据以及极化SAR融合数据中，是SAR海面溢油区域自动检测与分割的通用方法。

附图说明

图1是基于极化SAR数据融合的海面溢油分割方法流程图；

图2是仿真油膜数据的总功率图；

图3是真实的ENVISAT/ASAR传感器获取的海面溢油形状示意图；

图4是三个通道强度融合数据的初始化活动轮廓曲线以及水平集函数计算结果图；

图5是三个通道强度融合数据的统计活动轮廓模型的水平集函数最终演化结果图；

图6是三个通道强度融合数据的统计活动轮廓模型最终分割结果图；

图7是单极化强度数据的统计活动轮廓模型分割结果图；其中，(a)是HH通道强度数据分割结果图，(b)是HV通道强度数据分割结果图，(c)是VV通道强度数据分割结果图；

图8是双极化强度融合数据的统计活动轮廓模型分割结果图；其中，(a)是HH通道+HV通道的数据分割结果图，(b)是HV通道+VV通道的数据分割结果图；

图9是三个通道强度数据的监督ML分类结果图。

具体实施方式

下面结合附图，对优选实施例作详细说明。应该强调的是，下述说明仅仅是示例性的，而不是为了限制本发明的范围及其应用。

实施例1

图1是基于极化SAR数据融合的海面溢油分割方法流程图。图1中，本发明提供的基于极化SAR数据融合的海面溢油分割方法包括：

步骤1：选取不同合成孔径雷达数据源的海面溢油图像监测数据。

步骤2：确定所述图像的散射相干矩阵。

合成孔径雷达(SAR)通过发射一种电磁波，对目标观测的测量值就包含在雷达后向散射波中，由于受到相干斑噪声(Speckle)的影响，单视的SAR强度图像服从负指数分布，多视的SAR强度图像σ服从Gamma分布，记为

G(·)为Gamma函数，且有

$P\left(\mathrm{\σ}\right)={\left(\frac{L}{u_I}\right)}^L\frac{1}{\mathrm{\Γ}\left(L\right)}{\mathrm{\σ}}^{L-1}\exp (-\frac{\mathrm{L\σ}}{u_I})$

其中，Γ(L)＝(L-1)！，L表示视数，u_I是统计均值，

表示分布的标准方差。

全极化SAR以及双极化SAR系统可以同时获得目标的幅度以及相位信息，设目标在多极化方式下，其复散射矢量

具有p个自由度：

$\stackrel{\→}{k}={\left[\begin{array}{cccc}{k}_1& k_2& ...& k_p\end{array}\right]}^T$

其中，[·]^T表示转置，在后向散射互易的条件下，全极化SAR数据的p＝3，双极化SAR数据的p＝2。因此，多极化SAR数据的后向散射特性一般用复相关(或相干)矩阵表示：

$C=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\stackrel{\→}{k}}_i\×{\stackrel{\→}{k}}_i^H$

其中，H表示转置共轭。散射相关矩阵C是半正定的Hermitian矩阵，单视的C矩阵在自由度为p×p时服从复Gaussian分布，L视的C矩阵在自由度为p×p时服从复Wishart分布，记为C∈W_c(p，L，V)，其概率密度分布为

$P\left(C\right|V)=\frac{L^{\mathrm{pL}}{\left|C\right|}^{L-p}\exp (-\mathrm{LTr}\left(V^{-1}C\right))}{K(L,p){\left|V\right|}^L}$

其中，V＝E(C)表示求解散射相干矩阵C的均值，Tr(·)表示求矩阵的迹，K(L，P)是归一化系数，K(L，p)＝π^0.5p(p-1)Γ(L)…Γ(L-p+1)。从上式中可以看出，当p＝1时，复Wishart分布即变成单通道SAR强度数据的Gamma分布。因此，可以将多极化SAR数据外积的统计分布看成是一种统一的概率密度分布函数。

假设L视的散射相关矩阵T_x服从均值为V的复Wishart分布，根据最大似然函数准则可以得到目标T_x到其类中心的复Wishart距离为，

d_W(T_x，V)＝ln(|V|)+Tr(V^-1T_x)

此统计距离可以被应用到单极化SAR强度数据、双极化SAR数据以及全极化SAR数据的散射相干矩阵中。当同时具有多个波段、或者多时相的极化SAR数据、或者多通道的SAR强度数据，则可以将散射相关矩阵T统一表述成如下形式，

其中，C_i为1×1、2×2或3×3的半正定Hermitian矩阵，i＝1，2，...，n，n为半正定Hermitian矩阵的个数。此时，T是分块对角阵，但将不再服从复Wishart分布或者Gamma分布。设极化融合数据的散射相关矩阵形式如上所示，T是分块对角阵，C_i是其子矩阵。假设C_i之间是相互独立的，则极化SAR融合数据T的似然函数可以表示成，

P(T|V)＝P₁(C₁|V₁)P₂(C₂|V₂)...P_n(C_n|V_n)

其中，V＝E(T)，E(·)表示求解散射相干矩阵T的均值，P_i(·)是矩阵C_i的概率密度分布函数，V_i＝E(C_i)是矩阵C_i的均值。利用矩阵的迹以及行列式的性质，可以得到极化SAR融合矩阵T的统计距离形式为，

$d(T,V)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\ln \left(\right|V_i\left|\right)+\mathrm{Tr}\left(V_i^{-1}{C}_i\right)=\ln \left(\right|V\left|\right)+\mathrm{Tr}\left(V^{-1}T\right)$

上式与复Wishart距离具有相同的形式。这说明SAR数据、极化SAR数据以及极化SAR融合数据虽然具有不同的概率密度分布函数，但是它们的统计距离形式相同。因此，利用此性质可以融合不同的SAR数据源。

步骤3：建立区域统计活动轮廓能量函数模型，并对所述图像进行分割，其具体步骤如下：

步骤101：利用活动轮廓曲线c将所述图像分割成2个区域。

步骤102：利用公式g＝(1+k)/(1+k(R/T_f)²)计算每个分割区域的边缘截止函数；其中，g是边缘截止函数，R是由恒虚警率边缘检测算子得到的所述图像边缘密度，T_f由虚警率得到的边缘检测阈值，k是设定常数。

假设将整个图像被分割成m≥2个区域，区域均值V_i(i＝1，2，...，m)服从Gibbs先验分布，则有P(V_i)＝Z^-1exp(-uφ(V_i))。其中，u和Z是归一化的常数。可以将分割区域函数φ(V_i)表述成分割边界的函数，即∑_{i＝1，...，m}φ_i＝g(x，y)∮_cds。其中，∮_cds表示分割边界的长度，g(x，y)是目标边缘的函数，c表示目标区域的分割边界。设由恒虚警率(Constant False Alarm Rate，CFAR)边缘检测算子得到的SAR图像边缘密度为R(x，y)，g(x，y)可以定义为g＝(1+k)/(1+k(R/T_f)²)。其中，k是大于1的数，T_f是由虚警率得到的边缘检测阈值。g(x，y)是边缘停止函数，用它对分割曲线进行加权可以使得水平集函数的扩散呈现各向异性的扩散，即在非边缘区域的扩散速度加快，在边缘区域的扩散速度减慢。

步骤103：根据公式

$E_{\mathrm{POL}}(V_1,V_2,c)=\mathrm{ulength}\left(c\right)\·g(x,y)+{\mathrm{\λ}}_1{\∫}_{\mathrm{inside}\left(c\right)}\left(\log \right|V_1|+\mathrm{Tr}(V_1^{-1}\·T(x,y)))\mathrm{dxdy}+$

${\mathrm{\λ}}_2{\∫}_{\mathrm{outside}\left(c\right)}\left(\log \right|V_2|+\mathrm{Tr}(V_2^{-1}\·T(x,y)))\mathrm{dxdy}$

计算所述图像分割后的分割能量函数。

假设图像被活动轮廓曲线c分成了2个部分，根据最大后验概率(MAP)准则构造2相位的分割能量函数F(c)，使得F(c)最大，

$F\left(c\right)=\ln \left(P\left(T\right|c)\right)-\mathrm{u\φ}\left(c\right)-\ln \left(P\left(T\right)\right)+\ln \left(Z\right)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\ln \left(P\left(T\right|V_i)\right)-\mathrm{u\φ}\left(c\right)-\ln \left(P\left(T\right)\right)+\ln \left(Z\right)$

去掉与活动轮廓曲线c无关的量，则最大化F(c)可以写成最小化如下能量函数形式：

$E_{\mathrm{POL}}(V_1,V_2,c)=\mathrm{ulength}\left(c\right)\·g(x,y)+{\mathrm{\λ}}_1{\∫}_{\mathrm{inside}\left(c\right)}\left(\log \right|V_1|+\mathrm{Tr}(V_1^{-1}\·T(x,y)))\mathrm{dxdy}+$

${\mathrm{\λ}}_2{\∫}_{\mathrm{outside}\left(c\right)}\left(\log \right|V_2|+\mathrm{Tr}(V_2^{-1}\·T(x,y)))\mathrm{dxdy}$

其中，c为活动轮廓曲线，u、λ₁和λ₂为设定常数，V₁和V₂分别是分割区域的均值，length(c)表示分割边界的长度，∫_inside(c)(·)dxdy和∫_outside(c)(·)dxdy分别表示对活动轮廓曲线c内和活动轮廓曲线c外所有的像素点进行积分，Tr(·)表示求矩阵的迹，(x，y)是图像中的像素点的坐标，T(x，y)表示像素点的散射相干矩阵。

步骤104：将所述分割能量函数嵌入到水平集函数φ(x，y)中，令λ＝λ₁＝λ₂，得到水平集分割能量泛函为

$E_{\mathrm{POL}}(V_1,V_2,\mathrm{\φ})=U{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)g(x,y)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}$

$+\mathrm{\λ}\left(\begin{array}{c}{\∫}_{\mathrm{\Ω}}\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+\\ {\∫}_{\mathrm{\Ω}}\left(\log \right|V_2|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}\end{array}\right)$

其中，Ω是待分割的图像区域，δ(·)表示Dirac函数，

表示求φ(x，y)的梯度，H(·)是Heaviside函数，当x＞1时，H(x)＝1；当x＜0时，H(x)＝0；当x＝0时，H(x)无对应的值。水平集函数φ(x，y)的定义如下：

${\mathrm{\φ}}_0(x,y)=\left\{\begin{array}{cc}d& (x,y)\∈{\mathrm{\Ω}}_0\\ 0& (x,y)\∈{\∂\mathrm{\Ω}}_0\\ -d& (x,y)\∈\mathrm{\Ω}-{\mathrm{\Ω}}_0\end{array}\right.$

Ω₀是活动轮廓曲线c内部的区域，d为大于0的常数。

步骤105：初始化水平集函数并设定水平集函数的发展方程的参数。

油膜对海面的后向衰减特性至少要有1.4dB，因此利用1.4dB的油膜衰减特性对水平集函数φ进行初始化，设定初始水平集函数为φ₀，且令

${\mathrm{\φ}}_0(x,y)=\left\{\begin{array}{cc}d& (x,y)\∈{\mathrm{\Ω}}_0\\ 0& (x,y)\∈{\∂\mathrm{\Ω}}_0\\ -d& (x,y)\∈\mathrm{\Ω}-{\mathrm{\Ω}}_0\end{array}\right.$

另外，设定水平集函数的发展方程的参数t＝0、λ＝10、k＝9，其中，t为迭代初值。由于水平集发展方程中的参数u和λ对最终方程演化结果的影响很大，为了取得较好的演化结果，与u和λ有关的统计量必须成比例，这里取经验值0.25，并用比值

对参数进行确定，具体是：设定λ的值后，选取一个均匀区域Ω_h，分别计算

和的值；依据

计算u的大小。

步骤106：计算每个分割区域的均值；对水平集分割的发展方程求解，先保持函数φ_t(x，y)不动，求关于V₁和V₂的最小能量函数E_POL(V₁，V₂，φ)。根据最大似然检测准则可以得到：

$V_1=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}$ 和 $V_2=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}}$

步骤107：利用公式

更新水平集函数；其中，t＝0，1，2...为迭代次数，dt为步长并设定dt＝0.5。

是根据水平集分割能量泛函和每个分割区域的均值，确定的水平集函数发展方程。即先保持V₁和V₂的值不变，根据Euler-Lagrange方程最小化关于φ_i的能量函数E_POL(V₁，V₂，φ_i)，结合负梯度下降法，得到水平集函数的发展方程如下

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\δ}\left(\mathrm{\φ}\right)(u\mathrm{div}(g\·\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|})-\mathrm{\λ}(\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T\right))-\left(\log \right|V_2|+\mathrm{Tr}\left(V_2^{-1}T\right))))$

其中，div(·)表示求向量的散度。上式是窄带偏微分演化方程，为了加快演化以及防止陷入演化局部阈值，将δ(φ)去掉用

代替，可以得到水平集函数φ的发展方程

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(u\mathrm{div}(g\·\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|})-\mathrm{\λ}(\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T\right))-\left(\log \right|V_2|+\mathrm{Tr}\left(V_2^{-1}T\right)))).$

步骤108：判断|φ_t+1-φ_t|＜ε是否成立，如果|φ_t+1-φ_t|＜ε成立，则执行步骤109；否则，令t＝t+1，返回步骤106；其中，ε为设定阈值。

步骤109：分割结束。

通过上述说明可以看出，本发明提供的方法是一种改进的统计活动轮廓法，它适用于SAR数据以及SAR的数据融合，实现了海面溢油图像的自动分割。

实施例2

下面以一个实际的海面溢油图像仿真数据说明本发明的实施过程。

步骤201：图2给出了仿真油膜数据的总功率图，图像的大小是571×421，海面和油膜的后向散射系数分别为

[|S_HH|² 2|S_HV|² |S_VV|²]_oil＝[0.0031 5.3598×10^-4 0.0047]

                                                        (1)

[|S_HH|² 2|S_HV|² |S_VV|²]_sea＝[0.0048 5.7513×10^-4 0.0093]

此仿真图像的油膜形状选自真实ASAR传感器获取的溢油形状，具体溢油形状图见图3所示。在本实施案例中，所应用的数据源分别为三组单极化SAR强度数据，|S_HH|²，|S_HV|²和|S_VV|²；两组双极化SAR观测模式的强度融合数据，|S_HH|²+|S_HV|²和|S_HV|²+|S_VV|²；以及一组全极化SAR观测模式的三个通道强度融合数据，|S_HH|²+|S_HV|²+|S_VV|²。

先以三个通道强度的融合数据为例，给出的具体实施例过程如下：

步骤202：由三个通道强度数据组成的对角矩阵T为，

$T=\left[\begin{array}{ccc}{\left|S_{\mathrm{HH}}\right|}^2& 0& 0\\ 0& {\left|S_{\mathrm{HV}}\right|}^2& 0\\ 0& 0& {\left|S_{\mathrm{VV}}\right|}^2\end{array}\right]---\left(2\right)$

T将不再服从复Wishart分布或者Gamma分布，但是对于复Wishart分布以及Gamma分布的统计测距公式仍旧是适用的。先考察它利用CFAR边缘检测算子得到的边缘特性R，设定虚警率为5％，并由虚警率根据假设检验理论计算出边缘阈值T_f。则边缘截止函数为，

$g=\frac{1+k}{1+k{(R/T_f)}^2}---\left(3\right)$

其中，k是常数。

步骤203：依据衰减比求解初始化的水平集函数，依据C波段星载SAR油膜的后向散射衰减比，在中低海况下，VV通道的衰减比在2.2dB-7dB之间，HH通道的衰减比在1.4dB-5dB之间。因此，设定衰减比阈值为1.4dB，初始化水平集函数为如下三段式模式，

${\mathrm{\φ}}_0=\left\{\begin{array}{cc}d& (x,y)\∈{\mathrm{\Ω}}_0\\ 0& (x,y)\∈{\∂\mathrm{\Ω}}_0\\ -d& (x,y)\∈\mathrm{\Ω}-{\mathrm{\Ω}}_0\end{array}\right.---\left(4\right)$

测试数据的初始化水平集函数以及活动轮廓曲线如图4所示。

步骤204：本实施方案中，除系数u之外，其余参数的设定如下：λ＝10，dt＝0.5，k＝9。选取海面一块均匀区域，统计(∑_Ωg)/λ∑_ΩTr(V^-1T)的值，由不同的测试数据得到的值是不同的。对于此三个通道融合的强度数据，计算得到的u＝0.7788。

步骤205：固定测试数据的初始化水平集函数φ₀，求取能量函数E_POL(V₁，V₂，φ₀)的关于区域统计量V₁和V₂的最小值，根据最大似然比检测准则，可以得到区域统计均值如下：

$V_1=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·H\left({\mathrm{\φ}}_0(x,y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}H\left({\mathrm{\φ}}_0(x,y)\right)\mathrm{dxdy}};$ $V_2=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·(1-H\left({\mathrm{\φ}}_0(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H\left({\mathrm{\φ}}_0(x,y)\right))\mathrm{dxdy}}---\left(5\right)$

步骤206：固定求取的V₁和V₂的值，利用Euler-Lagrange方程以及负梯度下降法求解水平集函数的发展方程，

${\mathrm{\φ}}_{t+1}={\mathrm{\φ}}_t+$

$\mathrm{dt}\×\left|{\▿\mathrm{\φ}}_t\right|(u\mathrm{div}(g\·\frac{{\▿\mathrm{\φ}}_t}{{\▿\mathrm{\φ}}_t})-\mathrm{\λ}(\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T\right))-\left(\log \right|V_2|+\mathrm{Tr}\left(V_2^{-1}T\right))))---\left(6\right)$

其中，t表示迭代次数，对于初始的情况，t＝0。

步骤207：判断|φ_t+1-φ_t|是否大于设定阈值ε，若|φ_t+1-φ_t|是否大于设定阈值ε，则以φ_t+1取代φ_t，重复步骤205-步骤207；否则，曲线演化过程结束。

利用此实施过程对三个通道强度融合数据进行分割，最终的水平集演化结果如图5所示。从图5中可以看出，最终能量泛函的零水平集函数与油膜区域的边缘很接近，并且水平集函数近似为符号距离函数，说明此发明所提出的目标分割模型有效。求取最终水平集演化结果的φ_n＞0的区域，得到最终海面暗区域的分割结果，见图6所示。

前述步骤介绍了发明的统计活动轮廓模型法在极化SAR海面溢油检测中的应用。为了进一步说明融合数据在SAR海面溢油检测中的性能以及此发明的有效性，下面分别考察单极化以及双极化数据源的情况。

单极化数据的分割结果如图7所示。从图中可以看出VV极化方式以及HH极化方式可以将油膜区域分割出来，而HV通道的数据对油膜分割的结果完全失效。HV通道的数据中，油膜和海面之间的后向衰减比仅为0.3dB，根据SAR图像辐射分辨率的定义可知，对于此SAR数据，为使两类目标能够被分开则要求目标之间的对比度至少为1.208dB。因此，此HV通道的数据对海面暗区域的分割失效。

双极化强度融合数据的分割结果见图8，两种双极化SAR强度融合数据均能够将油膜区域分割出来，但是HV通道+VV通道的融合数据分割效果较好。本发明的一个主要特点是此活动轮廓模型可以适用于各种SAR数据中。因此，可以利用此模型评价不同的SAR数据在海面溢油分割中的性能，见表1所示。从表1的数据对比中可以看出，在海面溢油的监测任务中，多极化方式的检测概率较高；若仅利用一种极化方式观测，VV极化方式的SAR数据对海面油膜检测效果较好。

表1利用本发明方法的极化SAR融合数据在海面溢油检测中的性能比较

为说明本发明方法在溢油分割中的有效性，下面将基于Gamma分布的监督最大似然(ML)分类方法与此发明方法做对比。监督ML分类器的结果见图9所示，图中的两个方框区域为训练样本区域。设一共有w₁，...，w_n类目标，各类目标的均值为V₁，...，V_m。则监督的ML分类器判别准则为，

如果d(T_i，V_m)≤d(T_i，V_k)，w_m≠w_k，则T_i∈w_m.                (7)

从图9以及图6的对比中可以看出，本发明方法在分割区域的完整性以及连续性，分割边界的完整性以及连续性等方面的保持效果都优于监督的ML分类方法。

通过上述不同数据源以及方法的实施案例对比，可以看出本发明提出的方法相较于目前现有的海面溢油检测算法，其特点主要体现在：

1、同时应用了海面溢油的衰减特性、统计特性、以及空间纹理特性；

2、是一种海面暗区域自动检测分割模型；

3、证明了本发明可以应用于任何极化SAR数据以及极化SAR融合数据中，是极化SAR数据基于统计活动轮廓法的统一最小泛函分割模型；

4、通过给出的经验比值，本发明可以自动确定活动轮廓模型中的演化参数。

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

Claims (3)

1.一种基于极化SAR融合数据的海面溢油分割方法”，其特征是所述方法包括：

步骤1：选取不同合成孔径雷达数据源对海面溢油图像的监测数据；

步骤2：根据选取的监测数据确定所述图像的散射相干矩阵；

步骤3：依据活动轮廓法对所述图像进行分割。

2.根据权利要求1所述的方法，其特征是所述散射相干矩阵为分块对角矩阵，具体形式为

其中，C_i为1×1、2×2或3×3的半正定Hermitian矩阵，i＝1，2，...，n，n为半正定Hermitian矩阵的个数，并且T满足P(T|V)＝P₁(C₁|V₁)P₂(C₂|V₂)...P_n(C_n|V_n)，P(·)是矩阵T的概率密度分布函数，V＝E(T)是散射相干矩阵T的均值，P_i(·)是矩阵C_i的概率密度分布函数，V_i＝E(C_i)是矩阵C_i的均值。

3.根据权利要求1所述的方法，其特征是所述步骤3具体是：

步骤101：利用活动轮廓曲线c将所述图像分割成2个区域；

步骤102：利用公式g＝(1+k)/(1+k(R/T_f)²)计算每个分割区域的边缘截止函数；其中，g是边缘截止函数，R是由恒虚警率边缘检测算子得到的所述图像边缘密度，T_f由虚警率得到的边缘检测阈值，k是设定常数；

步骤103：根据公式

$E_{\mathrm{POL}}(V_1,V_2,c)=\mathrm{ulength}\left(c\right)\·g(x,y)+{\mathrm{\λ}}_1{\∫}_{\mathrm{inside}\left(c\right)}\left(\log \right|V_1|+\mathrm{Tr}(V_1^{-1}\·T(x,y)))\mathrm{dxdy}+$

${\mathrm{\λ}}_2{\∫}_{\mathrm{outside}\left(c\right)}\left(\log \right|V_2|+\mathrm{Tr}(V_2^{-1}\·T(x,y)))\mathrm{dxdy}$

计算所述图像分割后的分割能量函数；其中，c为活动轮廓曲线，u、λ₁和λ₂为设定常数，V₁和V₂分别是分割区域的均值，length(c)表示分割边界的长度，∫_inside(c)(·)dxdy和∫_outside(c)(·)dxdy分别表示对活动轮廓曲线c内和活动轮廓曲线c外所有的像素点进行积分，Tr(·)表示求矩阵的迹，(x，y)是图像中的像素点的坐标，T(x，y)表示像素点的散射相干矩阵；

步骤104：将所述分割能量函数嵌入到水平集函数φ(x，y)中，令λ＝λ₁＝λ₂，得到水平集分割能量泛函为

$E_{\mathrm{POL}}(V_1,V_2,\mathrm{\φ})=u{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)g(x,y)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}$

$+\mathrm{\λ}\left(\begin{array}{c}{\∫}_{\mathrm{\Ω}}\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+\\ {\∫}_{\mathrm{\Ω}}\left(\log \right|V_2|+\mathrm{Tr}\left(V_1^{-1}T(x,y)\right))(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}\end{array}\right)$

其中，Ω是待分割的图像区域，δ(·)表示Dirac函数，

表示求φ(x，y)的梯度，H(·)是Heaviside函数；

步骤105：初始化水平集函数并设定水平集函数的发展方程的参数；

所述初始化水平集函数具体是设定 ${\mathrm{\φ}}_0(x,y)=\left\{\begin{array}{cc}d& (x,y)\∈{\mathrm{\Ω}}_0\\ 0& (x,y)\∈{\∂\mathrm{\Ω}}_0\\ -d& (x,y)\∈\mathrm{\Ω}-{\mathrm{\Ω}}_0\end{array}\right.;$

所述设定水平集函数的发展方程的参数具体是设定t＝0、λ＝10、k＝9并且利用公式

计算u的值，Ω_h为选定的均匀区域；其中，t为迭代初值；

步骤106：计算每个分割区域的均值；其计算公式分别为

$V_1=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}H\left({\mathrm{\φ}}_t(x,y)\right)\mathrm{dxdy}}$ 和 $V_2=\frac{{\∫}_{\mathrm{\Ω}}T(x,y)\·(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}}{{\∫}_{\mathrm{\Ω}}(1-H\left({\mathrm{\φ}}_t(x,y)\right))\mathrm{dxdy}};$

步骤107：利用公式

${\mathrm{\φ}}_{t+1}={\mathrm{\φ}}_t+\mathrm{dt}\×\frac{\∂\mathrm{\φ}}{\∂t}$

更新水平集函数；其中，t＝0，1，2...为迭代次数，dt为步长并设定dt＝0.5；为水平集函数的发展方程

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(u\mathrm{div}(g\·\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|})-\mathrm{\λ}(\left(\log \right|V_1|+\mathrm{Tr}\left(V_1^{-1}T\right))-\left(\log \right|V_2|+\mathrm{Tr}\left(V_2^{-1}T\right))))$

其中，div(·)表示求向量的散度；

步骤108：判断|φ_t+1-φ_t|＜ε是否成立，如果|φ_t+1-φ_t|＜ε成立，则执行步骤109；否则，令t＝t+1，返回步骤106；其中，ε为设定阈值；

步骤109：分割结束。

Images