基于DEM的地形因子分析与岩性分类

王 婷， 潘 军， 蒋立军， 邢立新， 于一凡， 王鹏举

(吉林大学地球探测科学与技术学院，长春 130026)

0 引言

岩性的识别与分类在遥感地质领域具有重要的研究意义。在遥感地质解译的目视识别工作中，除利用色调、几何形状等基本标志以外，更多通过水系和纹理图案等标志来进行岩性的解译[1-4]，而水系和纹理图案正是不同地形地貌特征在遥感影像上的平面表征，由此可知在岩性识别与分类中地形地貌特征信息是重要的参考数据之一。而地形地貌特征可以借助由数字高程模型(digital elevation model，DEM)数据得到的高程因子及派生因子来进行数字表达[5-6]，因此将地形因子参数应用到遥感岩性识别与分类中对岩性识别精度的提高具有重要的意义。

近年来，已有研究将地形因子参数用于岩性的分类，分类精度得到了显著提高[7-10]。为了明确分析地形因子对识别每一类岩性所起的具体作用，本文将地形因子对岩性分类的有效性做了总结，并结合地形因子对岩性分类的结果图进行分析，以期得到对不同类别岩性识别的最佳地形因子组合。

1 研究区概况及数据源

研究区位于我国黑龙江西北部，大兴安岭南部松岭区(图1)。

图1 松岭区遥感影像

该区主要由伊勒呼里山绵延南伸的2条低山丘陵组成，西北高、东南低，海拔高度为400～800 m。地质调查表明，区内出露11类岩性单元，主要以花岗岩为主，如图2所示(图中岩性单元代号0为未识别岩性)。研究采用的地形数据为30 m空间分辨率的ASTER GDEM数据，能很好地反映地物的形态特征。研究区范围为N51°00′～51°20′，E124°00′～124°30′，空间范围涉及4幅DEM影像，在ENVI软件中将其进行镶嵌和裁剪等制图处理，结果如图3所示。

0-未识别岩性； 1-吉祥峰组： 深灰、灰绿色流纹岩和流纹质凝灰岩； 2-正长花岗岩； 3-二长花岗岩； 4- 花岗闪长岩； 5-上库力组： 流纹岩、流纹质凝灰岩、流纹质凝灰熔岩； 6-大网子组： 变英安质熔结凝灰岩、变英安岩、变砂岩； 7-小古里河组： 强片理化流纹岩、流纹质熔岩和凝灰岩； 8-低河漫滩堆积层： 黄褐色砂砾石、砂、粘土和淤泥； 9-碱长花岗岩； 10-石英闪长岩； 11-红水泉组： 变质砾岩、变质砂岩、含早石炭世昆虫

图2松岭区岩性单元分布

Fig.2LithologicunitdistributioninSongling

图3 松岭区 DEM影像

2 地形因子简介

在DEM数据中每个像元的DN值为高程。通过中心像元与邻近像元DN值的运算或一定区域内所有像元DN值的统计，即可得到反映地表地貌特征的各个地形因子参数。本文采用的地形因子参数包括高程、坡度、剖面曲率、平面曲率、纵向曲率、横向曲率、地形起伏度、地表切割深度、地表粗糙度和高程变异系数。其中： ①高程是地面点沿铅垂线到大地水准面的距离； ②坡度反映曲面的倾斜程度，用垂直高差和水平距离的比值表示； ③剖面曲率是对地表曲面在垂直方向上高程变化率的度量； ④平面曲率是对地表曲面在水平方向上扭曲变化的度量； ⑤纵向曲率是沿下坡方向的坡度变化率； ⑥横向曲率是沿下坡方向的垂直方向上的坡度变化率； ⑦地形起伏度是一个特定分析区域内高程最大值和最小值之差； ⑧地表切割深度是一个特定分析区域内高程均值与最小值之差； ⑨地表粗糙度是地形的曲面面积与其在水平面上的投影面积之比，即坡度的余弦值的倒数； ⑩高程变异系数是区域内高程标准差与均值之比[11-13]。在实际研究中，经常用因子①—⑥各自的均值来反映非点状要素信息。

3 地形因子的分类有效性评价

