随机森林协同Sentinel-1/2的东营市不透水层信息提取

时间：2024-07-28

刘春亭，冯权泷，金鼎坚，史同广，刘建涛，朱明水

(1.山东建筑大学测绘地理信息学院，济南 250101；2.中国农业大学土地科学与技术学院，北京 100083；3.中国自然资源航空物探遥感中心，北京 100083)

0 引言

不透水层是指诸如沥青或水泥道路、屋顶、停车场等阻止水下渗的天然或人工地表覆盖[1]。当前城市化正处在快速发展的阶段，随着城市化进程的不断加快，城市空间扩张和人口的快速增长，城市热岛、环境恶化、洪涝灾害、交通拥堵等一系列问题愈发严重，给生态环境及城市发展带来巨大挑战。城市化进程的显著表现之一就是不透水层的扩张。在城市化过程中，土壤、水域、植被等自然景观被大范围的不透水层替代，阻碍了地表的蒸散作用，破坏了地表热平衡，城市热岛效应使得城市增温显著。不透水层的广泛存在还导致水循环的失衡，造成地下水补给减少、水质下降、水污染、暴雨流量和洪水内涝灾害频发[2-3]，严重影响人类生产生活。研究表明，不透水层的面积和分布是城市化对环境造成影响的重要因子[4]，因此及时准确地获取不透水层信息对于城市化及生态环境保护具有重要意义[5]。

遥感技术由于大范围同步观测、重复周期短、动态、迅速的优势，成为从局部到全球尺度估算不透水层的主要技术手段。近年来不同数据源、不同尺度下提取不透水层的研究已经取得了很多成果，利用遥感数据提取不透水层的方法主要有混合光谱分析[6-8]、指数模型[9-11]、基于决策树[12-13]、支持向量机[14-17]、神经网络[18-19]等方法的影像分类，以及多源遥感数据结合的方法[20-22]等。但由于城市景观的复杂性，准确及时地估算不透水层空间分布信息仍然是一项具有很强挑战性的任务[23-25]，如建筑物及树木的阴影容易与反射率较低的不透水层相混淆，裸土、盐田等也易错分为不透水层。

通过文献综述发现目前关于黄河三角洲地区的不透水层的信息提取多依靠光学遥感影像。随着遥感技术的发展，多平台、多传感器、高空间分辨率、高光谱分辨率的遥感数据不断涌现，融合多源数据被认为是改善不透水层提取精度的一种有效方法。光学遥感容易受到云、雨天气的影响，而雷达数据具有全天时、全天候的优势，能够获得稳定的周期性数据，可以弥补光学影像特殊天气情况下无法及时监测的缺陷[26-28]。协同使用光学和雷达数据进行地物识别已经成为近年来的研究热点之一，部分研究表明使用合成孔径雷达遥感数据补充光学遥感数据，可在一定程度上提高不透水层的提取精度[29-32]。

近几十年来，石油开采[33]、旅游开发[34]、人工建筑[35-36]等人类活动导致黄河三角洲地区生态系统发生明显改变，生态环境逐渐恶化，引起国内外学者的广泛关注。但是通过文献调研发现结合主被动遥感对黄河三角洲地区进行不透水层提取，尤其是高空间分辨率不透水层提取的研究较少，因此本研究以东营市为研究区域，探索机器学习方法协同Sentinel-2光学数据和Sentinel-1雷达数据提取不透水层信息的有效性。具体来说，本文协同使用地表反射率、纹理、后向散射系数等特征，利用随机森林算法对研究区不透水层进行提取，探究和定量分析光学遥感与主被动遥感结合在信息提取中的精度差异，为黄河三角洲区域遥感监测、开发与保护提供科学依据。

1 研究区概况与数据源

1.1 研究区概况

黄河三角洲有着我国暖温带地区最年轻、最广阔、生物多样性最丰富的湿地生态系统，在维护区域生态安全和可持续发展方面有重要意义[37]。根据《黄河三角洲高效生态经济区发展规划》中划定的黄河三角洲范围，本文以山东省东营市为研究区域(图1)。东营市是黄河三角洲地区的中心城市，地理位置横跨N36°55′～38°10′，E118°07′～119°15′，属于典型的温带大陆性季风气候。其东部和北部濒临渤海，境内河流较多，水资源丰富，黄河在东营市垦利县流入渤海，泥沙淤积现象严重，后备土地资源量大，开发空间充足；土地类型多，可利用性较强。受河流、海洋的影响，盐碱地分布的范围较大，土壤发育不稳定，在一定程度上影响了当地的发展。东营市有丰富的石油资源、旅游资源，其独特的环境条件使得该地区土地利用/覆盖变化剧烈。

