CN102880856A - 基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法 - Google Patents

Abstract

本发明公开了一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法。步骤如下：1）对航空高光谱遥感影像的每个像元计算归一化表观反射率；2）对归一化表观反射率影像的每个像元计算油膜指数；3）基于油膜指数影像构建油膜指数梯度影像和油膜指数梯度角度影像；4）对边缘梯度值进行非极大值抑制；5）利用非极大值抑制后的油膜指数梯度影像形成油膜边缘影像；6）对油膜指数影像的每个像元识别薄油膜，形成薄油膜分布影像；7）在油膜边缘影像、薄油膜分布影像的基础上提取海表溢油全分布遥感影像。本发明较好的实现了航空高光谱遥感海面溢油的自动化识别，可用于指导海洋溢油的清理评估工作，对海洋环保具有重要意义。

Description

基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法

技术领域

本发明涉及遥感技术领域，尤其涉及一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法。

背景技术

随着全球工业化进程的快速推进，人类社会对石油的需求量与日俱增，其开采规模迅速扩大，海上油井数量和海上石油运输量急剧增加。由海上作业、运输和设备等原因造成的海面溢油事件也在不断增加，这加剧了海洋石油污染，进而危害了海洋生物的生存环境，严重威胁着海洋生态环境和渔业资源。相比船测，利用遥感手段监测溢油信息具有实时、全面、低成本优势，能有效地指导清理、评估工作，降低溢油对海洋生态环境的破坏。因此，研究溢油遥感监测机理和技术对海洋生态环境保护具有重要意义。

张永宁^[1]在黄海洁净海域选用重柴油、轻柴油、润滑油、煤油、中东原油和大庆原油进行了油膜光谱实验。通过对实验数据的细致分析，为TM影像的遥感解译提供了波段特征。赵冬至^[2]通过油膜光谱实验，分析了原油、柴油及润滑油的反射率峰值位置及反射率峰值随厚度的变化情况，并将其应用于卫星遥感溢油监测。陆应诚^[3]使用辽东湾海水样本与辽河油田原油，在400-2500nm波段范围内对多种厚度的油膜光谱特性进行测量。结果表明，在400-1150nm波段范围内，油膜光谱反射率与厚度呈很强的幂函数负相关关系。这些都是对溢油监测基础的研究，在卫星遥感上也进行了应用，但未在航空遥感方面开展研究。

通过对此类文献的分析总结可知，当前绝大多数提取算法都针对已知有溢油分布的影像，未实现有效的自动化提取。目前，光学影像的溢油信息提取方法主要分为以下两类：

1．基于统计的图像处理方法

此类方法主要分两种。

a．统计分类法

基于统计分类的方法有可分为两个类别，一种为非监督分类，不需要先验样本，但分类后需要人工判定类别，确定是否为溢油。如Sanchez^[4]采用ISODATA聚类提取油膜信息。另一种为监督分类，主要有如下几个步骤：提取训练样本，选择分类判据，确定判别准则。样本提取的主要方法有：目视解译^[5]、DancingPixel^[6]、光谱曲线库等。分类判据主要采用光谱角匹配^[5,6,7]，判别准则采用最大似然法^[5]、BP、GANN^[8]等。

b．纹理识别法

计算纹理特征，结合阈值分割提取油膜信息^[9,10]。这类提取方法主要应用于单波段影像或作为辅助方法用于多波段影像的油膜信息提取。

基于统计法的溢油监测需要符合以下条件：影像中各类别光谱差异较明显，溢油区与周边水体边界较窄。而大多数遥感影像中，较难使用单一或少数离散光谱代表不同水体信息，如存在耀斑、水体光谱连续变化等，因此，该方法局限性较大。

2．彩色合成法

对各波段进行直方图拉伸、组合运算、主成分分析^[11-14]，使油、水差异最大化，最终将差异最大的三个结果进行彩色合成，目视判别油膜分布，或人工确定阈值进行分割。

基于光谱特征差异的海表溢油遥感监测方法目前较少，而且尚未形成成熟的自动监测方法。而本专利通过相关的前期研究，针对油、水光谱差异特征（幅值、谱型），提出了一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法。

参考文献：

[1]张永宁,丁倩,李栖筠.海上溢油污染遥感监测的研究.大连海事大学学报.1999.25(3).1-5