在DEM影像上以专家遥感目视解译成果图(图2)为基准划分岩性单元，将各个岩性单元(这一部分不考虑低河漫滩堆积层)的地形因子数据提取或计算出来时，由于地形因子间往往存在着重叠信息并且有些因子区分岩性的效果有一定的相似性，本文首先通过分析地形因子间的相关系数，以及各因子区分岩性的相似性和差异性来对地形因子进行筛选，剔除相关性高、区分岩性效果相似的地形因子，然后深入分析筛选出的地形因子适宜识别区分的岩性。

每一岩性区的地形因子提取分为高程因子提取、微观地形因子提取和宏观地形因子提取3个方面，其中微观地形因子包括坡度、剖面曲率、平面曲率、纵向曲率和横向曲率； 宏观因子包括地形起伏度、地表切割深度、地表粗糙度和高程变异系数。高程因子通过DEM数据直接获取，微观地形因子利用ENVI软件的地形工具获取，再统计高程因子和微观地形因子的均值数据。对于宏观地形因子，地形起伏度由每一岩性区高程的最大值减最小值得到； 地表切割深度由每一岩性区高程的均值减最小值得到； 高程变异系数由每一岩性区高程的标准差与均值的比值得到； 每一岩性区的地表粗糙度表示为

(1)

式中：R为地表粗糙度；slopei为每一像元的坡度；n为分析区域内的总像元数。

表1示出各因子的区分情况，每一地形因子对应的区分情况和分开数量是以上述获取到的每一岩性区的所有地形因子数据为基础得到的。具体方法为： 首先，根据每一岩性区的10种地形因子数据统计每一类岩性单元的所有地形因子数据范围，再根据范围大小来区分岩性，从而得到每一地形因子对应的10类岩性两两分开的区分情况和分开数量(表1中区分情况部分的每个数字代表每一类岩性，与图2中的岩性编号一致)； 然后，依据10个地形因子之间的相关系数得到纵向曲率与剖面曲率的相关性强，以及坡度、高程变异系数、地形起伏度和地表切割深度4个因子之间的相关性强，并对相关性强的地形因子结合地形因子区分岩性的相似性和差异性来对其进行取舍。

表1 各因子的区分情况Tab.1 Distinction between the factors

①“区分情况”中的数字对应于图2的岩性编号。

从表1中可以看出纵向曲率与剖面曲率对于岩性的区分情况相似，并且纵向曲率能区分开的岩性剖面曲率也几乎都能分开，因此从两者中剔除纵向曲率而选择剖面曲率; 在地表切割深度、坡度、高程变异系数和地形起伏度中，地表切割深度的分开数量最佳，且与其他因子相比能将岩性5,6和7以及岩性1和6分开; 地形起伏度和高程变异系数所能区分开的岩性由高程、剖面曲率和地表切割深度也几乎都能分开，因此选择地表切割深度; 坡度区分岩性的效果与地表粗糙度相似，且能区分开的岩性种类少于地表粗糙度; 此外横向曲率能分开的岩性，平面曲率也能分开，因此剔除横向曲率。最后筛选出用于岩性分类的地形因子为高程、剖面曲率、地表切割深度、地表粗糙度和平面曲率5种。

对于以上筛选出的5个因子根据图4分析每个地形因子适宜识别区分的岩性。

(a) 高程 (b) 剖面曲率 (c) 平面曲率(d) 地表切割深度(e) 地表粗糙度

图4不同岩性的地形因子范围结果图

(注： 图中岩性标号对应于图2的岩性编号)

Fig.4Patternsoftopographicvariablesrangeofdifferentlithologies