图1 研究区Sentinel-2 B4(R),B3(G),B2(B)波段合成影像Fig.1 Image of Sentinel-2 B4(R),B3(G),B2(B)bands in the study area

1.2 数据源及预处理

本文所使用的影像来自欧州航天局“哥白尼计划”地球观测卫星系列的Sentinel系列卫星，数据通过欧州航天局数据共享网站(https://scihub.copernicus.eu/dhus/#/home)获得。为了提高不透水层与其他地物之间的可分性并结合研究区影像质量状况选取了2019年9月份的遥感影像开展不透水层的提取研究。

Sentinel-1号卫星星座拥用A和B 2颗卫星，可以达到每6 d对同一地点进行重复观测。本研究中选用的2景Sentinel-1数据为雷达干涉宽刈幅(interferometric wide swath，IW)成像模式的Level-1级别地距多视影像(ground range detected，GRD)，极化模式为VV和VH，影像的获取时间为2019年9月26日，空间分辨率10 m。

Sentinel-1数据预处理主要在SNAP Desktop软件中完成，处理步骤主要包括轨道校正、热噪声去除、辐射定标、斑点滤波、地理编码、影像拼接等。处理的具体流程为：①使用精密轨道文件进行细化，以获取准确地卫星位置和速度信息，提高影像配准和地理编码的精度；②根据元数据集提供的噪声查找表消除雷达系统内在的热噪声；③雷达数据存在辐射偏差，通过辐射定标后使雷达影像的像素值真正表示后向散射系数σ0，辐射定标公式为：

(1)

式中：DNi与Ai分别为像元i的灰度值与定标参数；④斑点噪声随机分布于影像中，使用Refined Lee滤波器进行斑点滤波，在降噪的同时能够更好地保留辐射和纹理信息[38-40]；⑤地理编码是将投影转换为地理坐标投影，并且纠正雷达成像过程中由于透视收缩、叠掩等引起的几何形变，采用30 m的SRTM DEM改正地物在雷达影像中发生的畸变，输出10 m空间分辨率的重采样数据。

Sentinel-2号卫星星座同样也是由A和B 2颗卫星组成，重访周期为5 d。本文选用的2景Sentinel-2数据是经过大气校正的地表反射率产品L2A，获取时间为2019年9月29日。Sentinel-2 MSI数据包含13个波段，具体参数如表1所示。在开展实验时，本文仅使用了10 m和20 m空间分辨率的波段。

表1 Sentinel-2波段信息Tab.1 Band information of Sentinel-2

为保证不同空间分辨率数据在空间上的匹配性，本文将20 m空间分辨率的数据重采样至10 m，且所有数据都统一转换为WGS-84 UTM 50N投影。

2 研究方法

不透水层信息提取的技术流程如图2所示，主要包括4个步骤：①影像下载及预处理，包括Sentinel-1影像获取后进行轨道校正、辐射定标、斑点滤波等处理得到后向散射系数，Sentinel-2影像重采样以及影像配准、镶嵌、裁剪等预处理工作；②特征选择，基于灰度共生矩阵计算纹理特征，协同地表反射率和后向散射系数构建多维特征空间；③确定分类体系及样本选择；④由选好的训练样本训练随机森林，执行分类；⑤根据验证样本得到混淆矩阵进行精度验证，并与不同特征组合、不同方法进行精度对比。

图2 技术流程Fig.2 Flowchart of research

2.1 特征提取

已有的研究表明纹理特征有助于改善地表信息提取的精度[41-44]。本文使用基于灰度共生矩阵(gray level co-occurrence matrix，GLCM)的纹理特征，选取了均值(MEA)、方差(VAR)、同质性(HOM)、相异性(DIS)、熵(ENT)、角二阶矩(ASM)6种彼此相关性弱的纹理特征，计算公式分别为：

(2)

(3)

(4)

(5)

(6)

(7)

式中：N为灰度级数目；P为N×N维归一化灰度共生矩阵；P(i,j)为矩阵中第i行、第j列的标准化灰度值。计算处理窗口设置为7像素×7像素[45]。另外考虑到光谱波段之间的相关性，为减少计算量，只使用红光波段计算纹理特征。

将纹理特征、后向散射系数与反射率特征相融合，共同构造不透水层提取的输入特征，即包含10个反射率特征波段、6个纹理特征波段、2个雷达数据的后向散射系数波段(VV和VH)。

2.2 样本选择

样本的选择直接影响识别模型的稳定性和分类精度。使用ENVI软件中ROI提取工具建立感兴趣区实现样本的选择，结合研究区的特点，将其土地利用/覆盖类型划分为：亮不透水层、暗不透水层、有作物耕地、空闲耕地、大棚用地、林地、水域、滩涂、盐田、未利用地共10类，其分类标准及数量如表2所示，训练样本和验证样本空间分布如图3所示。

表2 分类体系及样本数量Tab.2 Classification scheme and number of samples (个)

图3 样本点分布Fig.3 Distribution of sampling points

2.3 随机森林算法理论

随机森林算法是Breiman于2001年提出的[46]，目前已经广泛应用于遥感应用研究[47-50]。随机森林算法对数据的适应能力强，数据集无须规范化，随机森林分类受噪声和异常值的影响较小，且不会过拟合[46]，能够直接处理高维数据，对于缺省值也可以获得较好的结果。该算法的基本思想是首先利用Bootstrap方法从原始训练集中有放回的随机抽取n个样本并构建n个决策树构成森林，n个样本的个数约为总样本个数的2/3。假设训练样本数据中有m个特征，每次分裂时根据最优特征进行分裂，每棵树都充分生长，直到该节点的所有训练样本都属于同一类，让每棵决策树在不做任何修剪的前提下最大限度无限生长，最后将生成的多棵分类树组成随机森林，根据多棵树分类器投票决定最终分类结果。其中每次随机抽取的n个样本小于样本量总数用于替换训练样本，其余用来进行交叉验证，以此来提高随机森林方法的泛化能力。在使用随机森林分类时，通常需要设置生长树的数量和节点分裂时输入的特征变量数量。经过一系列实验发现，当N=500时，袋外误差(out of bag，OOB)趋于稳定，因此选择500棵树用于随机森林分类，特征变量数则由参与分类的特征数的平方根计算得到[51]。

2.4 精度评估

分类完成后，采用混淆矩阵对分类结果进行精度验证，由混淆矩阵派生的评价指标主要有总体分类精度、制图精度(producer’s accuracy，PA)、用户精度(user’s accuracy，UA)和Kappa系数。根据目视解译获得验证样本，在ENVI软件中利用Confusion Matrix Using Ground Truth ROIs 工具进行精度验证。

3 结果与分析

3.1 分类结果及精度验证

采用本文方法分类结果的混淆矩阵如表3所示，分类结果总体精度达到93.37%，Kappa系数达到0.925 8。从表3中可见，空闲耕地、滩涂分别有11个和9个样本点错分为亮不透水层，这些土地利用类型表面存在的盐碱颗粒与不透水层光谱性质相似，造成了错分，导致不透水层的高估现象。研究区内的未利用地主要是油田附近被压实的裸土以及正在建设的区域，因此未利用地易与亮不透水层或暗不透水层混淆，其中有36个未利用地样本点错分为亮不透水层，13个暗不透水层错分为未利用地。有作物耕地与林地也发生了混淆现象，其中68个林地样本被错分为有作物耕地，同时21个有作物耕地错分为林地，这类错分可能与样本的选择有关，同时也受到林地密度、树木高度及季节等因素的影响。

表3 分类结果混淆矩阵Tab.3 Confusion matrices for classification

将分类结果的暗不透水层和亮不透水层合并为不透水层，其他地类合并为非不透水层，得到东营市不透水层分布情况(图4)。结合多光谱影像及Google Earth影像发现建城区、农村居民点、道路等不透水层与现实分布状况吻合，取得了较好的提取效果。总体来看，研究区内不透水层多为暗不透水层，主要分布在东营区城区、垦利区城区，河口区、广饶县、利津县的建成区和分散式分布的农村居民地及道路。

图4 不透水层提取结果Fig.4 Extracted result of impervious surface

3.2 不同特征组合分类结果比较

为探索Sentinel卫星的光谱和雷达数据对不透水层提取精度的影响，本文设计了3种组合方案进行对比(表4)。

表4 Sentinel数据组合方案Tab.4 Combination schemes of Sentinel data

通过对上述3种特征组合方案进行随机森林分类(表5,图5)，发现仅使用后向散射系数(方案B)分类时精度较低，总体精度为40.76%，Kappa系数为0.335 8。使用地表反射率特征及纹理特征(方案A)进行随机森林分类时已经可以得到较为理想的分类精度(总体精度为92.04%、Kappa系数为0.911 0)，而且暗不透水层的PA和UA均大于90%，亮不透水层的PA和UA也达到87%以上。在方案A基础上引入后向散射系数辅助特征后，总体精度和Kappa系数分别达到93.37%和0.925 8，暗不透水层的PA和UA分别提高了1.80和2.72百分点，亮不透水层的UA提高0.47百分点，说明协同Sentinel-1和Sentinel-2数据可以有效改善不透水层的估算精度。

表5 分类精度统计Tab.5 Statistics of classification accuracy

(a)方案A (b)方案B (c)方案C

由于研究区域较大且使用10 m空间分辨率，图5中不同特征组合方案分类结果差异不太明显，尤其是方案A和C。因此结合混淆矩阵发现，后向散射系数对各类地物的提取精度都不高，但能提供光谱数据的互补信息，使得方案C分类精度有所提高。结合后向散射系数后，除大棚用地和盐田的分类精度在正常波动范围内略有降低外，其他地类提取精度均有不同程度的提高。数据组合A方案中，暗不透水层、未利用地和空闲耕地的混淆程度较高，有作物耕地与林地也存在一定程度的混淆。在方案C中这种情况被明显改善，后向散射系数提高了不透水层与空闲耕地、未利用地等裸露土地的区分度，地物错分减少，提取精度获得提高。

同时可以看到，后向散射系数(方案B)的分类结果除了精度较低以外，还有明显的椒盐现象，图像上存在大量的离散斑块和孤立像元。基于光谱和纹理(方案A)的分类结果椒盐现象有所缓解，孤立像元数据有所减少。光学、纹理和雷达结合(方案C)的分类结果中分类精度提高，地物边界分类清晰，椒盐现象被较好抑制。

总的来说，利用Sentinel-2地表反射率、纹理特征和Sentinel-1雷达影像后向散射系数能够提高不同地物类型的可分性，尤其在减少不透水层与裸露土地的混淆程度上有显著作用，并且对于基于像素分类常出现椒盐现象有明显抑制作用，因此这种特征组合在获取高空间分辨率不透水层信息方面有较大的优势。

3.3 与其他方法比较

为了进一步评估随机森林方法在提取不透水层的优势，本文选取了支持向量机、决策树等遥感图像分类领域较常用的机器学习方法进行了对比实验。为了保证不同方法所得结果的可比性，本文使用了相同的数据(方案C)和样本，3种方法的比较结果如表6所示。

表6 随机森林、支持向量机、决策树分类精度比较Tab.6 Comparison of classification accuracy of RF,SVM and CART

从表6中可以看出，随机森林分类的精度最高，支持向量机分类精度略低(总体精度为93.19%，Kappa系数为0.923 8)，决策树分类精度最低(总体精度为87.79%，Kappa系数为0.863 5)。随机森林与支持向量机分类方法精度都较好，体现了Sentinel数据对于不同方法的良好适用性。但是与随机森林分类算法相比，支持向量机分类方法需要设置多个参数来执行分类，计算量大，耗费时间长，这也正是随机森林得到广泛应用的原因之一。

由此可见，特征集合与分类方法是影响不透水层提取精度的关键因素。基于多源Sentinel数据，运用随机森林算法提取不透水层信息可以获得较高的精度，说明Sentinel卫星光谱数据与雷达数据结合在提取不透水层提取方面的应用潜力，也验证了随机森林算法在分类方面具有的优势。