[2]赵冬至,丛丕福.油膜波谱特征分析与遥感监测溢油.遥感技术与应用.2000.15(3).160-164.

[3]陆应诚,田庆久,王晶晶等.海面油膜光谱响应实验研究.科学通报.2008.53(9).1085-1088.

[4]Sanchez G.,Roper W.E.,Gomez R..Detection and monitoring of oil spills usinghyperspectral imagery.Proc.of SPIE.2003.5097.233-240

[5]Salem F.,Kafatos M..Hyperspectral Image Analysis for Oil Spill Mitigation,http://www.crisp.nus.edu.sg.

[6]Salem F.,Kafatos M.,El-Ghazawi K.,etc..Hyperspectral image assessment ofoil-contaminated wetland.International Journal of Remote Sensing.2005.26(4).811-821.

[7]卜志国.航空高光谱遥感赤潮与油污染水体信息提取.中国海洋大学硕士学位论文.2004.

[8]F.

and R.J.Duro.A hyperspectral based multisensor system formarine oil spill detection.Intelligent Information and Engineering Systems.2004.3213.669-676.

[9]石立坚.SAR及MODIS数据海面溢油监测方法研究.中国海洋大学博士学位论文.2008.

[10]Shi L.,Zhang X.,Seielstad G.,etc.;Oil Spill Detection by MODIS Images usingFuzzy Cluster and Texture Feature Extraction.OCEANS 2007–Europe.2007.1567-1571.

[11]李栖筠,陈维英,肖乾广.老铁山水道漏油事故卫星监测.环境遥感,1994,9(4):256～262.

[12]丁倩.海洋溢油卫星遥感图像处理.大连海事大学硕士学位论文.2000.

[13]杨娜.NOAA_AVHRR数据在事故溢油中应用研究.中国海洋大学硕士毕业论文.2006.

[14]靳伟.基于波谱特征的溢油识别和监测.大连海事大学硕士学位论文.2008.

发明内容

本发明的目的是克服现有技术的不足，提供一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法。

基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法步骤如下：

1）对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像；

2）对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像；

3）基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像；

4）基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像；

5）针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像；

6）对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像；

7）在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像。

所述的步骤1）为：对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像：为了减少不同光照条件引起的辐亮度差异，利用高光谱遥感器自带的下行辐亮度光纤传感器FODIS，Fiber opticdownwelling irradiance sensor，测量的向下光谱辐照度，在此基础上对航空高光谱遥感影像的每个像元点计算每个波段的表观反射率T_rs，计算公式如下：

$T_{\mathrm{rs}}\left(\mathrm{\λ}\right)=\frac{\mathrm{Radiance}\left(\mathrm{\λ}\right)}{\mathrm{FODIS}\left(\mathrm{\λ}\right)}---\left(1\right)$

上式中Radiance（λ）表示像元点高光谱遥感器中心波长为λ的波段的入瞳辐射，FODIS（λ）表示同一像元点下行辐亮度光纤传感器测得的中心波长为λ的波段的太阳辐射，这样就形成表观反射率影像；

然后基于上述表观反射率影像对每个像元点计算每个波段的归一化表观反射率T_normrs，计算公式如下：

$T_{\mathrm{normrs}}\left(\mathrm{\λ}\right)=\frac{T_{\mathrm{rs}}\left(\mathrm{\λ}\right)*\mathrm{FWHM}\left(\mathrm{\λ}\right)}{\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}{T}_{\mathrm{rs}}\left({\mathrm{\λ}}_i\right)*\mathrm{FWHM}\left({\mathrm{\λ}}_i\right)}---\left(2\right)$

上式中T_rs(λ)是像元点中心波长为λ的波段的表观反射率，FWHM（λ）是同一像元点中心波长为λ的波段的半波长宽度，

代表对高光谱数据在该点所有波段的T_rs(λ_i)*FWHM(λ_i)进行求和积分，i=1，2，3，…，n，n为高光谱遥感器的波段数，这样对每个像元点进行计算后就形成了航空高光谱归一化表观反射率影像。

所述的步骤2）为：对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像：对航空高光谱归一化表观反射率影像的每个像元点计算油膜指数I_oil，公式如下：

