湿地系统是全球自然生态环境的重要组成部分，与海洋、森林等并称为全球三大生态系统，其在调节气候、补充地下水、蓄洪防旱、降解污染、土壤侵蚀控制等方面发挥着重要作用^[1]。黄河三角洲地处陆海相互作用的敏感地带，其河口堆积来自陆域径流的大量泥沙，形成了许多新生湿地，同时又受到海水波流动力侵蚀影响，湿地资源变化较为频繁^[2]。

目前，黄河三角洲湿地面临断流、污染、海水入侵等威胁，生态破损严重，湿地面积严重萎缩，生态功能逐渐退化^[3]。因此，对黄河口三角洲湿地进行分类、识别，了解黄河三角洲湿地的动态变化及其影响因素，对该地区湿地资源的合理利用、开发保护具有极为重要的意义。

合成孔径雷达（synthetic aperture radar, SAR）数据具有全天候、全天时、高空间分辨率等特点，重访观测周期通常可达一两周至一两个月左右，空间上可以连续获取，其侧视成像、多极化能力可以反映地物纹理特征、地表植被散射机制，上述特点有利于开展湿地分类、识别与监测。然而，SAR作为一种相干测量系统，其信号受到相干斑噪声的干扰，进而会影响其地物识别精度。

由于光学数据能反映地物的光谱特征信息，在一定程度上可以弥补极化SAR数据在地物识别方面的缺陷。与单一数据源相比，将光学与雷达数据结合用于地物分类，能减少或消除目标地物可能存在的多义性、不完全性、不确定性和差异性^[4]。

已有研究表明，结合光学与微波雷达影像各自的优点，可增加数据所含的有用信息，以增强对地物的识别能力。例如，文献[5]融合RADARSAT-1 SAR和Landsat专题制图仪（thematic mapper, TM）数据，与原始TM影像分类结果相比，分类精度有较大的提高；文献[6]结合环境卫星多光谱影像与Envisat ASAR单极化影像，协同分类精度比单独利用多光谱数据提高约5%；文献[7]利用高分三号（GF-3）极化SAR和Sentinel-2A光学数据对洞庭湖盆地一年一熟、一年两熟水稻等作物进行分类，其总体精度达到了85.3%，且高于两者单独分类的总体精度。

本文利用黄河三角洲地区Sentinel-2B光学和GF-3雷达影像，结合该湿地地区7类典型地物，在辐射校正、配准等预处理基础上，通过主成分变换、波段指数计算、极化分解等方式提取不同类型的特征协同集成，开展最大似然法（maximum likelihood, ML）、基于分类回归树（classification and regression tree, CART）算法的决策树（decision tree, DT）、支持向量机（support vector machine, SVM）方法的有监督分类，并对分类结果进行分析与评价。

1 研究区域与数据 1.1 研究区域

黄河三角洲位于山东省东营市黄河入海口，东邻莱州湾，北靠渤海湾。黄河三角洲湿地主要指1976年5月黄河尾闾改道至清水沟流路入海至今形成的新生湿地，介于118°42′E~119°20′E、37°34′N~37°59′N之间，面积达1 144.8 km²。黄河三角洲位于北半球暖温带欧亚大陆东部，受海陆热力性质差异影响形成温带季风气候，四季分明，夏季高温多雨，冬季寒冷干燥^[8]。

本文研究区域范围如图 1中红色框所示，图 1中的底图为Landsat-8陆地成像仪（operational land imager, OLI）5/6/4波段合成影像。

1.2 研究数据

目前，黄河三角洲地区的湿地分类研究主要基于光学遥感卫星和国外SAR卫星（如RADARSAT-2、ALOS-2等）^[9-10]。

Sentinel-2B卫星由欧洲空间局于2017年3月发射，与2015年6月发射的Sentinel-2A组成星座，具有多光谱成像（13个波段）、宽幅（290 km）、高分辨率（10~60 m）等特性，侧重于监测全球陆地和沿海地区的地物覆盖变化。

国产GF-3 SAR卫星具备1~500 m空间分辨率、10~650 km成像幅宽、12种成像模式的观测能力，是我国实施海洋开发、陆地环境资源监测和防灾减灾的重要支撑^[11-13]。

