2022, Vol. 39 
2. 中国科学院大学电子电气与通信工程学院, 北京 100049
2. School of Electronic, Electrical and Communication Engineering, University of Chinese Academy of Sciences, Beijing 100049, China
合成孔径雷达(synthetic aperture radar, SAR)是一种高分辨率成像雷达，具备全天时、全天候、远距离观测目标的能力[1]。极化SAR(polarimetric SAR)是一种多通道多参数的SAR系统，与传统SAR系统相比可以获得更为丰富的目标散射特性[2]，在图像分类方面具备一定优势[3]。
地物分类是遥感图像解译的关键问题之一，在环境监测、灾害评估、城市规划等方面有着广泛应用。随着大数据处理技术的发展，学者们开始尝试将机器学习等算法应用在SAR数据处理中。文献[4]用深度卷积神经网络(convolutional neural network, CNN)对极化SAR进行分类，可以得到较高的分类精度，但其依赖大量的有标记训练样本，在数据量较少时分类效果不佳，而在此情况下机器学习等传统算法表现通常优于深度卷积神经网络[5]。
一般来说，特征参数种类越丰富，所描述的地物目标信息便越全面，然而由于特征之间冗余度和相关性的存在，不加区分的特征叠加往往会导致过高的计算复杂度，而分类效果却提升有限，同时SAR图像的相干斑噪声也一直影响着分类精度。针对这些问题，文献[6]提出一种基于特征筛选的二级分类结构对建筑物进行提取，作者利用随机森林(random forest, RF)[7]对特征进行筛选，并得到初级分类结果，将筛选后的特征输入到支持向量机(support vector machine, SVM)[8]得到二级分类结果，两种结果融合后获得最终结果。该方法在建筑物提取上效果较好，但由于RF没有区分决策树的分类能力，所有决策树在分类时默认权重相同，而实际上决策树分类能力有所差异，对权重适当调整可以提升分类器的总体精度。因此陈伟民等[9]针对该问题提出一种基于自适应提升(adaptive boosting, AdaBoost)改进RF算法，并在高光谱图像分类上验证了该算法的有效性。
为提升极化SAR图像地物分类精度，减少相干斑噪声对分类的影响，本文设计了一种基于AdaBoost改进RF(adaBoost random forest, ADA_RF)和SVM的二级分类算法。该方法采用AdaBoost提升RF性能的方法[9], 同时还引入了二级分类结构[6]。
该方法首先从极化SAR数据中提取特征数据集，利用分层采样法选取训练数据集和测试数据集，先得到ADA_RF分类器和特征重要性排名，再根据排名选择部分特征求出SVM的二级分类结果，通过邻域投票法融合两级结果，得出最终结果。
1 特征描述在极化SAR图像的分析过程中常用极化散射矩阵S来描述地物的极化散射特性，S矩阵经过变换之后又可以得到极化协方差矩阵C和极化相干矩阵T。
| $ \mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{S}}_{{\rm{HH}}}}}&{{\mathit{\boldsymbol{S}}_{{\rm{HV}}}}}\\ {{\mathit{\boldsymbol{S}}_{{\rm{VH}}}}}&{{\mathit{\boldsymbol{S}}_{{\rm{VV}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\left| {{\mathit{\boldsymbol{S}}_{{\rm{HH}}}}} \right|{{\rm{e}}^{{\rm{j}}\phi {\rm{HH}}}}}&{\left| {{\mathit{\boldsymbol{S}}_{{\rm{HV}}}}} \right|{{\rm{e}}^{{\rm{j}}\phi {\rm{HH}}}}}\\ {\left| {{\mathit{\boldsymbol{S}}_{{\rm{VH}}}}} \right|{{\rm{e}}^{{\rm{j}}\phi {\rm{VH}}}}}&{\left| {{\mathit{\boldsymbol{S}}_{{\rm{VV}}}}} \right|{{\rm{e}}^{{\rm{j}}\phi {\rm{VV}}}}} \end{array}} \right], $ | (1) |
| $ \mathit{\boldsymbol{T}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {\left\langle {{{\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} + {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right|}^2}} \right\rangle }&{\left\langle {\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} + {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right){{\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} - {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right)}^*}} \right\rangle }&{2\left\langle {\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} + {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right){\mathit{\boldsymbol{S}}}_{{\rm{HV}}}^*} \right\rangle }\\ {\left\langle {\left( {{\mathit{\boldsymbol{S}}_{{\rm{HH}}}} - {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right){{\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} + {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right)}^*}} \right\rangle }&{\left\langle {{{\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} - {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right|}^2}} \right\rangle }&{2\left\langle {\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} - {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right){\mathit{\boldsymbol{S}}}_{{\rm{HV}}}^*} \right\rangle }\\ {2\left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}{{\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} + {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right)}^*}} \right\rangle }&{2\left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}{{\left( {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}} - {{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right)}^*}} \right\rangle }&{4\left\langle {\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}} \right|} \right\rangle } \end{array}} \right], $ | (2) |
| $ {\mathit{\boldsymbol{C}}} = \left[ {\begin{array}{*{20}{c}} {\left\langle {{{\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}}} \right|}^2}} \right\rangle }&{\sqrt 2 \left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}^*} \right\rangle }&{\left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}}{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}^*} \right\rangle }\\ {\sqrt 2 \left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}^*} \right\rangle }&{2\left\langle {{{\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}} \right|}^2}} \right\rangle }&{\sqrt 2 \left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}}{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}^*} \right\rangle }\\ {\left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}{\mathit{\boldsymbol{S}}}_{{\rm{HH}}}^*} \right\rangle }&{\sqrt 2 \left\langle {{{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}{\mathit{\boldsymbol{S}}}_{{\rm{HV}}}^*} \right\rangle }&{\left\langle {{{\left| {{{\mathit{\boldsymbol{S}}}_{{\rm{VV}}}}} \right|}^2}} \right\rangle } \end{array}} \right]. $ | (3) |
式(2)和式(3)中：〈·〉表示集合平均，*表示共轭。在单站后向散射体制下，假定S满足互易性，即有SHV= SVH，因此通常定义极化SAR系统测量的散射总功率为
| $ \mathit{\boldsymbol{Span}} = {\mathop{\rm Tr}\nolimits} \left( {\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{S}}^{*{\rm{T}}}}} \right) = {\left| {{\mathit{\boldsymbol{S}}_{{\rm{HH}}}}} \right|^2} + 2{\left| {{\mathit{\boldsymbol{S}}_{{\rm{HV}}}}} \right|^2} + {\left| {{\mathit{\boldsymbol{S}}_{{\rm{VV}}}}} \right|^2}. $ |
式中: Tr(·)表示求矩阵的迹，T作为上角标时表示转置。
不同的极化目标分解方法描述地物散射特性的角度不同。极化目标分解方法主要分为3种，分别是基于散射模型的分解方法、基于特征矢量或特征值分析的分解方法和相干分解方法[2]。基于散射模型的目标极化分解方法有Freeman-Durden分解、Yamaguchi分解等。基于特征矢量目标极化分解方法有Cloude-Pottier分解、Holm分解和VanZyl分解等。相干分解是将S矩阵分解为目标散射机制的组合，有Pauli分解和Krogager分解等方法。考虑到地物目标一般较为复杂，单一的分解方法难以适用于不同场景，故本文从3种分解方法中各选取了几个典型特征组成极化特征集合。
纹理是应用于遥感图像分类的典型特征。它通过灰度空间重复变化反映图像的灰度统计信息、空间分布信息和结构信息。纹理特征中的灰度共生矩阵和半方差矩阵是具有代表性的两个特征，前者是应用最为广泛的纹理描述方法之一，后者则能有效表示建筑区域的特征。因此本文选择灰度共生矩阵和半方差矩阵作为纹理特征集合。
综上所述，本文选取如表 1所示的特征集合，分别是直接测量得到的极化数据、Yamaguchi分解特征[10]、Cloude-Pottier极化分解特征[11]，还有灰度共生矩阵[12]和半方差函数[13]等纹理特征。
|
表 1 特征集合 Table 1 Feature set |
其中，直接测量得到的极化数据主要有极化散射总功率Span和极化相干矩阵的对角线元素T11, T22, T33。Span表示极化SAR系统的回波强度，T11代表目标对称因子，T22代表目标非对称因子，T33代表目标非规则因子。
Yamaguchi分解将地物目标散射分解成表面散射PS、二次散射PD、体散射PV和螺旋散射PH。螺旋散射信息PH的存在是Yamaguchi分解优于Freeman-Durden分解之处，该成分在具有复杂结构的目标或城区建筑中较为常见。
Cloude-Pottier分解参数包括极化熵H、散射角α与各向异性度A，它们都来自于相干矩阵T的分解，均有旋转不变性。极化熵H用来衡量散射过程的随机性，散射角α与平均物理散射机制密切相关，各向异性度A用来描述特征分解得到的第2、3个特征值的相对大小。
纹理特征选取灰度共生矩阵和半方差函数。灰度共生矩阵包含熵、对比度、同质性、均值4种特征，分别用符号ENT、CON、HOM、MEA表示。ENT描述图像信息量的随机性，CON反映图像灰度的局部差异，HOM体现图像局部均匀程度，MEA体现纹理的规则化程度。半方差函数V[14]描述像素与其邻域像素的空间相关性。将图像中的灰度值作为区域化变量，在以像素点为中心的邻域窗口内对0°、45°、90°、135°等4个方向的半方差函数取平均，最终得到中心像素的半方差纹理特征值。
2 分类方法 2.1 分类器本文构建了由ADA_RF和SVM组成的二级分类结构。
AdaBoost是一种集成学习算法，它能根据分类器的性能赋予权重，得到分类器的加权组合，提升整体分类性能。具体来说便是根据各个分类器的分类正确率高低来分配权重，正确率高的获得高权重，低的则获得低权重。RF[7]是一种集成学习模型，决策树是它的基本分类器。它利用自助采样法随机采集样本组成子集，每个样本子集训练出一个决策树，多颗决策树组成随机森林。各个决策树通过简单投票法决定最终结果。
本文用AdaBoost算法调整随机森林中决策树的投票权重，分类能力越强权重越高，最后通过最大投票准则获得分类结果。AdaBoost算法是面向二分类的分类算法，因此本文采用“一对多”拆解法，将多类别分类问题转化为多个二分类问题。对像素分类时，每颗决策树单独给出分类结果，随机森林对投票结果加权统计，将得票最多的预测结果作为最终分类结果。
SVM是一种典型的有监督分类器[15]。它的基本思想是找到能够正确划分训练数据集并且间隔最大的超平面，该超平面要有足够间隔将离其最近的两类样本点分开，所以可将超平面的求解表示为一个凸二次优化问题
| $ \begin{array}{c} {\min _{w, b}}\frac{1}{2}\left\| \mathit{\boldsymbol{w}} \right\|, \\ {\rm{ s}}{\rm{. t}}{\rm{. }}{y_i}\left( {{\mathit{\boldsymbol{w}}^{\rm{T}}}{x_i} + b} \right) \ge 1, i = 1, 2, \cdots , m{\rm{. }} \end{array} $ |
其中：w为权重系数，b为位移项，(xi, yi)为样本集中的训练样本。SVM具有比较完备的数学理论支撑，但在大样本和高维度场景下计算效率不高。因此可以考虑利用随机森林的特征筛选功能，提前去除冗余特征，提高运行效率。
由ADA_RF和SVM得到的两级分类结果通过邻域投票法进行融合，其原理是当两级分类结果在像素点(i，j)处不同时，以该点为中心取一个宽度为d的正方形窗口，窗口中每个像素点都可以对分类结果进行投票，将中心像素的类别标记为得票最多的类别C
| $ \begin{array}{c} C = {\max _c}\left( {{\mathop{\rm SUM}\nolimits} \left( {{E_n}(k, l) = = c} \right)} \right), \\ c = 1, 2, 3, \cdots ;n = 1, 2. \end{array} $ |
其中：En代表初级分类结果E1、二级分类结果
本文算法流程如图 1所示。
|
Download:
|
| 图 1 本文算法流程 Fig. 1 Flow chart of the algorithm | |
步骤1 输入极化SAR数据，进行预处理。
步骤2 提取数据的极化和纹理特征，构成维度为M的原始特征参数集X。
步骤3 分析极化SAR图像，查找或绘制地面真值图。利用分层采样法选取训练集和测试集。
步骤4 通过输入训练数据集和标签到ADA_RF，使得随机森林的每个决策树都有各自的权重。用训练好的分类器对整幅图像进行预测，得到初级预测结果E1，同时得到各个特征的平均重要性排名V1。
步骤5 根据重要性排名V1，先将重要性最高的特征向量作为特征参数集X1，得到预测结果A1，并计算分类精度，再添加重要性次之的特征向量到X2中，计算得到对应的分类精度A2…依次添加重要性次之的特征向量组成特征参数集合，直到加完所有特征向量。最后得到{X1，X2, …, XM}和{A1，A2, …, AM}组成的折线图，观察折线图，选择使得X的维度N尽可能小且满足
步骤6 利用筛选后的特征集合X训练SVM分类器，并对整幅图像进行预测，得到二级分类结果E2。
步骤7 利用邻域投票法，将E1和E2进行融合，得到最终地物分类结果。
3 实验与分析本文通过对AIRSAR采集自美国旧金山和荷兰弗莱福兰省的L波段极化数据进行分类实验，得到RF、SVM、ADA_RF和本文算法4种分类器的分类结果，通过对比来验证本文算法的有效性。
3.1 特征筛选ADA_RF可以根据各个特征的贡献率得出重要性排名，作为特征筛选的依据。按照步骤5得到重要性排名，依次加入特征并记录精度的提升程度，得到图 2所示的折线图。
|
Download:
|
| 图 2 特征向量个数对分类精度的影响 Fig. 2 Effect of the number of features on classification accuracy | |
根据图 2确定两个实验的特征筛选个数N=10, 即选择重要性排名前10的特征集，分别是Span, PS, PD, PV, PH, H, α, A, ENT, V。
3.2 分类对比实验第1个对比试验利用的是美国旧金山海湾地区的数据。数据大小为900×1 024，PauliRGB图和地面真值图像如图 3(a)、3(b)所示，其中地面真值图像的绘制参考了文献[16]和谷歌地图。将地物类别分为5类，分别是山地、海洋、城市、植被、有角度城市，城市对应的是普通建筑区域，有角度城市对应的是高层建筑区域。
|
Download:
|
| 图 3 旧金山地区数据分类效果图 Fig. 3 Classification images of San Francisco data | |
在实验中，RF和ADA_RF的决策树各设置为200棵，SVM采用的是径向基函数核。训练和测试数据集通过分层采样法获得，二者比例为3∶7，最终分类效果如图 3(c)~3(f)所示，分类精度与耗时数据如表 2所示。
|
|
表 2 旧金山地区数据的分类精度与耗时指标 Table 2 Classification accuracy and time-consuming indexes of San Francisco data |
从分类效果和分类精度数据来看，本文算法的分类精度最高，达到92.72%。尤其是山地和城区的精度相比RF和SVM有较大提升。从分类效果看，本文算法得到的结果相比ADA_RF和SVM结果更加光滑连续，说明邻域投票法可有效抑制相干斑噪声。由图 3(c)、3(d)可知，RF和SVM在区分城市建筑物高度问题上效果欠佳，而图 3(e)、3(f)表明经过AdaBoost改进的随机森林可以更好地区分普通城区和高建筑物城区。在算法的耗时上，RF耗时最短，因为ADA_RF引入了决策树权重的计算，故相比于RF耗时有所增加，但本文算法的耗时小于ADA_RF和SVM耗时的直接加和，表明特征筛选在减少计算量上发挥了作用。
为进一步验证本文算法的有效性，选取Flevoland数据进行二次验证。数据大小为750×1 024，包含15类地物，图 4(b)地面真值图像的绘制参考了文献[17]和谷歌地图。分类器中RF和ADA_RF的决策数仍然设为200棵，训练数据集和测试数据集通过分层采样获取，比例设置为3∶7。分类结果如图 4所示，各项分类指标如表 3所示。
|
Download:
|
| 图 4 荷兰弗莱福兰地区数据分类图 Fig. 4 Classification images of Flevoland-Netherlands data | |
|
|
表 3 弗莱福兰地区数据的分类精度与耗时指标 Table 3 Classification accuracy and time-consuming indexes of Flevoland data |
从分类效果和分类精度上看，ADA_RF在分类精度上相比于RF和SVM有所提升，分别提升6.71和8.31个百分点。本文算法又在ADA_RF的基础上进一步提升了分类精度，达到94.68%。由此可以证明ADA_RF算法在提升分类精度上具有一定效果，经过邻域投票法之后可以有效抑制相干斑噪声，使得分类精度进一步提升。
4 结论本文设计了一种基于ADA_RF和SVM的二级分类结构。通过对旧金山和弗莱福兰两个地区极化SAR数据的实验，结果表示本文所提算法的分类精度高于RF、SVM和ADA_RF。证明了本文算法的有效性，且ADA_RF算法能够根据RF中决策树的分类能力给决策树赋予权重，提升整体分类精度，证明特征筛选之后再训练SVM分类器能够减少计算量与时间消耗。证明通过邻域投票法将两级分类结果融合，能够在保持空间一致性的前提下抑制相干斑噪声的影响。
| [1] |
张澄波. 综合孔径雷达: 原理、系统分析与应用(影印本)[M]. 北京: 科学出版社, 1989: 1-3.
|
| [2] |
Lee J S, Pottier E. 极化雷达成像基础与应用[M]. 洪文, 李洋, 尹嫱, 等译. 北京: 电子工业出版社, 2013.
|
| [3] |
Maitre H. 合成孔径雷达图像处理[M]. 孙洪, 译. 北京: 电子工业出版社, 2005: 124.
|
| [4] |
Chen S W, Tao C S. PolSAR image classification using polarimetric-feature-driven deep convolutional neural network[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(4): 627-631. Doi:10.1109/LGRS.2018.2799877 |
| [5] |
Liu H Y, Min Q, Sun C, et al. Terrain classification with polarimetric SAR based on deep sparse filtering network[C]//2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). July 10-15, 2016, BeiJing, China. IEEE, 2016: 64-67. DOI: 10.1109/IGARSS.2016.7729007.
|
| [6] |
张妙然, 刘畅. 基于特征筛选和二级分类的极化SAR建筑提取算法[J]. 中国科学院大学学报, 2018, 35(1): 89-95. Doi:10.7523/j.issn.2095-6134.2018.01.012 |
| [7] |
Du P J, Samat A, Waske B, et al. Random forest and rotation forest for fully polarized SAR image classification using polarimetric and spatial features[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 105: 38-53. Doi:10.1016/j.isprsjprs.2015.03.002 |
| [8] |
Salehi M, Sahebi M R, Maghsoudi Y. Improving the accuracy of urban land cover classification using radarsat-2 PolSAR data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(4): 1394-1401. Doi:10.1109/JSTARS.2013.2273074 |
| [9] |
陈伟民, 张凌, 宋冬梅, 等. 基于AdaBoost改进随机森林的高光谱图像地物分类方法研究[J]. 遥感技术与应用, 2018, 33(4): 612-620. |
| [10] |
Sun X, Song H J, Wang R, et al. High-resolution polarimetric SAR image decomposition of urban areas based on a POA correction method[J]. Remote Sensing Letters, 2018, 9(4): 363-372. Doi:10.1080/2150704X.2017.1418989 |
| [11] |
Kumar A, Panigrahi R K. Entropy based reconstruction technique for analysis of hybrid-polarimetric SAR data[J]. IET Radar, Sonar & Navigation, 2019, 13(4): 620-626. Doi:10.1049/jet-rsn.2018.5338 |
| [12] |
Gola J, Webel J, Britz D, et al. Objective microstructure classification by support vector machine (SVM) using a combination of morphological parameters and textural features for low carbon steels[J]. Computational Materials Science, 2019, 160: 186-196. Doi:10.1016/j.commatsci.2019.01.006 |
| [13] |
Zhao L J, Zhou X G, Kuang G Y. Building detection from urban SAR image using building characteristics and contextual information[J]. EURASIP Journal on Advances in Signal Processing, 2013, 2013(1): 56. Doi:10.1186/1687-6180-2013-56 |
| [14] |
Dekker R J. Texture analysis and classification of ERS SAR images for map updating of urban areas in The Netherlands[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(9): 1950-1958. Doi:10.1109/TGRS.2013.814628 |
| [15] |
陈刚, 王宏琦, 孙显. 基于核函数原型和自适应遗传算法的SVM模型选择方法[J]. 中国科学院研究生院学报, 2012, 29(1): 62-69. Doi:10.7523/j.issn.2095-6134.2012.1.009 |
| [16] |
Xiao D L, Liu C, W Q, et al. PolSAR image classification based on dilated convolution and pixel-refining parallel mapping network in the complex domain[EB/OL]. arXiv: 1909. 10783. (2020-01-20)[2021-03-08]. http://arxiv.org/abs/1909.10783.
|
| [17] |
滑文强, 王爽, 郭岩河, 等. 基于邻域最小生成树的半监督极化SAR图像分类方法[J]. 雷达学报, 2019, 8(4): 458-470. |