$I_{\mathrm{oil}}=\frac{T_{\mathrm{normrs}}\left(894\right)}{T_{\mathrm{normrs}}\left(516\right)}---\left(3\right)$

上式中T_normrs(894)为像元点中心波长为894nm的波段的归一化表观反射率，T_normrs(516)为同一像元点中心波长为516nm的波段的归一化表观反射率，这样对每个点计算后，就形成油膜指数影像；然后采用高斯滤波模板对油膜指数影像进行高斯滤波，滤波方法采用3×3高斯滤波模板在油膜指数影像上移动计算实现，滤波模板为 $\frac{1}{16}\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right],$ 具体计算方法如下：在油膜指数影像中逐点移动高斯滤波模板，通过计算高斯滤波模板系数与油膜指数影像对应位置的像素值的乘积之和确定油膜指数影像滤波后的响应值，油膜指数影像在像元点(x,y)处响应值R计算公式为：

R=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+...+w(0,0)f(x,y)+...

+w(1,0)f(x+1,y)+w(1,1)f(x+1,y+1)                              （4）

式中，w(i，j)为高斯滤波模板在(i，j)处的取值，i=-1，0，1，j＝-1，0，1，w(0,0)与f(x,y)重合，f(x,y)为油膜指数影像在像元点(x,y)处的取值；这样就形成了滤波后的航空高光谱油膜指数影像。

