我国海域辽阔，水产养殖业发展空间巨大^[1]。近年来，海水养殖的规模和密度不断加大^[2]，由此带来的近海养殖污染等问题也日益严峻。因此，科学准确地监管近海水产养殖业具有重要的意义和价值。传统的养殖区监测通常采用海上实地考察的方式，易受气象条件等因素限制，并且耗时、耗力。近年来，卫星遥感技术不断发展，其数据具有宏观性、实时性和长时序性等特点，在海洋资源和环境调查等领域已展现出其特有的技术优势。因此，遥感技术可以高效率、低成本的实现养殖区信息的精确提取，是辅助水产养殖业健康发展的重要技术手段。

当前利用遥感数据的近海养殖区提取研究大多基于高分辨率影像，虽然数据空间分辨率较高，但光谱通道有限，幅宽较小，并且影像获取成本较高，适合局部、小区域尺度养殖区域的高精度提取，不利于大面积的动态监测^[3-4]。相比之下，长时序、中分辨率卫星数据的获取更加方便快捷，并且提供了更多通道的光谱信息，有助于提高养殖区信息的解译能力。在利用遥感影像进行近海养殖区识别提取的研究中，Hazarika等^[5]利用Landsat数据，通过目视解译对泰国东海岸的虾池养殖区进行了监测。马艳娟等^[6]针对ASTER影像光谱信息，提出了构建水体指数和波段运算函数的提取方法。然而，受海流、悬浮泥沙等因素的影响，近海养殖区常存在同物异谱或同谱异物等现象，若仅利用影像光谱信息，无法克服“椒盐”噪声的影响，导致养殖区像元的误判。武易天等^[7]基于Landsat数据，在增强养殖区光谱信息的基础上，采用局部标准差自适应分割法结合纹理信息，分别提取了山东桑沟湾和福建三都澳养殖区的数据。陈思明^[8]用Landsat数据，结合特征指数与面向对象检测法完成了三沙湾海水养殖区数据的提取。结合纹理特征或者基于面向对象等框架的空间信息后处理方法，虽然在一定程度上可以降低光谱差异的影响，但是受空间分布和尺度等不确定性因素的限制，提升效果往往有限。

综上，本文针对哨兵二号中分辨率影像，以辽东半岛东部的海王九岛养殖区作为典型研究区，提出一种基于谱空结合的中分辨率遥感影像近海养殖区提取方法，通过综合利用影像光谱信息和空间信息实现近海养殖区的精准提取。

1   材料与方法

1.1   研究区概况

海王九岛位于辽东半岛东部海域（123°02′13″E－123°06′23″E，39°25′32″N－39°30′49″N），海域面积约为600 km²，由大王家岛、小王家岛、元宝岛、寿龙岛等组成，具体位置见图1。海王九岛自然条件优越，是国家级海岛森林公园，具有良好的旅游发展前景。

图  1  辽东半岛东部海域研究实验区图

Fig.  1  The experimental area of east Liaodong Peninsula

下载: 全尺寸图片 幻灯片

海王九岛海域饵料资源充足，各种贝类和经济鱼类资源十分丰富。但近年来，随着水产养殖业的迅速发展，无序化开发活动日益加剧，海域生态系统健康面临着严重的威胁。对该海域进行养殖区的精确提取可以为相关部门对养殖用地的科学规划和管理提供参考依据，为海王九岛长期健康发展提供辅助信息。

通过实地勘察，笔者了解到研究区以浮筏、网箱、参圈3种养殖模式为主，养殖种类有刺参、扇贝和牡蛎等。3种养殖模式的遥感影像及现场图详见表1。

表  1  研究区养殖类型说明

Tab.  1  Description of types of mariculture in the study area

养殖类型	现场图	影像图	说明
浮筏养殖			浮筏由浮力球、绳索搭建而成， 在遥感影像中呈规则的黑色条带状分布； 养殖内容以海湾扇贝、牡蛎等贝类为主
网箱养殖			网箱由塑料、木材等材料 组成，在影像中呈规则的亮白色 分布；养殖内容为各种经济鱼类
参圈养殖			参圈为工人在海边用石头围砌 而成的养殖区或人工挖掘出的池塘，在 影像中呈暗色；养殖内容一般为虾或刺参

下载: 导出CSV 

| 显示表格

1.2   数据源