如表 1、表 2所示，本文选取2017-10获取的Sentinel-2B多光谱影像与GF-3全极化条带模式SAR影像作为实验数据。

2 数据预处理 2.1 Sentinel-2B预处理

Sentinel-2B L1C级数据是经过几何精校正的正射影像，但是没有进行辐射定标和大气校正，可以通过欧洲空间局专门用于实现Sentinel-2 L1C级到L2A级产品数据转换的Sen2Cor插件获得经过辐射定标和大气校正的大气底层反射率数据。考虑到GF-3分辨率为8 m，本文仅提取Sentinel-2B分辨率为10 m的4个波段（波段2/3/4/8，分别对应蓝光、绿光、红光、近红外波段）用于分类特征的提取。

本文采用的Sentinel-2B预处理及地物分类详细流程如图 2所示，图 2中的PCA1、PCA2为主成分变换（principal component analysis, PCA）后的首要两个主成分分量，I_NDVI为归一化植被指数，I_BI为土壤亮度指数，I_NDWI2为第二类归一化水体指数。

2.2 GF-3预处理

为了获得较好的分类特征，同样需要对GF-3全极化（HH、HV、VH、VV）单视复数影像（single look complex, SLC）数据进行预处理。本文首先进行辐射定标获得后向散射系数，然后使用改进型Lee滤波进行极化滤波处理，以减弱相干斑等噪声的影响，以便获得准确的强度、纹理、极化散射等特征信息。本文采用的GF-3预处理及地物分类详细流程如图 3所示，为了将Sentinel-2B预处理结果与GF-3预处理结果进行对比，将两者预处理结果放在了一起，结果如图 4所示。

3 特征提取与分类 3.1 分类算法

本文采用目前常用的3种经典监督分类算法进行分类^[14]。

1）基于统计模型的最大似然法。该方法建立在贝叶斯判决规则上，通过似然函数判断类别。由于遥感信息统计分布具有高度的复杂性和随机性，使用该方法常会产生较大误差。

2）基于知识的决策树分类。该方法通过对训练样本的特征属性和目标变量建立适当的树状分类结构进行分类，具有灵活、直观、效率高等特点。常用的决策树算法有ID3、C4.5、CART等，本文使用CART决策树分类算法。

3）基于统计学习理论的支持向量机法^[15]。该方法利用非线性变换将输入空间的线性不可分问题转化到高维空间，根据间隔最大化准则构造最优分类超平面^[16]。由于在处理非线性数据和高维特征空间识别方面具有独特优势，该方法已广泛用于遥感影像分类和多源遥感影像信息融合等领域。

3.2 Sentinel-2B特征提取与分类

如图 2所示，Sentinel-2B数据经过预处理后的4个波段之间信息并不独立。为了消除波段间的相关性，减少信息冗余，通过PCA提取两个独立的主成分分量PCA1、PCA2（贡献率大于99%）。结合研究区植被、裸地、水体多的特点，利用提取得到的4个波段计算3种指数特征，分别是I_NDVI、I_BI和I_NDWI2。其中，I_NDVI可以反映植被的密度和健康程度；I_BI对土壤的亮度敏感，与土壤的湿度及盐度相关，可在一定程度上反映裸土的分布；I_NDWI2可以检测陆地表面水的范围，即使在湿地环境下，仍有一定的识别水体能力。

3种指数特征的计算方式如下：

${{I}_{\text{NDVI}}}=\frac{{{\rho }_{\text{NIR}}}-{{\rho }_{\text{RED}}}}{{{\rho }_{\text{NIR}}}+{{\rho }_{\text{RED}}}}=\frac{{{B}_{8}}-{{B}_{4}}}{{{B}_{8}}+{{B}_{4}}}$

(1)

${{I}_{\text{BI}}}=\sqrt{\frac{{{\rho }_{\text{GREEN}}}+{{\rho }_{\text{RED}}}}{2}}=\sqrt{\frac{{{B}_{3}}+{{B}_{4}}}{2}}$

(2)

$ {{I}_{\text{NDWI}2}}=\frac{{{\rho }_{\text{GREEN}}}-{{\rho }_{\text{NIR}}}}{{{\rho }_{\text{GREEN}}}+{{\rho }_{\text{NIR}}}}=\frac{{{B}_{3}}-{{B}_{8}}}{{{B}_{3}}+{{B}_{8}}} $

