1. 引言
山西省吕梁市柳林县是山西省滑坡发生最多的地区 [1]，达70多处，占山西省滑坡总数的31.53%。目前，越来越多的学者倾向于利用高空间、高光谱及高时间分辨率的遥感数据解译、识别大区域、大面积滑坡时空分布特征及规律 [2]。然而针对山西省黄土高原滑坡提取的研究还比较少，且主要的提取方法是实地调查。目前针对山西省黄土滑坡的研究有：周潇朗 [3]，乔鹏腾等 [4] 选取了吕梁山区较为典型的滑坡分析了滑坡形成机理；欧阳林辉等 [5] 分析了山西省临汾至吉县高速公路乡宁_吉县段地质病害特征；王昌明 [6] 进行了山西吕梁黄土滑坡成灾模式及降雨预警模型研究，通过野外实际调查记录了柳林县25处典型的黄土滑坡；胡胜等 [7] 通过无人机和三维激光扫描仪等先进测量技术开展了黄土高原滑坡野外调查，遗憾的是此次调查在柳林县仅采取两个样本。相关研究表明，实地调查需要进行大量的野外试验，然而野外试验更多的是收集一些验证样本，不适合大区域的滑坡提取，基于遥感影像的滑坡提取可以弥补这一缺陷，节省大量的时间和人力 [8]。目前基于遥感影像进行柳林县黄土滑坡提取的实验还比较少，为丰富相关领域的研究，本文选择了黄土高原滑坡提取时，最佳波段组合这一不确定性问题进行研究。
研究选用地类提取常用的遥感影像Landsat-8 OLI [9] 和新兴的用于灾害支援的Sentinel-2A [10] 数据进行实验。波段组合实验在土地利用类型提取中展开的较多，李秀梅等 [11] 经过光谱统计，提出波段742为渤海湾海岸线变化信息提取的最佳波段组合；晁增福 [12] 经过最佳指数模型(optimum index factor, OIF)计算，发现波段457为新疆维吾尔自治区阿拉尔市土地利用分类的最佳波段组合；周旭等 [13] 在光谱统计和OIF计算的基础上，对研究区的地类绘制了光谱特征曲线，为最佳波段的选择拓宽了研究思路。现阶段，针对黄土滑坡提取研究中通过波段组合增强影像质量的实验还比较少，研究从光谱特征统计、最佳指数因子和光谱特征曲线出发，分析遥感影像提取黄土滑坡的最佳波段组合。
2. 研究区及数据
2.1. 研究区概况
柳林县地处山西省吕梁市西部，西邻黄河，东靠吕梁山，坐标范围在东经110˚37'~111˚6'，北纬37˚9'~37˚39'之间(图1)，境内皆被黄土覆盖，地形破碎，坡陡沟深，平地较少，黄土地貌主要为黄土梁谷和黄土塬谷 [14]。柳林县内水系发育，跨越多条河流，地质构造复杂，地层岩性主要为泥岩、砂岩、粉质粘土、奥陶系灰岩和松散的黄土堆积层。特殊的黄土地貌，复杂的地质构造、极其脆弱的地质环境，导致区域内滑坡发育率较高。
柳林县是山西省滑坡发生频次最高的地区 [15]，研究根据文献描述，选择柳林县境内滑坡空间分布相对集中的区域作为研究区，坐标范围在东经110˚51'~110˚55'，北纬37˚26'~37˚35'之间，该区域也为黄河河岸带向吕梁山过度的区域，满足滑坡发生的高程、坡度等基本的孕灾条件 [16]。
Figure 1. Location map of the study area
图1. 研究区位置图
2.2. 数据及预处理
研究采用的影像数据为谷歌2 m空间分辨率数据(免费来源于谷歌地图下载器BIGEMAP http://www.bigemap.com/)、Sentinel-2A 10 m多光谱数据(免费来源于欧空局http://scihub.copernicus.eu/)、Landsat-8 OLI多光谱数据(免费来源于中国地理空间数据云http://www.gscloud.cn/sources/)。其中，谷歌2 m数据只包含RGB三个波段，作为滑坡验证样本数据来源；哨兵2A数据共包含13个波段，其中红、绿、蓝和近红外4个波段为10 m空间分辨率，为了保持数据空间分辨率的一致性，研究将20米分辨率的11、12波段在SNAP中上采样为10米空间分辨率 [17]；Landsat-8 OLI数据选用30米空间分辨率的波段。所选数据都是经过几何精校正的影像产品，研究只对Landsat-8 OLI数据进行大气校正处理。卫星影像及其参数见表1。
研究总体技术路线如图2所示，研究主要包括三个方面，首先是对Landsat-8 OLI和Sentinel-2A做基础处理，包括大气校正、拼接和裁剪等；其次是进行最佳波段选择实验，从单波段标准差、多波段相关系数、最佳指数因子以及地类光谱特征曲线四个方面进行分析；最后对最佳波段组合的结果进行验证，分别进行目视判别的定性验证和监督分类的定量验证。
Table 1. Satellite images and their parameters
表1. 卫星影像及其参数
3. 滑坡提取最佳波段选择实验
3.1. 选取最佳波段的必要性
为了获得更好的滑坡目视解译效果，研究根据遥感的滑坡信息提取应用，对影像进行图像增强处理，增强的目的是将滑坡突出显示。目前文献中多是通过计算NDVI [18]、波段加权等 [19] 波段计算突出植被的光谱信息，增加植被与非植被区的对比差异，以此来增加滑坡的识别度。然而，波段合成是滑坡灾害遥感解译中一项关键的不确定性问题 [20]，对于提取黄土滑坡这一特殊的地质信息，正确的波段合成可以提高滑坡的识别度。本研究重点对Sentinel-2A和Landsat-8 OLI数据进行波段合成实验，探讨适于黄土滑坡提取的最佳波段组合。王治华教授 [21] 在讨论中国滑坡遥感及进展时指出，在基于遥感影像进行滑坡提取时，有必要进行多波段合成，选取合成波段时尽可能遵循两个原则，一是选取的波段之间的相关性达到相对最小，二是选取能反映岩石特征的波段。基于该原则，研究需要进行波段光谱特征分析、波段组合最佳指数因子分析、滑坡等地物光谱曲线分析等，逐步筛选滑坡提取最佳波段。
3.2. 波段光谱特征分析
3.2.1. 单波段标准差分析
为了保持数据波段的一致性，研究根据Sentinel-2A和Landsat-8 OLI数据波段的特征，同时选取两种数据的Blue、Green、Red、NIR、SWIR1和SWIR2，7个波段。分别计算两种数据各波段的标准差并对结果进行排序，结果见表2。结果表明，两类数据的排序结果一致，均为：SWIR1 > SWIR2 > NIR > Red > Green > Blue，可知波段SWIR1所含信息量最大，波段Blue所含信息量最小。依据所选波段信息量最大的原则，应优先选择SWIR1波段，Landsat-8 OLI对应波段为6波段，Sentinel-2A为11波段。然而，根据本文滑坡识别的主要目的以及王治华教授提出的原则，结合Landsat-8 OLI中SWIR2主要用于识别岩石以及矿物 [22]，Sentinel-2A中SWIR2主要用于地质监测 [23]，以及SWIR2标准差仅次于SWIR1的排序结果，选择Landsat-8 OLI的波段7，Sentinel-2A的波段12用于研究。
Table 2. Spectral characteristics statistics of single band of remote sensing data in the study area
表2. 研究区遥感数据单波段光谱特征统计
3.2.2. 多波段相关系数分析
多波段间相关系数矩阵表明(Landsat-8 OLI表3，Sentinel-2A 表4)，Landsat-8 OLI中Blue、SWIR1、SWIR2，Sentinel-2A中NIR、SWIR1、SWIR2，与各波段间的相关度比较低，若单以波段间相关系数为筛选标准，则这几个波段可作为候选波段，但分析发现呈现的相关系数整体偏高，反映出这些波段之间的信息重复率较高 [24]，不能将地物明显区分，需要进一步深入分析，加以筛选。
Table 3. Multi-band spectral correlation coefficients of Landsat-8 OLI data
表3. Landsat-8 OLI数据多波段光谱相关系数
Table 4. Multi-band spectral correlation coefficients of Sentinel-2A data
表4. Sentinel-2A数据多波段光谱相关系数
3.3. 最佳指数因子分析
3.3.1. 最佳指数因子法
最佳指数因子(Optimum Index Factor, OIF)是为了定量评估波段组合后图像所包含信息量大小而提出的指标。该指数可由单波段标准差和多波段相关系数计算而得，OIF值越大，表明波段组合后的图像所含的信息量越多，即波段组合的效果越好。OIF数学表达式为，
(1)
式中：SDi表示第i个波段的标准差；Rij表示第i、j两波段的相关系数。
3.3.2. 最佳指数因子分析
为了减小计算范围，根据上述光谱特征分析结果，以及选择能反映岩石特征光谱段的原则，选择含有SWIR2 (Landsat-8 OLI，波段7；Sentinel-2A，波段12)的波段组合计算OIF，各组合OIF排序结果见表5 (Landsat-8 OLI)、表6 (Sentinel-2A)。Landsat-8 OLI计算结果显示，排序前4的波段组合均含有6、7波段，然而，多波段间相关系数结果表明，Landsat-8 OLI中6、7波段的相关系数高达0.99，光谱信息重叠度非常高，所以研究中不选择包含6、7波段的组合方案。从剩下的组合中看，波段4-5-7组合的OIF值最大，最适合进行遥感信息的提取。Sentinel-2A计算结果显示，排序前5的波段组合(除4-8-12)均含有11、12波段，同样，Sentinel-2A中11、12波段的相关系数高达0.99，应予排除，其余组合中4-8-12组合的OIF值最大，表明该组合为最佳候选方案。但由于研究是针对黄土高原滑坡的提取，所选取的波段是否易于区分滑坡，仍需进一步分析。研究进一步绘制地类光谱特征曲线，以筛选出更适合黄土高原滑坡提取的波段组合方案。
Table 5. OIF of Landsat-8 OLI Remote Sensing Image
表5. Landsat-8 OLI数据OIF
Table 6. OIF of Sentinel-2A Remote Sensing Image
表6. Sentinel-2A数据OIF
3.4. 基于地类光谱特征曲线的最佳波段选择
3.4.1. 地类光谱曲线的绘制
为了实现多数地物或者某些地物容易区分，结合研究区实际的土地利用情况，选取7种地物分别计算各类用地在波段上的光谱响应值，各类样本数量为滑坡25、耕地30、城乡住宅50、工矿用地14、交通用地16、其他建设用地6、山地22。各类地物的光谱特征曲线计算结果见下图(Landsat-8 OLI图3，Sentinel-2A图4)。
Figure 3. 7 Kinds of land use type spectrum characteristic curve of landsat-8 OLI
图3. Landsat-8 OLI 7种土地利用类型的光谱特征曲线
Figure 4. 7 Kinds of land use type spectrum characteristic curve of Sentinel-2A
图4. Sentinel-2A 7种土地利用类型的光谱特征曲线
3.4.2. 基于土地利用类型光谱曲线的最佳波段选择
图2表明，Landsat-8 OLI数据中7种土地利用类型的光谱响应值在2-6波段均呈上升趋势，在6波段处形成明显的拐点，6-7波段呈下降趋势，其中其他建设用地在所有波段均可以明显区分；7种土地利用类型的光谱响应值在2、3、4波段上走势相似，且不易区分；到5波段差异初显，其中耕地、城乡住宅、其他建设用地以及工矿用地已经可以明显区分，山地、滑坡和交通用地也得到区分；6波段，除了山地和滑坡不能区分，其余用地均可以得到区分；7波段，滑坡和山地得到了区分。由上述分析可判定，2、3、4波段地类光谱响应差异不大；4、6波段上，山地、滑坡和工矿用地的区分性不佳；5、7波段，地类差异最大，可分性好，为Landsat-8 OLI数据提取黄土滑坡的最佳候选波段。
图3显示，Sentinel-2A 数据中7种土地利用类型的光谱响应值在2、3波段上不易区分；在4波段地类差异比较明显，7种土地利用类型完全可以区分，其中，工矿用地在2-4波段呈下降趋势，4波段后呈上升趋势，4波段处形成拐点；8波段处7种地类的区分性较好；到11波段，山地、交通用地和滑坡难以区分，7种土地利用类型在11波段处均呈明显的拐点；在12波段，7种土地利用类型均可以得到区分。通过上述分析认为，Sentinel-2A数据中，2、11波段地类差异不大，可分性最差；3波段处，除了工矿用地和城乡住宅不能区分，其余用地类型可以得到区分；4、8、12波段，7种土地利用类型数值差异明显，可以区分大部分地类，可分性最好，优先考虑组成最佳波段进行滑坡提取。
3.5. 滑坡提取的最佳波段选择结果
根据最佳波段选择原则，对于Landsat-8 OLI数据，单波段标准差分析结果显示，信息量较大的波段7为候选的最佳波段；多波段相关系数表明2、6、7波段为最佳候选波段；最佳指数因子结果显示4-5-7波段组合的OIF值最大；地类光谱特征曲线表明5、7波段为最佳候选波段；综合上述分析结果，确定含有2、5、7波段的组合为最佳波段。对于Sentinel-2A数据，单波段标准差分析结果显示，信息量较大的波段12为候选的最佳波段；多波段相关系数结果显示8、11、12波段为最佳候选波段；最佳指数因子结果显示4-8-12波段组合的OIF值最大；地类光谱特征曲线表明4、8、12波段为最佳候选波段；综合上述分析结果，确定含有4、8、12波段的组合为最佳组合。
4. 最佳波段组合验证
4.1. 基于目视判别效果的定性验证
上述分析内容表明波段2、5、7 (Landsat-8 OLI)、波段4、8、12 (Sentinel-2A)是最适合研究区滑坡提取的波段，对所有波段组合情况进行全部解析，组合后影像如图5 (Landsat-8 OLI)、图6 (Sentinel-2A)所示。
经过目视判别，Landsat-8 OLI 6种波段组合结果中，波段组合2-5-7、2-7-5、5-2-7以及7-2-5与实际地物的颜色差异较大，应予排除；波段组合5-7-2的建筑物偏绿，不易与植被区分；波段组合7-5-2的颜色层次分明，建筑用地可以和滑坡、耕地等得到明显的区分，耕地和山地也可以明显区分开来，利于滑坡的提取，故认为波段7-5-2的组合顺序为Landsat-8 OLI数据进行黄土滑坡提取的最佳波段组合顺序。
图5显示，Sentinel-2A 6种波段组合结果中，波段4-12-8、8-4-12、8-12-4和12-4-8组合结果的颜色与实际地物的颜色差异较大，给予排除；波段4-8-12组合的结果显示建筑物呈现蓝色，滑坡、山地颜色也偏蓝，建筑物和山地、滑坡等无法得到明显的区分；波段12-8-4组合结果显示滑坡呈红褐色，山地呈绿色，建筑用地呈蓝色，各类用地类别与真实地物最接近，适合滑坡的提取，认为波段12-8-4的组合顺序为Sentinel-2A数据进行黄土滑坡提取的最佳波段组合方案。
4.2. 基于监督分类滑坡提取的定量验证
为了定量分析所选取的方案是否为适合黄土滑坡提取的最佳波段组合，研究选用支持向量机(support vector machine, SVM)分类方法对Landsat-8 OLI和Sentinel-2A两种数据的最佳波段组合进行监督分类的
(a) B2-5-7 (b) B2-7-5 (c) B5-2-7 (d) B5-7-2 (e) B7-2-5 (f) B7-5-2
Figure 5. Images of 6 Kinds of Landsat-8 OLI of Optimal Band Combination
图5. Landsat-8 OLI 6种最佳波段排列组合结果
(a) B4-8-12 (b) B4-12-8 (c) B8-4-12 (d) B8-12-4 (e) B12-4-8(f) B12-8-4
Figure 6. Images of 6 Kinds of Sentinel-2A of Optimal Band Combination
图6. Sentinel-2A 6种最佳波段排列组合结果
滑坡提取，滑坡提取结果如图7所示。Landsat-8 OLI和Sentinel-2A监督分类结果均显示，滑坡主要分布在贯穿研究区南北的“三大线–柳结线”公路两旁，该区域也是工矿用地和城乡住宅的聚集区，分类结果中将很大部分的工矿用地和城乡住宅划分为了滑坡，滑坡的错分误差较大。在ENVI中对两类数据的滑坡提取结果做误差分析，发现Landsat-8 OLI滑坡提取的用户精度为64.77%，Sentinel-2A滑坡提取的用户精度为50%。
从分类结果中滑坡所在的空间位置来看，滑坡与工矿用地、城乡住宅用地的混淆度较高，通过解译2 m空间分辨率的谷歌影像，发现位于交通主干线两旁的部分滑坡已经经过人工处理，不再具有滑坡的光谱特征，被错分为工矿用地和住宅用地。结果表明滑坡多发生在人类工程活动较为频繁的地区，该结果从侧面反映黄土高原固有的岩土物化结构不会轻易导致滑坡发生，相反人类工程活动对滑坡的贡献度更大。
Figure 7. Landslide extraction map of the study area
图7. 研究区滑坡提取图
5. 小结
柳林县是山西省黄土滑坡的高发区，分析发现其境内复杂的地形地貌条件是导致柳林县黄土滑坡难以提取的主要原因。研究基于山西省柳林县滑坡高发区Landsat 8 OLI和Sentinel-2A影像数据，通过分析单波段标准差、多波段光谱特征、最佳指数因子排序以及地物光谱特征曲线逐步筛选满足条件的最佳波段组合方案。经过定性和定量验证，发现波段7-5-2 (Landsat 8 OLI)、波段12-8-4 (Sentinel-2A)为适合研究区黄土滑坡提取的最佳波段组合顺序，同时，选取的波段均满足相关性相对最小、反映岩石信息以及滑坡与山地、耕地易于区分的最佳波段选择原则。
基金项目
国家自然科学基金青年基金(41701537)；自然资源部地理国情监测重点实验室开放基金(2020NGCM07)。