哨兵二号载有覆盖13个谱段的多光谱成像仪，包含A、B两颗卫星，双星每5天更新一次成像数据，各波段的空间分辨率分别为10 m（Band₂、Band ₃、Band ₄、Band ₈）、20 m（Band₅、Band₆、Band₇、Band_8A、Band₁₁、Band₁₂）、60 m（Band₁、Band₉、Band₁₀）。实验影像为2020年3月14日获取的哨兵二号L2A级数据，通过欧空局哥白尼数据中心下载（https://scihub.copernicus.eu/dhus/#/home）。

利用高分二号影像进行目视解译，结合现场勘查，以此制作模型训练和精度评价的验证数据。高分二号影像源于辽宁省国土资源厅2020年2月18日4 m分辨率的研究区影像。实地勘查于2021年4月在海王九岛进行，获取浮筏、网箱、参圈验证像元数量分别为132个、26个、14个。结合高分二号影像及现场勘查的标注结果制作验证数据，在研究区标注浮筏、网箱、参圈3类样本分别为65497个、2073个、2860个。

1.3   研究方法

1.3.1   技术路线

首先，利用SNAP和ENVI软件完成研究区域的重采样、影像裁剪等预处理操作；其次，为进一步增强养殖区目标的特征表达能力，在原始数据基础上进行指数计算，补充原始数据的光谱和纹理特征；再次，对研究区浮筏、网箱、参圈、海水4类地物进行分类，从各类地物中选取一定的训练样本输入支持向量机模型（support vector machine，SVM）进行训练，得到光谱信息的初分类结果；最后，结合马尔科夫随机场模型（markov random field，MRF），通过空间信息处理对初分类结果进行修正和调整。此外，为了评估本文方法的准确性和优势，选取同为监督分类方法的最大似然分类法（maximum likelihood estimation，MLE）以及未添加空间后处理的SVM算法作为对比实验，实验中各算法均采用多组交叉验证的方式进行运算测试。此外，结合现场实际勘查得到的照片和高分辨率影像，根据数据判读、专家解译和实地调研对提取结果进行精度评价，具体技术路线如图2所示。

图  2  技术路线

Fig.  2  Technology roadmap

下载: 全尺寸图片 幻灯片

1.3.2   数据预处理

本文采用的哨兵二号L2A级数据已经过大气校正和辐射校正。首先，将哨兵二号20 m、60 m分辨率的单波段数据在SNAP中重采样为10 m分辨率，再利用ENVI软件进行波段合成、影像裁剪等操作。此外，为排除非感兴趣区海岛的影响，本文利用改进的归一化差异水体指数（modification of normalized difference water index，MNDWI）^[9]对研究区影像进行水陆分离和海岛掩膜。

1.3.3   指数计算

从图1和表1可以看出，浮筏和网箱养殖区均有较明显的空间纹理特征，因此，本文利用灰度共生矩阵（gray level co-occurrence martrix，GLCM）^[10]对养殖区空间纹理信息进一步增强。选择对养殖区信息敏感的红、绿、蓝3个波段，进行灰度共生矩阵的计算。经多次试验发现，均值、对比度、同质性、非相似性这4种纹理特征参数组合的方式对养殖区纹理增强效果最好，因此选取该组合方式。针对纹理特征计算得到的12维数据，进一步利用主成分分析进行特征提取，以提高特征表达能力，减少数据冗余。本研究选取主成分分析得到的前3个主成分（包含97.12%的数据信息量）作为增加的空间纹理特征。

叶绿素浓度和悬浮泥沙浓度对近岸水域的反射特性影响较大，导致了近岸养殖区复杂的水色背景。为增强水产养殖区，本文在原始波段的基础上添加了叶绿素浓度指数（chlorophyll a concentration，CHL）、悬浮泥沙指数（total suspended matter concentration，TSM）^[11]。同时添加归一化的水域自动提取指数（automated water extraction index，AWEI）^[7]对非水产养殖区（海水）进行增强，使水产养殖区与背景海水光谱差异更明显。水体指数具体表达式详见表2，指数计算结果详见图3。为便于后文对比实验的说明，将其统一命名为WIs。将以上6个新特征与原始影像叠加，形成新的实验数据应用于后续养殖区提取实验中。

表  2  指数计算中水体指数表达式

Tab.  2  The expression of water index in index calculation

