基于GF-2卫星影像的冲积扇地貌识别方法

周 炜，关洪军

(陆军工程大学 野战工程学院， 南京 210007)

微地貌单元是重要的战场地理要素之一，也是影响景观格局的重要因素[1]。微地貌通常指规模相对较小的地貌形态，包含丰富的覆被属性信息、水源信息、地质信息等，对分析军事地质和特殊目标具有重要意义。近年来，微地貌的提取技术主要基于数字高程模型数据(Digital Elevation Model)实现。麻莉[2]基于DEM数据提出了一种针对洪泛平原湿地微地貌单元的数字分析方法。Baade等[3]应用TerraSAR-X雷达数据监测秘鲁南部海岸沙漠地区微地貌变化，比较了利用DEM提取微地貌的局限性。冲洪积扇微地貌广义上讲包括冲积扇与洪积扇两种。冲积扇沉积体具有构造“ 敏感性”，其复杂的沉积机制是解释盆地边缘构造活动的钥匙，也是“源-聚-汇” 系统中的“源” 位置[4]。其前沿多为潜水溢出带和浅埋地下水分布区，对野战给水保障具有重要意义。目前，冲洪积扇微地貌提取技术多基于遥感影像的光学特征，DEM多用于辅助筛选泥石流分布区的沟谷和纵坡分析[5]，这些微地貌识别方法缺少地学因子的参考，不适应复杂多变的微地貌边界的识别。本文以青藏高原西部，北纬32°26′，东经79°36′为研究区域，提出了一种基于CART分类器的冲洪积扇提取模型。该区域全年降雨量少，季节性强，气候寒冷干燥[6]。寒冷、昼夜温差大导致冰川发育广泛，形成了高山、沟谷、冰蚀、冰碛、冲积扇等地貌类型，为冲洪积扇提取提供了大量的实验对象。该模型以GF-2影像、Landsat 8影像、DEM影像为数据源，结合冲洪积扇微地貌的光学特征、几何特征、纹理特征、地学特征位置特征，能够准确界定冲洪积扇扇缘，实现冲洪积扇高精度提取。

1 数据处理

数据源包括Landsat8影像、GF-2卫星影像、数字高程模型。其中，数字高程模型根据1∶5万地形图生成；Landsat8影像数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)；GF-2影像为1 T级产品，参数如表1所示。

表1 GF-2卫星数据波段参数

数据处理包括几何校正、滤波降噪、图像融合3部分。几何校正主要针对GF-2遥感影像、DEM影像、Landsat 8影像进行图像配准，依靠重采样对图像亮度值进行插值计算，建立新的图像矩阵，结果误差控制在0.5个像元以内。滤波降噪通过中值滤波器实现，滤波窗口大小为5×5。中值滤波能够有效滤除椒盐噪声，尽量保护边缘信息的完整[7]。图像融合的主要目的在于借助GF-2全色段影像的高分辨率优势，提高洪积扇边缘的识别精度。融合采用主成分分析法，融合效果如图1所示。

2 冲洪积扇特征分析

冲洪积扇具有独特的地学特征和影像特征，深入研究这些特征因子是构建冲洪积扇提取模型的必要前提。

1) 地学特征分析

① 冲洪积扇顶位于河流出山口，地学特征表现为高程的快速下降，位置关系表现为邻接山体与平原。

② 冲洪积扇体在遥感影像上的几何形态通常呈现为近似扇形。其发育范围取决于山区河流的汇水面积、岩体风化特征、水文特征等。

③ 冲洪积扇体的物质组成取决于其河流上游汇水范围的岩体类型；不同的物质成分在遥感影像上呈现不同的光谱特征。

④ 冲洪积扇体的土壤颗粒较粗，主要为砂、砾，通常植被稀疏，不利于植物生长。扇体边缘和冲洪积扇前缘的土壤颗粒较细，一般为沙、粉沙及亚黏土，通常植被发育，植被覆盖率高。

⑤ 冲洪积扇在地形上呈现为向盆地倾斜的山前倾斜平原或倾斜台地，其坡度为5°～10°，扇体前缘的坡度为1°～2°。

2) 影像特征分析

冲洪积扇体的反射率随波长变化的趋势与冲洪积扇土体颗粒变化基本一致。由于扇体前缘土体颗粒变细，以沙、粉沙以及亚黏土为主，土体粒径约在0.45～0.25 μm之间，土壤孔隙率相对增高，土壤含水率也明显增高，其在各波段的反射率略低于冲洪积扇主体的反射率。冲洪积扇前缘地下水溢流带产生的“微阴影”以及相对高的土壤含水率，导致扇体前缘部分的反射率、植被指数和地表温度的突变，形成冲洪积扇前缘分界带。

3 冲洪积扇提取模型构建

3.1 特征指数构建

1) 植被特征指数

归一化差分植被指数[8](NDVI)常用于监测植被生长的状态、植被覆盖度，能够有效排除植被密集区，计算公式如下：

(1)

式(1)中:R表示影像的红波段，NIR表示影像的近红外波段。

2) 温度植被干旱指数

温度植被干旱指数[9](TVDI)可用于干旱检测以及研究干旱程度的空间变化特征，计算公式如下：

(2)