由图可知，高程因子能分开的岩性类别最多，由图4(a)还可以直观看出其对岩性1,2,6,9,11的区分性强，岩性1高程值最大； 岩性11几乎与除了岩性1以外的岩性完全区分开； 岩性2可与岩性6,7,9,10,11区分开； 岩性6与岩性2,7,9,11区分开，与岩性10也只有一点重叠； 岩性9与岩性2,6,11区分开。由图4(b)可以直观看出剖面曲率因子对岩性1,3,10,11的区分性强，相较于其他岩性，岩性1的剖面曲率值最大； 岩性3与岩性5,7,9区分开，且与岩性1和10只有一点重叠； 岩性10可与岩性1,5,6,7,9区分开，与岩性2,3,4也都只有一点重叠； 岩性11可与岩性1,5,7,9区分开。对于平面曲率因子，由图4(c)可以直观看出其对岩性1,10,11的区分性强，岩性1平面曲率值是最小的，岩性1与岩性10,11明显区分开； 岩性10与岩性1,5，9分开，与岩性2,7也只有一点重叠； 岩性11与岩性1,9分开，与岩性5,7也只有一点重叠。对于地表切割深度因子，由图4(d)可以直观看出其对岩性5,6,10的区分性强，岩性5几乎可以与岩性4,6,7,9,10,11区分开； 岩性6与岩性1,5完全分开，与岩性4,7,9也几乎可以分开； 岩性10几乎与岩性1,2,3,4,5,7,9都区分开，另外相比较其他因子其对岩性4的区分力是最好的，岩性4与岩性5,6,10,11都只有一点重叠。对于地表粗糙度因子，由图4(e)可以直观看出其对岩性7的区分性最强，另外岩性5的粗糙度值是所有岩性中最大的，岩性7除了与岩性5有点重叠以外与其余岩性全部区分开； 岩性5与岩性2,3,4,6,11也都只有一点重叠。

基于以上分析可以看出5个地形因子对于岩性分类效果显著，且每类岩性都有其特有的地形因子组合可以将其与其他类岩性区分开。

3.1 地形因子基础数据提取

研究区用于岩性分类的5个地形因子数据都是一定大小区域内的统计数据，可以通过ENVI软件以模板窗口的形式运算得到。由于不同窗口下的地形因子数据不同，因此需要考虑窗口问题，最佳窗口下得到的地形因子能准确反映地形地貌特征，将其用于岩性分类能得到最好的岩性分类效果(表2)。

表2 地形因子的单元大小与均值的对应情况Tab．2 Correspondence between the unit size of the terrain factor and the mean value

从表2可以看出随窗口由小到大变化，地表粗糙度均值、高程均值、平面曲率均值和剖面曲率均值的变化幅度都很小，说明窗口大小对其影响不大，因而选取3像素×3像素窗口下的这4种地形因子数据作为岩性分类的基础数据； 而地表切割深度随窗口的增大，其均值变化幅度较大，所以增加提取地表切割深度的窗口数量进一步研究其变化趋势(表3)。

表3 地形因子的单元大小与地表切割深度均值的对应情况Tab.3 Correspondence between the unit size of the terrainfactor and the mean value of the surface cutting depth

可以看出地表切割深度均值随窗口的增大在不断地增大，且增大幅度在逐渐减小。

将其拟合为对数方程，即

y=15.454lnx-74.809，

(2)

式中：x为单元面积；y为地表切割深度均值。

对应的拟合曲线如图5所示，其确定性系数R2=0.957 7，拟合程度较好，并通过了统计学检验。在图5中可以看到存在增大幅度由快变慢的拐点，这一点所对应的窗口大小即为最佳窗口。

图5 地表切割深度均值与单元面积对应关系拟合曲线

综合以上分析，可采用均值变点分析法来确定最佳窗口[14-15]。其具体步骤如下：

1)计算不同窗口大小下单位面积上的地表切割深度均值，并对其求对数得到要统计的数列{xK},K=1，…，N,N为样本数。

3)计算统计量SK和S，其中SK为2段样本的离差平方和之和,S为总的离差平方和。

变点的存在会使整个样本的变量S与样本分段后的变量SK之间的差距增大，最大差值对应的窗口大小即为最佳窗口单元。表4示出均值变点法的统计结果，可以看出当K=9时差值最大，对应的窗口大小为19，因此选用19像元×19像元窗口下的地表切割深度数据作为岩性分类的基础数据。

表4 均值变点法的统计结果Tab.4 Statistical results of mean change point method

3.2 分类研究与对比分析

为了有效地运用地形因子进行岩性分类，采用ISODATA迭代自组织的数据分析方法对研究区进行非监督分类，共分出8类岩石，结果如图6所示。