水体指数（WIs）	表达式
TSM	${\rm{TSM = }}0.028{b_2}{\rm{ + }}0.019{b_3} - 5.31{b_3}/{b_2} + 0.537 $ （1）
AWEI	$ {\rm{AWEI = (4}}({b_{\rm{3}}}{\rm{ - }}{b_{{\rm{11}}}}{\rm{) - (0}}{\rm{.25}}{b_{\rm{8}}}{\rm{ + 2}}{\rm{.75}}{b_{{\rm{12}}}}{\rm{)) / (}}{b_{\rm{3}}}{\rm{ + }}{b_{\rm{8}}}{\rm{ + }}{b_{{\rm{11}}}}{\rm{ + }}{b_{{\rm{12}}}}{\rm{)}}$ （2）
CHL	${\rm{CHL}} = 8.48{b_4}/{b_1} $ （3）
注：b1~12分别表示哨兵二号影像对应波段

下载: 导出CSV 

| 显示表格

图  3  研究区指数计算结果图

Fig.  3  The result of index calculation in the study area

下载: 全尺寸图片 幻灯片

1.3.4   养殖区谱空综合提取方法

针对影像中养殖区的多维小样本的数据特点，本文选择SVM分类器对养殖区进行初步提取。但研究区水色背景较复杂，单一谱分类器又易受噪声以及同物异谱等因素的影响，导致对水产养殖区提取效果不佳。因此，本文在SVM分类结果的基础上进一步添加MRF空间后处理，对初分类结果进行空间上的修正。通过综合利用光谱和空间信息，实现水产养殖区的精准分类提取。

（1）支持向量机

支持向量机是在统计学理论的基础上发展起来的一种经典分类器，为机器学习领域的标准工具。基本原理是寻找一个超平面使待分类的样本间隔最大化，它将分类问题转化为凸二次规划问题的求解，学习效率高、抗噪声性能强，适用于小样本监督分类^[12-13]。在近海养殖区提取研究中，SVM可依据少量的先验训练样本，提供快速有效的养殖区初分类结果。

如图4所示，以二分类问题为例，假定训练数据集为$\left[\left(x_{i}, y_{i}\right), i=1,2, \cdots, N\right]$，其中$x$是训练样本，$ y $是对应的类别标记（$ y \in\{+1,-1\} $）。分类面为将一个空间按照类别切分为两部分的平面，在线性分类中，分类面表达式为$ f(x)=w^{T} x+b $，其中$ \omega $为系数，$ b $定义了偏移量。SVM要求找到能将两类样本分开，且满足分类间隔最大的最优分类面，两类样本分类间隔大小可表示为：$ Margin = 2/||w||$。

图  4  SVM原理示意图

Fig.  4  SVM principle diagram

下载: 全尺寸图片 幻灯片

上述最优分类面的求解问题可以表示为一个优化问题，该问题通过拉格朗日（Lagrange）函数求解，基于Karush Kuhn Tucker互补条件求解得到的系数，代入线性分类面，得到线性SVM的目标函数形式为：

实际应用中的待分类图像通常为线性不可分的多维数据（如本文的研究区影像），此时类别之间的分类用非线性曲面才能更好判别。通常需要构造核函数将样本点空间映射到外一个特征空间，以达到更好的分类效果。本文采用的是高斯径向基（radial basis function，RBF）核函数，函数表达式为：

式中：$ \sigma $为控制函数径向参数范围的参数。

基于高斯径向核函数的SVM（Kernel SVM）的目标函数可以表示为：

式中：$ \delta(x) $为映射后的样本集。

（2）马尔科夫随机场

研究区影像中，各类养殖类型形状较为规则，相邻像元间存在着空间上的相互联系，邻近像元属于同一类别的概率较大。以哨兵二号影像中网箱养殖区为例（图5），MRF基本假设为与网箱像元在空间位置上相邻的像元大概率也属于网箱，通过利用这种邻域空间相关性信息对谱分类结果进行修正，提升养殖区识别分类效果^[14-15]。

图  5  邻域空间相关性示意图

Fig.  5  Schematic diagram of neighborhood spatial correlation

下载: 全尺寸图片 幻灯片

本文依据MRF与Gibbs分布的等价关系^[16]，并假设各个类别的分布符合多维高斯分布，通过最大似然概率进行估计，得出图像上任一像元基于马尔科夫随机场分类判别函数表达式为：

其中，$\delta ({\omega _{kx}},{\omega _{\partial x}})$为Kroneker函数，表达式为：

式中：$ m_{k} $和$\sum\nolimits_k {}$分别表示类别$ k $的均值矢量和协方差矩阵；$ \beta $为权重系数；${\omega _{kx}}$表示将像元$x$分到类别$ k $中；${\omega _{\partial x}}$表示像元x的邻域像元的类别。该式表达了加入像元空间上下文特征后的分类概率。

（3）支持向量机结合马尔科夫随机场模型