(3)

式中，ρ_NIR、ρ_RED、ρ_GREEN表示近红外、红、绿波段的反射率；B₃、B₄、B₈表示3、4、8波段的波长

结合文献[1-3, 17]及实地调查资料，黄河三角洲湿地的典型地物主要包括：

1）河流；

2）水塘，包括人工养殖塘、自然池塘、水库等；

3）居民地；

4）草地；

5）裸地，包括滩涂、盐碱地等空地；

6）林地；

7）庄稼地。

对应的训练样本分布如图 5所示。

根据图 2将Sentienl-2B得到的5个分类特征（PCA1、PCA2、NDVI、BI、NDWI2）和样本数据输入到ML、DT、SVM等3种监督分类器，分类结果如图 6（a）、6（b）、6（c）所示。

3.3 GF-3特征提取与分类

GF-3极化散射特征主要利用Freeman三分量分解法和H-A-α分解得到。

Freeman三分量分解法是以物理实际为基础，对极化协方差矩阵C₃建立3种基本散射机制模型，分别为适度粗糙表面散射、随机方向偶极子组成的体散射、由正交平面构成的二次散射^[18-19]。假设表面散射P_S、二次散射P_D和体散射分量P_V互不相关，由每种散射机制成分的权重值（f_S、f_D、f_V）可以得到目标极化总功率P_SPAN，即上述每种独立散射机制成分统计量的叠加之和^[20]：

$ \begin{array}{*{35}{l}} {{P}_{\text{SPAN}}}={{\left| {{S}_{\text{HH}}} \right|}^{2}}+2{{\left| {{S}_{\text{HV}}} \right|}^{2}}+{{\left| {{S}_{\text{VV}}} \right|}^{2}}= \\ \ \ \ \ \ \ \ {{P}_{S}}+{{P}_{D}}+{{P}_{V}} \\ \end{array} $

(4)

其中，

总功率P_SPAN是由HH、HV、VH、VV 4通道强度叠加得到，强度特征较为强烈，叠加后部分正负噪声干扰得到了消除或减弱。考虑到不同地物纹理特征存在差异，结合纹理特征可以提升地物分类精度，因此本文利用灰度共生纹理矩阵作用于极化总功率交换器端口分析器（switched port analyzer，SPAN）图像来提取纹理特征，共得到8个纹理特征^[21]。由于这些纹理特征并非相互独立，为了减少数据冗余量，通过7类地物感兴趣区域（region of interest, ROI）样本均值来提取纹理特征。另外，纹理特征单位或数量级存在差异，需要将其归一化才能区分差异。

图 7给出了GF-3 ROI样本纹理特征均值曲线。由ROI样本区域均值的统计分析可知，均值、均匀性、非相似性3个特征对7类样本的区分较为明显，故选择这3个纹理特征参与分类。

H-A-α分解利用散射矩阵S转换得到相干矩阵〈[T₃]〉，其中〈[T₃] 〉是半正定埃尔米特矩阵。根据矩阵分解理论，可将其扩展为3个独立目标之和，即用3个散射矩阵来表示3个独立目标^[22]，即：

$ \left\langle \left[ {{T}_{3}} \right] \right\rangle =\sum\limits_{j=1}^{3}{{{\lambda }_{j}}}{{\mathrm{u}}_{j}}\mathrm{u}_{j}^{\text{*T}}=\sum\limits_{j=1}^{3}{{{\lambda }_{j}}}{{T}_{j}} $

(5)

式中，〈[·]〉表示多视处理；上标*表示共轭；上标T表示转置；T_j、λ_i和u_j分别表示第j个独立目标、特征值、特征向量。

H-A-α分解的3个二阶参数是〈[T₃]〉的特征值和特征矢量函数，其定义如下。

1）极化散射熵H:

$ H=-\sum\limits_{k=1}^{3} P_{k} \log _{3}\left(P_{k}\right) $

(6)

其中，${{P}_{i}}={{\lambda }_{i}}/\sum\limits_{k=1}^{3}{{{\lambda }_{k}}},\sum\limits_{k=1}^{3}{{{P}_{i}}}=1$。

2）平均α角：