1-大网子组： 变英安质熔结凝灰岩、变英安岩、变砂岩； 2-花岗闪长岩； 3-正长花岗岩； 4-吉祥峰组： 深灰、灰绿色流纹岩和流纹质凝灰岩； 5-二长花岗岩； 6-上库力组： 流纹岩、流纹质凝灰岩、流纹质凝灰熔岩； 7-碱长花岗岩； 8-低河漫滩堆积层： 黄褐色砂砾石、砂、粘土和淤泥

图6本文方法图像分类结果

Fig.6Imageclassificationresult

专家遥感目视解译成果图(图2)共有11类岩石(包括低河漫滩堆积层)，为了精确了解分类后的8类岩石与岩性单元分布图上11类岩石的对应情况，在MapGIS中将2幅图进行区域判别分析，得到分类后的8类岩石与已知11类岩石面积对应的属性信息(表5)，根据对应面积对8类岩石定义类别如表6。

表5 岩性单元面积对应情况Tab．5 Corresponds result of the lithologic unit area (m2)

表6 岩性单元对应情况Tab．6 Corresponds result of lithologic unit

3.3 精度评价与分析

由表5和表6可知每类岩性都有误分和漏分，但几乎都被识别出来了，只有小古里河组、石英闪长岩和红水泉组没有分出来，原因是这3类岩性区域太小，在分类参数设置时更侧重于使整体分类效果达到最佳，导致有的小区域没有分出来，此外3类岩性都只分布于一小块区域内，可能代表性较差，没能充分表达出其特有的地形地貌特征。区内的第四纪，即低河漫滩堆积层，其高程明显偏低，因此利用高程因子就可以将其分出来，其误分情况严重的主要原因是岩性单元分布图上多以细条的树枝状呈现而被处理到了临近的岩性类中。吉祥峰组岩性的平面曲率偏小而剖面曲率偏大，因而最适宜将其与其他岩性区分开； 此外其高程值在所有岩性中最大。由于一部分上库力组和正长花岗岩岩性的高程值也偏大，因此对于区分吉祥峰组与除了上述2种以外的岩性还可以通过高程因子实现。正长花岗岩的高程值较大，与大部分岩性区分明显。上库力组的地表切割深度与地表粗糙度较大，是将其与其他岩性区分开的最佳因子。从分类结果中可以看出吉祥峰组、上库力组与正长花岗岩3种岩性之间有一定的混淆，主要原因是三者的上述因子值都有一定程度的重叠。二长花岗岩的剖面曲率偏小，因此将其与其他岩性区分开的效果最佳。对于大网子组，其高程和地表切割深度值都偏小，因此这2个因子为将其区分出来的最佳因子。从分类结果中还可以看出二长花岗岩与大网子组岩性混淆严重主要是因二者5个地形因子值都有一定程度的重叠。对于花岗闪长岩，地表切割深度因子是将其与其他岩性区分开的最佳因子，但是可以看出分类图上花岗闪长岩分类精度很低，主要是其在5个地形因子值上与其他岩性都有一定的重叠。对于碱长花岗岩，地表粗糙度因子将其区分出来的效果最佳，但是其分类的误分和漏分情况都比较严重，原因可能是其在研究区内岩性区范围较小。

4 结论

1)以研究区30 m空间分辨率的DEM数据为基础提取地形因子，通过比较分析地形因子区分岩性的有效性，得出地形因子中高程、剖面曲率、平面曲率、地表切割深度和地表粗糙度5个因子都有其最适宜识别区分的岩性，将5种因子用于岩性分类，得到的分类结果中主要岩性区均能被识别出来。

2)研究结果表明岩性与地形因子间存在一定的相关性，主要的岩性类都有其特有的地形特征都可以通过一定的地形因子组合而被识别出来。因此加入地形因子参数会明显提高岩性的识别能力，这也为有植被覆盖的岩性区解译精度的提高提供了一种可能。为了深入研究各类岩性的地形特征从而进一步提高岩性的解译精度，还需挖掘岩性的更多地形信息，因此下一步需要研究更多的地形因子。

参考文献(References)：

[1] 余海阔,李培军.运用LANDSAT ETM+和ASTER数据进行岩性分类[J].岩石学报,2010,26(1):345-351.

Yu H K,Li P J.Lithologic mapping using LANDSAT ETM+ and ASTER data[J].Acta Petrologica Sinica,2010,26(1):345-351.

[2] 王晓东.水系提取方法研究及其地质意义[D].长春:吉林大学,2015.