由于当前研究区的部分养殖区像元与海水像元光谱特征相近，存在同物异谱、同谱异物的现象，单一谱分类器往往很难精准识别，这时需要考虑综合像元的光谱与空间信息。MRF是一种有效利用光谱特征和空间特征的分类器。如公式（9）所示，整合马尔科夫的分类过程可以表示为由光谱项和空间项组成的概率表达形式，既能充分利用实验数据的光谱信息，又能结合空间信息对误判结果进行修正，从而提高养殖区识别效果。

式中：${a_i}(k)$定义为光谱项，是像元$ i $根据光谱特征判定属于类别$ k $的概率；${b_i}(k)$为空间项，是像元$ i $根据空间特征判定属于类别$ k $的概率；$ \beta $为权重系数。

本文利用SVM的初分类结果作为光谱项，再使用MRF进行空间后处理。结合SVM目标函数，利用Platt’s提出的后验概率改进形式，计算得到基于SVM的像元分类的后验概率分布为：

式中：A和B为参数，可以通过最小化交叉熵误差函数获取。

将公式（10）带入公式（9）中，可得SVM-MRF的目标函数为：

在对上述目标函数的优化过程中，将求解目标函数最大值转化为求解函数能量最小值问题，采用模拟退火算法计算。该算法根据Metropolis接受准则，使用冷却进度表参数控制算法进程，从而获得全局近似最优解^[17]。

1.3.5   精度评价方法

本文采用整体分类精度（overall accuracy，OA）、Kappa系数、用户精度（user’s accuracy，UA）对养殖区提取结果进行精度评价。

整体分类精度OA是指分类正确的样本个数占所有样本个数的比例，具体见式（12）。

Kappa分析是评价分类精度的多元统计方法，Kappa系数代表被评价分类比完全随机分类产生错误减少的比例，表达式见式（13）。

式中：c为类别数；${M_{j,i}}$表示类别j被识别为类别i的样本数量；n为测试样本的数量。

OA和Kappa系数越高，表示研究区的整体分类效果越好。为更好的评价各类养殖区的提取准确性，本文选择用户精度作为评价指标计算。

UA是表示分类结果中某一随机样本与地面实际类型相同的条件概率，表达式见（14）。

式中：${x_{i,i}}$为误差矩阵中第i行、第i列的观测点个数；${x_{i + }}$表示第i行的和。

2   结果与讨论

2.1   研究区提取结果

本文结合高分二号影像以及现场勘查的标注结果进行验证数据的制作。在研究区共标注浮筏、网箱、参圈三类样本分别为65497、2073、2860个，得到的地面真实参考图像如图6（a）所示。研究共设置五种算法的对比：MLE算法、SVM算法、添加水体指数的SVM算法（SVM_WIs）、添加水体指数和纹理信息的SVM（SVM_WIs+GLCM）、本文提出的谱空综合提取算法（SVM_WIs+GLCM-MRF）。为达到更好地对比效果，四组对比实验与本实验使用相同的训练样本数据集。对比试验及本文方法的提取结果图详见图6（b－f）。

图  6  水产养殖区提取识别结果

Fig.  6  Recognition results of mariculture area

下载: 全尺寸图片 幻灯片

对比地面参考图像可以发现，MLE提取效果较差，大量海水像元被误判为浮筏，岛屿周围的海水像元被误判为网箱和参圈，部分浮筏养殖区连成一片，无法辨别出筏区条带状的空间养殖特征；SVM整体识别效果较好，误判的海水像元虽有明显减少，但仍存在较多海水被错分为浮筏的现象（图中的黄色散点），同时养殖区内部少量像元被误判为海水；添加水体指数后，养殖区的识别效果虽稍有提升，但错分现象并没有得到很大改善；进一步添加纹理信息后，水产养殖区信息得到增强，识别出的水产养殖区纹理更为清晰，浮筏养殖区提取效果有明显提升，然而错分和误判现象仍然较明显；与此相比，本文添加MRF空间后处理后，前四种方法中出现的大量误判的海水像元几乎全部得到修正，各类地物的识别效果均有较大提升，尤其是浮筏养殖区，纹理更为清晰。整体而言，本文方法提取的养殖区更接近地面参考图像。