所述的步骤3）为：基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像：本步骤采用Sobel边缘检测方法提取梯度信息，这里的Sobel边缘检测方法使用两个有向算子，其中一个水平有向算子 $s_1=\left(\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right),$ 一个垂直有向算子 $s_2=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right),$ 分别用水平有向算子和垂直有向算子对航空高光谱油膜指数影像进行卷积，得到两个矩阵G_x(x，y)、G_y(x,y)，他们分别表示航空高光谱油膜指数影像的偏导数矩阵，在不考虑边界的情形下也是和航空高光谱油膜指数影像同样大小，然后把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数平方后相加并开根号得到一个新的矩阵G(x,y)，G(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度分布，把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数相除并进行反正切运算得到一个新的矩阵θ(x,y)，θ(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度角度分布，这样就得到了航空高光谱油膜指数梯度影像和油膜指数梯度角度影像；

所述的步骤4）为：基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像：油膜指数梯度影像的像元值越大说明其为边缘可能性也越大，但为了精确定位边缘，必须细化可能的边缘带，保留梯度值局部最大的点，这就是边缘梯度非极大值抑制过程，本步骤在3×3窗口内进行非极大值抑制，首先将油膜指数梯度角度影像中各像元的梯度角度分为4类，其中，-22.5°~22.5°、157.5°~180°、-180°~-157.5°设为一类编号为0，22.5°~67.5°、-157.5°~-112.5°设为一类编号为1，67.5°~112.5°、-112.5°~-67.5°设为一类编号为2，-67.5°~-22.5°、112.5°~157.5°设为一类编号为3。这4类分别对应横、竖及两个斜对角方向。根据角度类别，对应的油膜指数梯度影像中每一个像元与相应方向上的两个像元进行梯度比较，若其梯度不是最大，则将像元梯度值置为0，从而完成对边缘梯度值的非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像。

所述的步骤5）为：针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像：通过高T_high、低阈值T_low对上述非极大值抑制结果进行边缘提取，具体判断规则公如下：

$\left\{\begin{array}{cc}\left|G(x,y)\right|\≥T_{\mathrm{high}}& \→I_1(x,y)=1\\ T_{\mathrm{high}}\≥\left|G(x,y)\right|\≥T_{\mathrm{low}}& \→I_2(x,y)=1\end{array}\right.---\left(5\right)$

式中G(x,y)表示航空高光谱油膜指数梯度影像，x、y分别代表横纵坐标，T_high为高阈值，T_low为低阈值，I₁、I₂为边缘标志图像，初始置0，置1表示为初始边缘，由于I₁图像采用高阈值提取，故其标识的边缘确定为所需边缘，但最终提取结果还将包括I₂图像中标识的部分或全部边缘，本步骤通过梯度值图像的直方图统计，自适应获取最大梯度值|G|_max，根据这一参数构建高、低阈值，其计算公式如下：

$\left\{\begin{array}{c}{T}_{\mathrm{high}}=0.9\×{\left|G\right|}_{\max }\\ T_{\mathrm{low}}=0.85\×{\left|G\right|}_{\max }\end{array}\right.---\left(6\right)$

然后进行边缘连接，I₁标识的边缘存在间断，需要进行连接修补，具体步骤为：a．扫描I₁图像，当I₁(x₀,y₀)＝1时，在I₁中跟踪起始点为S(x₀,y₀)的边缘，直到该边缘的终点E(x₁,y₁)；b．查看I₂中与E点对应的E′点周围8邻域内是否有点R'(x',y')，使得I₂(x',y')＝1，若有，则将I₁中与R′对应的R点置为1，即I₁(x',y′)＝1；再重复a步骤，继续查找最终边缘点，直到无法从I₂中提取最终边缘点，从而完成一条边缘检测，并标记为已检测；c．当完成一条边缘检测后继续在I₁中查找新的边缘起始点，并重复a、b步骤。最终完成所有边缘检测，形成航空高光谱遥感油膜边缘影像；虽然这个步骤提取了大部分油膜边缘，但利用这些边缘直接提取油膜分布难度较大。海面波纹、亮度渐变使得提取的油膜边缘发生间断，部分地区油膜边缘因过于细小而被忽略，因此需要结合阈值法提取的薄油膜分布遥感影像获得全油膜分布，参见下一步骤。

所述的步骤6）为：对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像：通过阈值法识别薄油膜分布，具体判断规则公式如下：

上式中Thred为判断阈值，I_oil(x,y)为航空高光谱油膜指数影像中点(x,y)处的值；我们的方法采用0.75这一阈值提取薄油膜分布，形成二值影像；然后采用窗口统计方式对薄油膜分布进行滤波，其具体方法为：在薄油膜分布图中，遍历每个像元，每次遍历采用如下判据判断像元是否为真正的薄油膜像元：当像元不为薄油膜时，则维持原来的判断；当像元为薄油膜时，以该像元为中心统计9×9邻域窗口内的薄油膜像元数，若总数超过9个，则将其判断为薄油膜，反之，则不为薄油膜；这样就形成航空高光谱薄油膜分布遥感影像。

所述的步骤7）为：在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像：首先，将油膜边缘叠加到薄油膜分布图上，然后通过四连通方式对非薄油膜像元进行连通标记，由于清洁水体像元基本是大面积连续分布，厚油膜分布在薄油膜内部，其斑块面积远小于清洁水体，因此，可将连通区域大小在阈值300以下的区域判为厚油膜，最后综合厚、薄油膜分布得到海表溢油全分布遥感影像。

本发明与现有技术相比具有的有益效果：通过航空高光谱遥感影像中油、水的光谱差异特征，构建了油膜指数，并结合阈值分割、边缘提取、连通分析等方法，提取了油膜分布，充分发挥航空高光谱在海洋溢油监测中的作用，将为溢油事件的清理、评估工作提供及时、全面、准确的科学依据。

附图说明

图1是基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法技术路线图；

图2是原始的航空高光谱影像；

图3是航空高光谱油膜指数影像；

图4是基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像对边缘梯度值进行非极大值抑制时，油膜指数梯度角度影像中各像元分类的示意图；

图5是航空高光谱遥感油膜边缘影像；

图6是薄油膜分布遥感影像，左侧为缩略图、右上方为原始分辨率图像、右下方为16倍局部放大图像；

图7是海表溢油全分布遥感影像。

具体实施方式

基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法流程参见附图1，步骤如下：

1）对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像；

2）对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像；

3）基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像；

4）基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像；

5）针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像；

6）对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像；

7）在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像。

所述的步骤1）为：对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像：为了减少不同光照条件引起的辐亮度差异，利用高光谱遥感器自带的下行辐亮度光纤传感器FODIS，Fiber opticdownwelling irradiance sensor，测量的向下光谱辐照度，在此基础上对航空高光谱遥感影像的每个像元点计算每个波段的表观反射率T_rs，计算公式如下：

$T_{\mathrm{rs}}\left(\mathrm{\λ}\right)=\frac{\mathrm{Radiance}\left(\mathrm{\λ}\right)}{\mathrm{FODIS}\left(\mathrm{\λ}\right)}---\left(1\right)$

上式中Radiance（λ）表示像元点高光谱遥感器中心波长为λ的波段的入瞳辐射，FODIS（λ）表示同一像元点下行辐亮度光纤传感器测得的中心波长为λ的波段的太阳辐射，这样就形成表观反射率影像；

然后基于上述表观反射率影像对每个像元点计算每个波段的归一化表观反射率T_normrs，计算公式如下：

$T_{\mathrm{normrs}}\left(\mathrm{\λ}\right)=\frac{T_{\mathrm{rs}}\left(\mathrm{\λ}\right)*\mathrm{FWHM}\left(\mathrm{\λ}\right)}{\underset{{\mathrm{\λ}}_{\min }}{\overset{{\mathrm{\λ}}_{\max }}{\∫}}{T}_{\mathrm{rs}}\left({\mathrm{\λ}}_i\right)*\mathrm{FWHM}\left({\mathrm{\λ}}_i\right)}---\left(2\right)$

上式中T_rs(λ)是像元点中心波长为λ的波段的表观反射率，FWHM（λ）是同一像元点中心波长为λ的波段的半波长宽度，

代表对高光谱数据在该点所有波段的T_rs(λ_i)*FWHM(λ_i)进行求和积分，i=1，2，3，…，n，n为高光谱遥感器的波段数，这样对每个像元点进行计算后就形成了航空高光谱归一化表观反射率影像。

所述的步骤2）为：对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像：对航空高光谱归一化表观反射率影像的每个像元点计算油膜指数I_oil，公式如下：

$I_{\mathrm{oil}}=\frac{T_{\mathrm{normrs}}\left(894\right)}{T_{\mathrm{normrs}}\left(516\right)}---\left(3\right)$

上式中T_normrs(894)为像元点中心波长为894nm的波段的归一化表观反射率，T_normrs(516)为同一像元点中心波长为516nm的波段的归一化表观反射率，这样对每个点计算后，就形成油膜指数影像；然后采用高斯滤波模板对油膜指数影像进行高斯滤波，滤波方法采用3×3高斯滤波模板在油膜指数影像上移动计算实现，滤波模板为 $\frac{1}{16}\left[\begin{array}{ccc}1& 2& 1\\ 2& 4& 2\\ 1& 2& 1\end{array}\right],$ 具体计算方法如下：在油膜指数影像中逐点移动高斯滤波模板，通过计算高斯滤波模板系数与油膜指数影像对应位置的像素值的乘积之和确定油膜指数影像滤波后的响应值，油膜指数影像在像元点(x,y)处响应值R计算公式为：

R=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+...+w(0,0)f(x,y)+...

+w(1,0)f(x+1,y)+w(1,1)f(x+1,y+1)                              （4）

式中，w(i，j)为高斯滤波模板在(i，j)处的取值，i=-1，0，1，j＝-1，0，1，w(0,0)与f(x,y)重合，f(x,y)为油膜指数影像在像元点(x,y)处的取值；这样就形成了滤波后的航空高光谱油膜指数影像，参见附图3。

所述的步骤3）为：基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像：本步骤采用Sobel边缘检测方法提取梯度信息，这里的Sobel边缘检测方法使用两个有向算子，其中一个水平有向算子 $s_1=\left(\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right),$ 一个垂直有向算子 $s_2=\left(\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right),$ 分别用水平有向算子和垂直有向算子对航空高光谱油膜指数影像进行卷积，得到两个矩阵G_x(x，y)、G_y(x,y)，他们分别表示航空高光谱油膜指数影像的偏导数矩阵，在不考虑边界的情形下也是和航空高光谱油膜指数影像同样大小，然后把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数平方后相加并开根号得到一个新的矩阵G(x,y)，G(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度分布，把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数相除并进行反正切运算得到一个新的矩阵θ(x,y)，θ(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度角度分布，这样就得到了航空高光谱油膜指数梯度影像和油膜指数梯度角度影像；

所述的步骤4）为：基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像：油膜指数梯度影像的像元值越大说明其为边缘可能性也越大，但为了精确定位边缘，必须细化可能的边缘带，保留梯度值局部最大的点，这就是边缘梯度非极大值抑制过程，本步骤在3×3窗口内进行非极大值抑制，首先将油膜指数梯度角度影像中各像元的梯度角度分为4类，其中，-22.5°~22.5°、157.5°~180°、-180°~-157.5°设为一类编号为0，22.5°~67.5°、-157.5°~-112.5°设为一类编号为1，67.5°~112.5°、-112.5°~-67.5°设为一类编号为2，-67.5°~-22.5°、112.5°~157.5°设为一类编号为3，参见附图4。这4类分别对应横、竖及两个斜对角方向。根据角度类别，对应的油膜指数梯度影像中每一个像元与相应方向上的两个像元进行梯度比较，若其梯度不是最大，则将像元梯度值置为0，从而完成对边缘梯度值的非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像。

所述的步骤5）为：针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像：通过高T_high、低阈值T_low对上述非极大值抑制结果进行边缘提取，具体判断规则公式如下：

$\left\{\begin{array}{cc}\left|G(x,y)\right|\≥T_{\mathrm{high}}& \→I_1(x,y)=1\\ T_{\mathrm{high}}\≥\left|G(x,y)\right|\≥T_{\mathrm{low}}& \→I_2(x,y)=1\end{array}\right.---\left(5\right)$

式中G(x,y)表示航空高光谱油膜指数梯度影像，x、y分别代表横纵坐标，T_high为高阈值，T_low为低阈值，I₁、I₂为边缘标志图像，初始置0，置1表示为初始边缘，由于I₁图像采用高阈值提取，故其标识的边缘确定为所需边缘，但最终提取结果还将包括I₂图像中标识的部分或全部边缘，本步骤通过梯度值图像的直方图统计，自适应获取最大梯度值|G|_max，根据这一参数构建高、低阈值，其计算公式如下：

$\left\{\begin{array}{c}{T}_{\mathrm{high}}=0.9\×{\left|G\right|}_{\max }\\ T_{\mathrm{low}}=0.85\×{\left|G\right|}_{\max }\end{array}\right.---\left(6\right)$

然后进行边缘连接，I₁标识的边缘存在间断，需要进行连接修补，具体步骤为：a．扫描I₁图像，当I₁(x₀,y₀)＝1时，在I₁中跟踪起始点为S(x₀,y₀)的边缘，直到该边缘的终点E(x₁,y₁)；b．查看I₂中与E点对应的E′点周围8邻域内是否有点R'(x',y')，使得I₂(x',y')＝1，若有，则将I₁中与R′对应的R点置为1，即I₁(x',y′)＝1；再重复a步骤，继续查找最终边缘点，直到无法从I₂中提取最终边缘点，从而完成一条边缘检测，并标记为已检测；c．当完成一条边缘检测后继续在I₁中查找新的边缘起始点，并重复a、b步骤。最终完成所有边缘检测，形成航空高光谱遥感油膜边缘影像；虽然这个步骤提取了大部分油膜边缘，但利用这些边缘直接提取油膜分布难度较大。海面波纹、亮度渐变使得提取的油膜边缘发生间断，部分地区油膜边缘因过于细小而被忽略，因此需要结合阈值法提取的薄油膜分布遥感影像获得全油膜分布，参见下一步骤。参见附图5。

所述的步骤6）为：对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像：通过阈值法识别薄油膜分布，具体判断规则如下：

上式中Thred为判断阈值，I_oil(x,y)为航空高光谱油膜指数影像中点(x,y)处的值；我们的方法采用0.75这一阈值提取薄油膜分布，形成二值影像；然后采用窗口统计方式对薄油膜分布进行滤波，其具体方法为：在薄油膜分布图中，遍历每个像元，每次遍历采用如下判据判断像元是否为真正的薄油膜像元：当像元不为薄油膜时，则维持原来的判断；当像元为薄油膜时，以该像元为中心统计9×9邻域窗口内的薄油膜像元数，若总数超过9个，则将其判断为薄油膜，反之，则不为薄油膜；这样就形成航空高光谱薄油膜分布遥感影像，参见附图6。

所述的步骤7）为：在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像：首先，将油膜边缘叠加到薄油膜分布图上，然后通过四连通方式对非薄油膜像元进行连通标记，由于清洁水体像元基本是大面积连续分布，厚油膜分布在薄油膜内部，其斑块面积远小于清洁水体，因此，可将连通区域大小在阈值300以下的区域判为厚油膜，最后综合厚、薄油膜分布得到海表溢油全分布遥感影像，参见附图7。

Claims (8)

1.一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于它的步骤如下：

1）对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像；

2）对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像；

3）基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像；

4）基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像；

5）针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像；

6）对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像；

7）在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像。

2.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤1）为：对航空高光谱遥感影像的每个像元计算归一化表观反射率，形成航空高光谱归一化表观反射率影像：利用高光谱遥感器自带的下行辐亮度光纤传感器FODIS，Fiber optic downwelling irradiance sensor，测量的向下光谱辐照度，在此基础上对航空高光谱遥感影像的每个像元点计算每个波段的表观反射率T_rs，计算公式如下：

上式中Radiance（λ）表示像元点高光谱遥感器中心波长为λ的波段的入瞳辐射，FODIS（λ）表示同一像元点下行辐亮度光纤传感器测得的中心波长为λ的波段的太阳辐射，这样就形成表观反射率影像；

然后基于上述表观反射率影像对每个像元点计算每个波段的归一化表观反射率T_normrs，计算公式如下：

上式中T_rs(λ)是像元点中心波长为λ的波段的表观反射率，FWHM（λ）是同一像元点中心波长为λ的波段的半波长宽度， 

代表对高光谱数据在该点所有波段的T_rs(λ_i)*FWHM(λ_i)进行求和积分，i=1，2，3，…，n，n为高光谱遥感器的波段数，这样对每个像元点进行计算后就形成了航空高光谱归一化表观反射率影像。

3.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤2）为：对航空高光谱归一化表观反射率影像的每个像元计算油膜指数，形成航空高光谱油膜指数影像：对航空高光谱归一化表观反射率影像的每个像元点计算油膜指数I_oil，公式如下：

上式中T_normrs(894)为像元点中心波长为894nm的波段的归一化表观反射率，T_normrs(516)为同一像元点中心波长为516nm的波段的归一化表观反射率，这样对每个点计算后，就形成油膜指数影像；然后采用高斯滤波模板对油膜指数影像进行高斯滤波，滤波方法采用3×3高斯滤波模板在油膜指数影像上移动计算实现，滤波模板为具体计算方法如下：在油膜指数影像中逐点移动高斯滤波模板，通过计算高斯滤波模板系数与油膜指数影像对应位置的像素值的乘积之和确定油膜指数影像滤波后的响应值，油膜指数影像在像元点(x,y)处响应值R计算公式为：

R=w(-1,-1)f(x-1,y-1)+w(-1,0)f(x-1,y)+...+w(0,0)f(x,y)+...

+w(1,0)f(x+1,y)+w(1,1)f(x+1,y+1)                              （4）

式中，w(i，j)为高斯滤波模板在(i，j)处的取值，i=-1，0，1，j＝-1，0，1，w(0,0)与f(x,y)重合，f(x,y)为油膜指数影像在像元点(x,y)处的取值；这样就形成了滤波后的航空高光谱油膜指数影像。 

4.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤3）为：基于航空高光谱油膜指数影像构建航空高光谱油膜指数梯度影像和航空高光谱油膜指数梯度角度影像：本步骤采用Sobel边缘检测方法提取梯度信息，这里的Sobel边缘检测方法使用两个有向算子，其中一个水平有向算子

一个垂直有向算子

分别用水平有向算子和垂直有向算子对航空高光谱油膜指数影像进行卷积，得到两个矩阵G_x(x,y)、G_y(x,y)，他们分别表示航空高光谱油膜指数影像的偏导数矩阵，在不考虑边界的情形下也是和航空高光谱油膜指数影像同样大小，然后把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数平方后相加并开根号得到一个新的矩阵G(x,y)，G(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度分布，把矩阵G_x(x,y)、G_y(x,y)对应位置的两个数相除并进行反正切运算得到一个新的矩阵θ(x,y)，θ(x,y)表示航空高光谱油膜指数影像中各个像素的灰度的梯度角度分布，这样就得到了航空高光谱油膜指数梯度影像和油膜指数梯度角度影像。

5.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤4）为：基于航空高光谱油膜指数梯度影像和油膜指数梯度角度影像，对边缘梯度值进行非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像：本步骤在3×3窗口内进行非极大值抑制，首先将油膜指数梯度角度影像中各像元的梯度角度分为4类，其中，-22.5°~22.5°、157.5°~180°、-180°~-157.5°设为一类编号为0，22.5°~67.5°、-157.5°~-112.5°设为一类编号为1，67.5°~112.5°、-112.5°~-67.5°设为一类编号为2，-67.5°~-22.5°、112.5°~157.5°设为一类编号为3。这4类分别对应横、竖及两个斜对角方向。根据角度类别，对应的油膜 指数梯度影像中每一个像元与相应方向上的两个像元进行梯度比较，若其梯度不是最大，则将像元梯度值置为0，从而完成对边缘梯度值的非极大值抑制，形成非极大值抑制后的航空高光谱油膜指数梯度影像。

6.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤5）为：针对边缘梯度值进行非极大值抑制后的航空高光谱油膜指数梯度影像，通过双阈值方法进行边缘提取和连接，形成航空高光谱遥感油膜边缘影像：通过高T_high、低阈值T_low对上述非极大值抑制结果进行边缘提取，具体判断规则如下：

式中G(x,y)表示航空高光谱油膜指数梯度影像，x、y分别代表横纵坐标，T_high为高阈值，T_low为低阈值，I₁、I₂为边缘标志图像，初始置0，置1表示为初始边缘，由于I₁图像采用高阈值提取，故其标识的边缘确定为所需边缘，但最终提取结果还将包括I₂图像中标识的部分或全部边缘，本步骤通过梯度值图像的直方图统计，自适应获取最大梯度值|G|_max，根据这一参数构建高、低阈值，其计算公式如下：

然后进行边缘连接，I₁标识的边缘存在间断，需要进行连接修补，具体步骤为：a．扫描I₁图像，当I₁(x₀,y₀)＝1时，在I₁中跟踪起始点为S(x₀,y₀)的边缘，直到该边缘的终点E(x₁,y₁)；b．查看I₂中与E点对应的E′点周围8邻域内是否有点R'(x',y')，使得I₂(x',y')＝1，若有，则将I₁中与R′对应的R点置为1，即I₁(x',y′)＝1；再重复a步骤，继续查找最终边缘点，直到无法从I₂中提取最终边缘点，从而完成一条边缘检测，并标记为已检测；c．当完成一条边缘检测后继续在I₁中查找新的边缘起始点，并重复a、b步骤。最终完成所有边缘检测，形成航空高光谱遥感油膜边缘影像。

7.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤6）为：对航空高光谱油膜指数影像的每个像元识别薄油膜并采用窗口统计方式进行滤波，形成薄油膜分布遥感影像：通过阈值法识别薄油膜分布，具体判断规则如下： 

上式中Thred为判断阈值，I_oil(x,y)为航空高光谱油膜指数影像中点(x,y)处的值；本方法采用0.75这一阈值提取薄油膜分布，形成二值影像；然后采用窗口统计方式对薄油膜分布进行滤波，其具体方法为：在薄油膜分布图中，遍历每个像元，每次遍历采用如下判据判断像元是否为真正的薄油膜像元：当像元不为薄油膜时，则维持原来的判断；当像元为薄油膜时，以该像元为中心统计9×9邻域窗口内的薄油膜像元数，若总数超过9个，则将其判断为薄油膜，反之，则不为薄油膜；这样就形成航空高光谱薄油膜分布遥感影像。

8.根据权利要求1所述的一种基于油、水光谱特征差异的航空高光谱遥感自动识别海表溢油的方法，其特征在于所述的步骤7）为：在航空高光谱遥感油膜边缘影像、航空高光谱薄油膜分布遥感影像的基础上采用连通分析法来提取海表溢油全分布遥感影像：首先，将油膜边缘叠加到薄油膜分布图上，然后通过四连通方式对非薄油膜像元进行连通标记，由于清洁水体像元基本是大面积连续分布，厚油膜分布在薄油膜内部，其斑块面积远小于清洁水体，因此，可将连通区域大小在阈值300以下的区域判为厚油膜，最后综合厚、薄油膜分布得到海表溢油全分布遥感影像。 

Images