$ \alpha ={{P}_{1}}{{\alpha }_{1}}+{{P}_{2}}{{\alpha }_{2}}+{{P}_{3}}{{\alpha }_{3}}$

(7)

式中，α₁、α₂、α₃的大小指示主要的散射机制

3）极化各向异性度A：

$ A=\frac{{{\lambda }_{2}}-{{\lambda }_{3}}}{{{\lambda }_{2}}+{{\lambda }_{3}}}$

(8)

根据图 3所示流程图，本文将GF-3 SAR影像得到的强度特征SPAN、3个Freeman三分量分解极化散射特征、3个H-A-α分解极化散射特征、3个纹理特征作为敏感特征与训练样本输入ML、DT、SVM等3种监督分类器，分类结果见图 6。

3.4 Sentinel-2B和GF-3协同分类

Sentinel-2B和GF-3的精确配准是实现两者协同分类的关键，本文通过野外实地调查及目视选择控制点实现Sentinel-2B和GF-3数据准确配准。将配准后的Sentinel-2B和GF-3影像特征（共15个）进行叠加，与训练样本一起输入到ML、DT、SVM等3种监督分类器，实现两者的协同分类，分类结果如图 6（g）、6（h）、6（i）所示。

3.5 结果分析与讨论

本文采用生产者精度、用户精度、总体精度、Kappa系数作为定量评估不同数据与分类算法的性能指标^{[21, 23]}，结合图 6和表 3可知：

1）由3种数据的分类结果来看，极大似然法的精度最低，决策树和支持向量机的精度相对较高，且两者差别微弱，但仅凭精度难以区分后两个算法的好坏，因为取不同参数时，不同分类算法的分类精度也会有变化。

2）从生产者精度和用户精度综合来看，居民地和庄稼地的分类精度较低，这主要是因为它们在图中范围分布数量较小（影像像素少），且往往是局部分布，容易受到其邻近地物拓扑关联（邻接、包含等）分布的影响。同样地，尽管两者分类精度较低，但它们占比少，对于总体精度、Kappa系数影响小，只要其他地类的分类精度足够高，仍然可以得到较高的总体精度和Kappa系数。

3）其他5类地物的分类精度相对较高，特别是裸地在所有分类方案下的精度普遍达到90%以上，这可能与该地物在该地区面积较大且连片分布有关。水塘和河流、草地和林地之间存在部分混淆，主要由于其特征较为相似且存在邻接边界。

4）如图 6（d）、6（e）、6（f）所示，GF-3大部分海域被分为裸土，这是由于裸地训练样本包含了南部潮间带地区，该地区裸地定期受涨潮影响而被海水覆盖，而训练样本没有考虑海域部分，并且两者均有较强的表面散射特征，因此海域可能被错分为裸土区域。

5）总体上，GF-3的总体分类精度和Kappa系数略低于Sentinel-2B分类，而其居民地分类精度要优于Sentinel-2B，这可能是由于居民地在GF-3影像二次散射较强烈，有利于将其区分出来。

6）两者协同分类后各类地物的分类精度比单独分类得到一定程度的提升，尤其是Sentinel-2B和GF-3分类精度都较低的居民地和庄稼地。协同分类的总体精度和Kappa系数均达到较高数值，主要是由于光学数据与SAR数据的协同分类得到的是“1+1 > 2”的信息融合与互补，更加有利于地物识别。

4 结语

本文利用Sentinel-2B多光谱数据和GF-3全极化SAR数据分别提取了黄河三角洲地区7类湿地地物的光谱、指数、极化散射、纹理等特征，分别利用最大似然法、决策树法、支持向量机3种算法实现了有监督协同分类，总体精度可达90%以上，均高于两者单独分类的总体精度，表明了基于两种不同种类的卫星遥感数据源进行湿地协同分类的可靠性和优势。

本文后续研究将考虑加入不同波段的多光谱数据和SAR数据（如L/X波段等）开展更细致的地物分类，以便准确反映该地区的湿地资源状况。结合黄河三角洲地区长时间序列多源遥感数据，有望实现该地区地物覆盖的动态分布监测。

致谢: 感谢国家卫星海洋应用中心提供高分三号SAR数据以及欧洲空间局提供Sentinel-2B多光谱数据。