Wang X D.The Research of Drainage Extraction Method and Its Geological Significance[D].Changchun:Jinlin University,2015.

[3] 于亚凤,杨金中,陈圣波,等.基于光谱指数的遥感影像岩性分类[J].地球科学,2015,40(8):1415-1419.

Yu Y F,Yang J Z,Chen S B,et al.Lithologic classification from remote sensing images based on spectral index[J].Earth Science,2015,40(8):1415-1419.

[4] 黄颖端,李培军,李争晓.基于地统计学的图像纹理在岩性分类中的应用[J].国土资源遥感,2003,15(3):45-49.doi:10.6046/gtzyyg.2003.03.11

Huang Y D,Li P J,Li Z X.The application of geostatistical image texture to remote sensing lithological classification[J].Remote Sensing for Land and Resources,2003,15(3):45-49.doi:10.6046/gtzyyg.2003.03.11

[5] 曾德耀.基于最佳地形因子组合的地貌形态类型划分研究[D].重庆:重庆交通大学,2015.

Zeng D Y.Classification of Relief Form Based on the Best Terrain Factor Combination[D].Chongqing:Chongqing Jiaotong University,2015.

[6] 杨晏立,何政伟,杨 斌,等.最佳因子复合的四川省地貌类型自动划分[J].陕西理工学院学报(自然科学版)2009,25(4):74-79.

Yang Y L,He Z W,Yang B,et al.Automatic classification of landform types in Sichuan Province with the optimum factors complex[J].Journal of Shaanxi University of Technology(Natural Science Edition)2009,25(4):74-79.

[7] 姜莎莎,李培军.基于ASTER图像和地形因子的岩性单元分类——以新疆木垒地区为例[J].地球信息科学学报,2011,13(6):825-832.

Jiang S S,Li P J.Lithologic unit mapping using ASTER data and topographic variables:A case study of Mulei area of XinJiang[J].Journal of Geo-Information Science,2011,13(6):825-832.

[8] Grebby S,Cunningham D,Naden J,et al.Lithological mapping of the Troodos ophiolite,Cyprus,using airborne LiDAR topographic data[J].Remote Sensing of Environment,2010,114(4):713-724.

[9] Grebby S,Naden J,Cunningham D,et al.Integrating airborne multispectral imagery and airborne LiDAR data for enhanced lithological mapping in vegetated terrain[J].Remote Sensing of Environment,2011,115(1):214-226.

[10] Li P J,Cheng T,Guo J C.Multivariate image texture by multivariate variogram for multispectral image classification[J].Photogrammetric Engineering and Remote Sensing,2009,75(2):147-157.

[11] 周启明,刘学军.数字地形分析[M].北京:科学出版社,2006:52-75.

Zhou Q M,Liu X J.Digital Terrain Analysis[M].Beijing: Science Press,2006:52-75.

[12] 刘少峰,王 陶,张会平,等.数字高程模型在地表过程研究中的应用[J].地学前缘,2005,12(1):303-309.

Liu S F,Wang T,Zhang H P,et al.Application of digital elevation model to surficial process research[J].Earth Science Frontiers,2005,12(1):303-309.

[13] 杨 昕,汤国安,刘学军,等.数字地形分析的理论、方法与应用[J].地理学报,2009,64(9):1058-1070.

Yang X,Tang G A,Liu X J,et al.Digital terrain analysis:Theory, method and application[J].Acta Geographica Sinica,2009,64(9):1058-1070.

[14] 赵斌滨,程永锋,丁士君,等.基于SRTM-DEM的我国地势起伏度统计单元研究[J].水利学报,2015,46(s1):284-290.

Zhao B B,Cheng Y F,Ding S J,et al.Statistical unit of relief amplitude in China based on SRTM-DEM[J].Journal of Hydraulic Engineering,2015,46(s1):284-290.

[15] 高 蜻,唐丽霞,谷晓平,等.基于ArcGIS的望谟河流域地势起伏度分析[J].中国水土保持科学,2015,13(4):9-14.

Gao Q,Tang L X,Gu X P,et al.Analysis of ArcGIS-based relief amplitude of the Wangmo River watershed in Guizhou[J].Science of Soil and Water Conservation,2015,13(4):9-14.