为进一步对比本文方法对养殖区提取效果的细节，取图6中白色框线内的区域进行局部放大，如图7所示。图7（a）和图7（b）为现场勘察影像，图7（c）为局部放大区域的遥感影像图。结合现场照片和遥感影像可以看出，由于受海流、海风、悬浮泥沙等自然因素以及养殖设施自身下沉等影响，部分水产养殖区与背景海水的光谱差异较小，两者易于混淆。图7（d）为仅考虑添加纹理信息和水体指数时的提取结果，浮筏养殖区条带之间出现明显粘连的现象，矩形网箱养殖区的形状也不规整。而添加MRF进行空间后处理后，从图7（e）中可以看出误判像元得到有效的修正，浮筏养殖区提取结果得到明显改善，空间纹理特征更清晰，条带之间的间隔更明显，网箱呈现出规则的矩形特征，与原始影像基本吻合。

图  7  养殖区现场图及局部提取结果

Fig.  7  Scene map of the mariculture area and partial extraction results

下载: 全尺寸图片 幻灯片

2.2   精度评价

研究区水产养殖区提取的精度评价结果详见表3。由于MLE在提取时将大量像元误判为海水，导致三类水产养殖区的提取精度较低，而SVM对养殖区的提取精度明显高于MLE算法。在SVM算法基础上添加水体指数后（SVM_WIs），养殖区提取精度有一定提升。其中，浮筏的精度提升约1.04%，网箱和参圈的提升约为0.24%和0.19%；添加纹理信息后（SVM_WIs+GLCM），相较于SVM_WIs，养殖区提取精度进一步提升。其中，浮筏的精度提升为1.51%，网箱和参圈的提升分别为0.51%与1.18%；与此相比，在利用本文谱空综合提取方法后，养殖区提取精度有了显著提升。三类养殖区的平均精度（浮筏、网箱、参圈精度的平均值）为83.48%，三类精度分别提升了5.55%、1.97%和4.84%。本文方法取得了最高的OA和Kappa系数。OA高达94.46%，较前四种方法分别提升22.20%、1.79%、1.52%和1.13%。Kappa系数为0.7604，较前四种方法分别提升0.4291、0.0624、0.0522、0.0378。

表  3  研究区提取精度评价表

Tab.  3  Evaluation table of extraction precision in study area

方法	MLE	SVM	SVMWIs	SVMWIs+GLCM	SVMWIs+GLCM-MRF
浮筏	30.78%	65.37%	66.41%	67.92%	73.47%
网箱	34.48%	79.18%	79.42%	79.93%	81.90%
参圈	28.57%	88.86%	89.05%	90.23%	95.07%
海水	99.49%	97.59%	97.69%	97.82%	97.81%
OA	72.26%	92.67%	92.94%	93.33%	94.46%
Kappa	0.3313	0.6980	0.7082	0.7226	0.7604

下载: 导出CSV 

| 显示表格

结合各算法的养殖区提取结果与精度评价（图6与表3）可知，水体指数的添加，使得海水和浮筏的混分现象得到了些许改善，这表明两种类别间少部分的“同谱异物”现象可以通过增加光谱信息进行区分，但是仍然存在诸多混分像元，需要进一步改进处理；由于影像中的养殖区呈规则的空间纹理特征、浮筏呈条带状、网箱呈现矩形的形状特点，增加了灰度共生矩阵计算的纹理特征后，网箱与浮筏与海水可以进一步被区分，海水和养殖区的混分现象得到了进一步的改善；但由于水体自然因素以及养殖设施状态等影响，依然有一部分混分像元无法通过上述方式进行有效区分。在结合马尔科夫随机场模型后，像元空间上的相关性弥补了光谱特征在区分“同谱异物”像元时的局限性，修正了绝大部分被误判为养殖区的海水像元，正确识别出的浮筏养殖区的像元也明显增多，三类养殖区均达到了最高的提取精度。

综上所述，本文提出的基于支持向量机与马尔科夫随机场模型的养殖区谱空综合提取方法能精确提取近海养殖区。谱分类器结合空间后处理模型后，能有效纠正初分类结果中误判的像元，从而提升海王九岛研究区各养殖类型的提取精度。

3   结 论

本文利用哨兵二号卫星遥感影像，针对近海养殖区光谱与空间特征的利用不足问题，提出了一种谱空结合的养殖区自动提取算法。该方法通过添加水体指数和纹理特征增强养殖区信息，利用支持向量机分类器结合马尔科夫随机场模型有效挖掘了影像的谱空信息。该方法在研究区达到了94.46%的整体分类精度，相较于仅利用光谱信息和纹理特征增强的方法，各类养殖区的精度得到显著提升。研究区的海水养殖类型大部分为浮筏养殖，但近岸水体组分较为复杂，受悬浮泥沙等因素的影响，部分筏区光谱特征不强，较难识别，后续研究将针对该问题进行深入探讨，以提升浮筏养殖区的提取精度。