式中:Ts表示给定像元的地表温度，Tsmin、Tsmax表示给定像元相同NDVI值下的最低温度和最高温度。根据Landsat 8热红外波段影像反演研究区温度植被干旱指数，指数值越接近0，土壤湿度越大。反演结果如图2所示。

3) 纹理特征指数

灰度共生矩阵是一种描述纹理特征的常用方法，是建立在估计图像的二阶组合条件概率密度函数基础上的统计方法，能够有效提取纹理特征[10]。通过灰度共生矩阵一定程度上能够识别冲洪积扇顶的河道纹理。选用7×7窗口提取5种统计量，分别为一致性(Homogeneity)、对比度(Contrast)、相异性(Dissimilarity)、熵(Entropy)、均值(Mean)。其计算公式如下：

(3)

式(3)中:P(i,j)表示i行j列像元的亮度值，u为7×7窗口内的像元亮度均值。

4) 坡度、坡向

坡度(Slope)、坡向(Aspect)能充分呈现高程的变化率以及变化方向，有效识别冲洪积扇顶的地学特征。

(4)

式(4)中:m表示点在x方向上高程变化率，n表示点在y方向上高程变化率。坡度、坡向计算一般是在最小像元的基础上，以5×5窗口为单位进行局部微分或曲面回归拟合得出。

5) 光谱特征参数

冲洪积扇的物质组成主要以河流携带的沙砾石和粘土为主，分选性较好，植被极少。其在多光谱影像上色调比较均匀、亮度值变化差异较小。对冲洪积扇及该区域其他地物进行样本采集分析，结果如表2所示。

表2 冲洪积扇及其他地物光谱值

3.2 分类器设计

决策树是一种经典的分类器，其扩展通过节点的类别情况，每个类别导致一个结果[11]。分类器使用改进于传统决策树的CART(Classification And Regression Tree)分类方法，即分类与回归树[12]。

CART为在给定输入x的条件下，输出随机变量y的条件概率分布。CART算法二分每个特征(包括标签特征以及连续特征)，通过最优二分特征及其最优二分特征值的选择、切分，生成二叉树并进行剪枝[13]。分类过程中，回归CART树用误差平方和准则选择特征，分类CART树用基尼系数准则选择特征，并递归调用构建二叉树过程生成CART树。回归树实现流程如图3所示。

3.3 冲洪积扇提取模型构建

1) 计算DEM影像的坡度值。以3°和10°为阈值，将坡度低于3°的区域视为冲洪积扇缘发育区；将坡度为3°～10°的区域视为冲洪积扇体发育区；将坡度高于10°的区域视为冲洪积扇顶发育区。

2) 计算目标区域多光谱影像的纹理特征指数。冲洪积扇顶的影像纹理特征与山体存在明显区别，条带状河流冲刷痕迹明显。为缩减灰度共生矩阵计算量，基于Matlab软件平台对影像的灰度级数进行压缩，即从影像最小像元值起，每5级压缩为1级，如1～5视为1，6～10视为2。

3) 计算NDVI、TDVI指数。研究区地处高原，多盐碱地，植被指数主要与土壤类型相关，参考表2的平原台地地光谱特征指数，以0作为区分阈值。根据冲洪积扇体表面植被覆盖度极低的原则，将NDVI指数大于0的区域视为植被存在区域，排除在冲洪积扇可能存在的区域之外。从冲洪积扇根部体向扇体边缘的土壤含水率逐渐增高，到扇前缘的地下水水溢出带土壤含水率有明显跃升，多泉水发育。根据TDVI指数，以0.3作为阈值提取影像中的高表层土壤含水率的区域。

4) 根据影像光谱特征参数，以蓝波段像元值720为阈值，区分盐碱地与其他地物，排除在冲洪积扇可能存在的区域之外。

5) 基于决策树法排除有明显特征的非冲洪积扇地物，界定山体与山间盆地的分界线。基于CART树获取最优分割特征值，实现不同源头冲洪积扇的提取与区分。

6) 基于冲洪积扇与山体之间的邻接关系判定冲洪积扇，实现冲洪积扇提取。

4 结果验证与分析

为检验模型效果，分别使用本文提出的模型、最大似然法、SVM法对冲洪积扇进行提取，并通过混淆矩阵对提取结果进行精度评价，结果如图4、表3所示。

注：图中白色表示冲积扇A，灰色表示冲积扇B，黑色表示洪积扇C

表3 冲积扇分类结果

从表3看，本文提出的模型提取效果明显优于最大似然法和SVM法，对影像中3种不同源头的扇体提取结果的Kappa系数分别为0.95、0.86、0.90，总体精度达到92.37%。结合图4所示，模型通过冲洪积扇与山体之间的拓扑邻接关系，有效滤除了 “同谱异物”的其他背景地物，极大地提高了提取精度。

5 结论

1) 以Landsat8影像、GF-2卫星影像、数字高程模型为数据源，基于CART分类器的冲洪积扇提取模型，实现了冲洪积扇微地貌单元的精确提取，为区域地貌环境的监测提供技术支持。该模型提取冲积扇的总精度为92.37%。

2) DEM数据是冲洪积扇提取的重要参数，结合光学、地学特征的模型提取精度远优于仅基于光学特征的最大似然法和支持向量机方法。

3) 温度植被干旱指数能有效表达土壤湿度，对识别冲洪积扇顶部、前缘具有重要作用，能够准确划分冲洪积扇前缘地下水溢出带。