基于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.