王 婷， 潘 軍， 蔣立軍， 邢立新， 于一凡， 王鵬舉

(吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026)

0 引言

巖性的識別與分類在遙感地質(zhì)領(lǐng)域具有重要的研究意義。在遙感地質(zhì)解譯的目視識別工作中，除利用色調(diào)、幾何形狀等基本標(biāo)志以外，更多通過水系和紋理圖案等標(biāo)志來進(jìn)行巖性的解譯[1-4]，而水系和紋理圖案正是不同地形地貌特征在遙感影像上的平面表征，由此可知在巖性識別與分類中地形地貌特征信息是重要的參考數(shù)據(jù)之一。而地形地貌特征可以借助由數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)得到的高程因子及派生因子來進(jìn)行數(shù)字表達(dá)[5-6]，因此將地形因子參數(shù)應(yīng)用到遙感巖性識別與分類中對巖性識別精度的提高具有重要的意義。

近年來，已有研究將地形因子參數(shù)用于巖性的分類，分類精度得到了顯著提高[7-10]。為了明確分析地形因子對識別每一類巖性所起的具體作用，本文將地形因子對巖性分類的有效性做了總結(jié)，并結(jié)合地形因子對巖性分類的結(jié)果圖進(jìn)行分析，以期得到對不同類別巖性識別的最佳地形因子組合。

1 研究區(qū)概況及數(shù)據(jù)源

研究區(qū)位于我國黑龍江西北部，大興安嶺南部松嶺區(qū)(圖1)。

圖1 松嶺區(qū)遙感影像

該區(qū)主要由伊勒呼里山綿延南伸的2條低山丘陵組成，西北高、東南低，海拔高度為400～800 m。地質(zhì)調(diào)查表明，區(qū)內(nèi)出露11類巖性單元，主要以花崗巖為主，如圖2所示(圖中巖性單元代號0為未識別巖性)。研究采用的地形數(shù)據(jù)為30 m空間分辨率的ASTER GDEM數(shù)據(jù)，能很好地反映地物的形態(tài)特征。研究區(qū)范圍為N51°00′～51°20′，E124°00′～124°30′，空間范圍涉及4幅DEM影像，在ENVI軟件中將其進(jìn)行鑲嵌和裁剪等制圖處理，結(jié)果如圖3所示。

0-未識別巖性； 1-吉祥峰組： 深灰、灰綠色流紋巖和流紋質(zhì)凝灰?guī)r； 2-正長花崗巖； 3-二長花崗巖； 4- 花崗閃長巖； 5-上庫力組： 流紋巖、流紋質(zhì)凝灰?guī)r、流紋質(zhì)凝灰熔巖； 6-大網(wǎng)子組： 變英安質(zhì)熔結(jié)凝灰?guī)r、變英安巖、變砂巖； 7-小古里河組： 強片理化流紋巖、流紋質(zhì)熔巖和凝灰?guī)r； 8-低河漫灘堆積層： 黃褐色砂礫石、砂、粘土和淤泥； 9-堿長花崗巖； 10-石英閃長巖； 11-紅水泉組： 變質(zhì)礫巖、變質(zhì)砂巖、含早石炭世昆蟲

圖2松嶺區(qū)巖性單元分布

Fig.2LithologicunitdistributioninSongling

圖3 松嶺區(qū) DEM影像

2 地形因子簡介

在DEM數(shù)據(jù)中每個像元的DN值為高程。通過中心像元與鄰近像元DN值的運算或一定區(qū)域內(nèi)所有像元DN值的統(tǒng)計，即可得到反映地表地貌特征的各個地形因子參數(shù)。本文采用的地形因子參數(shù)包括高程、坡度、剖面曲率、平面曲率、縱向曲率、橫向曲率、地形起伏度、地表切割深度、地表粗糙度和高程變異系數(shù)。其中： ①高程是地面點沿鉛垂線到大地水準(zhǔn)面的距離； ②坡度反映曲面的傾斜程度，用垂直高差和水平距離的比值表示； ③剖面曲率是對地表曲面在垂直方向上高程變化率的度量； ④平面曲率是對地表曲面在水平方向上扭曲變化的度量； ⑤縱向曲率是沿下坡方向的坡度變化率； ⑥橫向曲率是沿下坡方向的垂直方向上的坡度變化率； ⑦地形起伏度是一個特定分析區(qū)域內(nèi)高程最大值和最小值之差； ⑧地表切割深度是一個特定分析區(qū)域內(nèi)高程均值與最小值之差； ⑨地表粗糙度是地形的曲面面積與其在水平面上的投影面積之比，即坡度的余弦值的倒數(shù)； ⑩高程變異系數(shù)是區(qū)域內(nèi)高程標(biāo)準(zhǔn)差與均值之比[11-13]。在實際研究中，經(jīng)常用因子①—⑥各自的均值來反映非點狀要素信息。

3 地形因子的分類有效性評價

在DEM影像上以專家遙感目視解譯成果圖(圖2)為基準(zhǔn)劃分巖性單元，將各個巖性單元(這一部分不考慮低河漫灘堆積層)的地形因子數(shù)據(jù)提取或計算出來時，由于地形因子間往往存在著重疊信息并且有些因子區(qū)分巖性的效果有一定的相似性，本文首先通過分析地形因子間的相關(guān)系數(shù)，以及各因子區(qū)分巖性的相似性和差異性來對地形因子進(jìn)行篩選，剔除相關(guān)性高、區(qū)分巖性效果相似的地形因子，然后深入分析篩選出的地形因子適宜識別區(qū)分的巖性。

每一巖性區(qū)的地形因子提取分為高程因子提取、微觀地形因子提取和宏觀地形因子提取3個方面，其中微觀地形因子包括坡度、剖面曲率、平面曲率、縱向曲率和橫向曲率； 宏觀因子包括地形起伏度、地表切割深度、地表粗糙度和高程變異系數(shù)。高程因子通過DEM數(shù)據(jù)直接獲取，微觀地形因子利用ENVI軟件的地形工具獲取，再統(tǒng)計高程因子和微觀地形因子的均值數(shù)據(jù)。對于宏觀地形因子，地形起伏度由每一巖性區(qū)高程的最大值減最小值得到； 地表切割深度由每一巖性區(qū)高程的均值減最小值得到； 高程變異系數(shù)由每一巖性區(qū)高程的標(biāo)準(zhǔn)差與均值的比值得到； 每一巖性區(qū)的地表粗糙度表示為

(1)

式中：R為地表粗糙度；slopei為每一像元的坡度；n為分析區(qū)域內(nèi)的總像元數(shù)。

表1示出各因子的區(qū)分情況，每一地形因子對應(yīng)的區(qū)分情況和分開數(shù)量是以上述獲取到的每一巖性區(qū)的所有地形因子數(shù)據(jù)為基礎(chǔ)得到的。具體方法為： 首先，根據(jù)每一巖性區(qū)的10種地形因子數(shù)據(jù)統(tǒng)計每一類巖性單元的所有地形因子數(shù)據(jù)范圍，再根據(jù)范圍大小來區(qū)分巖性，從而得到每一地形因子對應(yīng)的10類巖性兩兩分開的區(qū)分情況和分開數(shù)量(表1中區(qū)分情況部分的每個數(shù)字代表每一類巖性，與圖2中的巖性編號一致)； 然后，依據(jù)10個地形因子之間的相關(guān)系數(shù)得到縱向曲率與剖面曲率的相關(guān)性強，以及坡度、高程變異系數(shù)、地形起伏度和地表切割深度4個因子之間的相關(guān)性強，并對相關(guān)性強的地形因子結(jié)合地形因子區(qū)分巖性的相似性和差異性來對其進(jìn)行取舍。

表1 各因子的區(qū)分情況Tab.1 Distinction between the factors

①“區(qū)分情況”中的數(shù)字對應(yīng)于圖2的巖性編號。

從表1中可以看出縱向曲率與剖面曲率對于巖性的區(qū)分情況相似，并且縱向曲率能區(qū)分開的巖性剖面曲率也幾乎都能分開，因此從兩者中剔除縱向曲率而選擇剖面曲率; 在地表切割深度、坡度、高程變異系數(shù)和地形起伏度中，地表切割深度的分開數(shù)量最佳，且與其他因子相比能將巖性5,6和7以及巖性1和6分開; 地形起伏度和高程變異系數(shù)所能區(qū)分開的巖性由高程、剖面曲率和地表切割深度也幾乎都能分開，因此選擇地表切割深度; 坡度區(qū)分巖性的效果與地表粗糙度相似，且能區(qū)分開的巖性種類少于地表粗糙度; 此外橫向曲率能分開的巖性，平面曲率也能分開，因此剔除橫向曲率。最后篩選出用于巖性分類的地形因子為高程、剖面曲率、地表切割深度、地表粗糙度和平面曲率5種。

對于以上篩選出的5個因子根據(jù)圖4分析每個地形因子適宜識別區(qū)分的巖性。

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

圖4不同巖性的地形因子范圍結(jié)果圖

(注： 圖中巖性標(biāo)號對應(yīng)于圖2的巖性編號)

Fig.4Patternsoftopographicvariablesrangeofdifferentlithologies

由圖可知，高程因子能分開的巖性類別最多，由圖4(a)還可以直觀看出其對巖性1,2,6,9,11的區(qū)分性強，巖性1高程值最大； 巖性11幾乎與除了巖性1以外的巖性完全區(qū)分開； 巖性2可與巖性6,7,9,10,11區(qū)分開； 巖性6與巖性2,7,9,11區(qū)分開，與巖性10也只有一點重疊； 巖性9與巖性2,6,11區(qū)分開。由圖4(b)可以直觀看出剖面曲率因子對巖性1,3,10,11的區(qū)分性強，相較于其他巖性，巖性1的剖面曲率值最大； 巖性3與巖性5,7,9區(qū)分開，且與巖性1和10只有一點重疊； 巖性10可與巖性1,5,6,7,9區(qū)分開，與巖性2,3,4也都只有一點重疊； 巖性11可與巖性1,5,7,9區(qū)分開。對于平面曲率因子，由圖4(c)可以直觀看出其對巖性1,10,11的區(qū)分性強，巖性1平面曲率值是最小的，巖性1與巖性10,11明顯區(qū)分開； 巖性10與巖性1,5，9分開，與巖性2,7也只有一點重疊； 巖性11與巖性1,9分開，與巖性5,7也只有一點重疊。對于地表切割深度因子，由圖4(d)可以直觀看出其對巖性5,6,10的區(qū)分性強，巖性5幾乎可以與巖性4,6,7,9,10,11區(qū)分開； 巖性6與巖性1,5完全分開，與巖性4,7,9也幾乎可以分開； 巖性10幾乎與巖性1,2,3,4,5,7,9都區(qū)分開，另外相比較其他因子其對巖性4的區(qū)分力是最好的，巖性4與巖性5,6,10,11都只有一點重疊。對于地表粗糙度因子，由圖4(e)可以直觀看出其對巖性7的區(qū)分性最強，另外巖性5的粗糙度值是所有巖性中最大的，巖性7除了與巖性5有點重疊以外與其余巖性全部區(qū)分開； 巖性5與巖性2,3,4,6,11也都只有一點重疊。

基于以上分析可以看出5個地形因子對于巖性分類效果顯著，且每類巖性都有其特有的地形因子組合可以將其與其他類巖性區(qū)分開。

3.1 地形因子基礎(chǔ)數(shù)據(jù)提取

研究區(qū)用于巖性分類的5個地形因子數(shù)據(jù)都是一定大小區(qū)域內(nèi)的統(tǒng)計數(shù)據(jù)，可以通過ENVI軟件以模板窗口的形式運算得到。由于不同窗口下的地形因子數(shù)據(jù)不同，因此需要考慮窗口問題，最佳窗口下得到的地形因子能準(zhǔn)確反映地形地貌特征，將其用于巖性分類能得到最好的巖性分類效果(表2)。

表2 地形因子的單元大小與均值的對應(yīng)情況Tab．2 Correspondence between the unit size of the terrain factor and the mean value

從表2可以看出隨窗口由小到大變化，地表粗糙度均值、高程均值、平面曲率均值和剖面曲率均值的變化幅度都很小，說明窗口大小對其影響不大，因而選取3像素×3像素窗口下的這4種地形因子數(shù)據(jù)作為巖性分類的基礎(chǔ)數(shù)據(jù)； 而地表切割深度隨窗口的增大，其均值變化幅度較大，所以增加提取地表切割深度的窗口數(shù)量進(jìn)一步研究其變化趨勢(表3)。

表3 地形因子的單元大小與地表切割深度均值的對應(yīng)情況Tab.3 Correspondence between the unit size of the terrainfactor and the mean value of the surface cutting depth

可以看出地表切割深度均值隨窗口的增大在不斷地增大，且增大幅度在逐漸減小。

將其擬合為對數(shù)方程，即

y=15.454lnx-74.809，

(2)

式中：x為單元面積；y為地表切割深度均值。

對應(yīng)的擬合曲線如圖5所示，其確定性系數(shù)R2=0.957 7，擬合程度較好，并通過了統(tǒng)計學(xué)檢驗。在圖5中可以看到存在增大幅度由快變慢的拐點，這一點所對應(yīng)的窗口大小即為最佳窗口。

圖5 地表切割深度均值與單元面積對應(yīng)關(guān)系擬合曲線

綜合以上分析，可采用均值變點分析法來確定最佳窗口[14-15]。其具體步驟如下：

1)計算不同窗口大小下單位面積上的地表切割深度均值，并對其求對數(shù)得到要統(tǒng)計的數(shù)列{xK},K=1，…，N,N為樣本數(shù)。

3)計算統(tǒng)計量SK和S，其中SK為2段樣本的離差平方和之和,S為總的離差平方和。

變點的存在會使整個樣本的變量S與樣本分段后的變量SK之間的差距增大，最大差值對應(yīng)的窗口大小即為最佳窗口單元。表4示出均值變點法的統(tǒng)計結(jié)果，可以看出當(dāng)K=9時差值最大，對應(yīng)的窗口大小為19，因此選用19像元×19像元窗口下的地表切割深度數(shù)據(jù)作為巖性分類的基礎(chǔ)數(shù)據(jù)。

表4 均值變點法的統(tǒng)計結(jié)果Tab.4 Statistical results of mean change point method

3.2 分類研究與對比分析

為了有效地運用地形因子進(jìn)行巖性分類，采用ISODATA迭代自組織的數(shù)據(jù)分析方法對研究區(qū)進(jìn)行非監(jiān)督分類，共分出8類巖石，結(jié)果如圖6所示。

1-大網(wǎng)子組： 變英安質(zhì)熔結(jié)凝灰?guī)r、變英安巖、變砂巖； 2-花崗閃長巖； 3-正長花崗巖； 4-吉祥峰組： 深灰、灰綠色流紋巖和流紋質(zhì)凝灰?guī)r； 5-二長花崗巖； 6-上庫力組： 流紋巖、流紋質(zhì)凝灰?guī)r、流紋質(zhì)凝灰熔巖； 7-堿長花崗巖； 8-低河漫灘堆積層： 黃褐色砂礫石、砂、粘土和淤泥

圖6本文方法圖像分類結(jié)果

Fig.6Imageclassificationresult

專家遙感目視解譯成果圖(圖2)共有11類巖石(包括低河漫灘堆積層)，為了精確了解分類后的8類巖石與巖性單元分布圖上11類巖石的對應(yīng)情況，在MapGIS中將2幅圖進(jìn)行區(qū)域判別分析，得到分類后的8類巖石與已知11類巖石面積對應(yīng)的屬性信息(表5)，根據(jù)對應(yīng)面積對8類巖石定義類別如表6。

表5 巖性單元面積對應(yīng)情況Tab．5 Corresponds result of the lithologic unit area (m2)

表6 巖性單元對應(yīng)情況Tab．6 Corresponds result of lithologic unit

3.3 精度評價與分析

由表5和表6可知每類巖性都有誤分和漏分，但幾乎都被識別出來了，只有小古里河組、石英閃長巖和紅水泉組沒有分出來，原因是這3類巖性區(qū)域太小，在分類參數(shù)設(shè)置時更側(cè)重于使整體分類效果達(dá)到最佳，導(dǎo)致有的小區(qū)域沒有分出來，此外3類巖性都只分布于一小塊區(qū)域內(nèi)，可能代表性較差，沒能充分表達(dá)出其特有的地形地貌特征。區(qū)內(nèi)的第四紀(jì)，即低河漫灘堆積層，其高程明顯偏低，因此利用高程因子就可以將其分出來，其誤分情況嚴(yán)重的主要原因是巖性單元分布圖上多以細(xì)條的樹枝狀呈現(xiàn)而被處理到了臨近的巖性類中。吉祥峰組巖性的平面曲率偏小而剖面曲率偏大，因而最適宜將其與其他巖性區(qū)分開； 此外其高程值在所有巖性中最大。由于一部分上庫力組和正長花崗巖巖性的高程值也偏大，因此對于區(qū)分吉祥峰組與除了上述2種以外的巖性還可以通過高程因子實現(xiàn)。正長花崗巖的高程值較大，與大部分巖性區(qū)分明顯。上庫力組的地表切割深度與地表粗糙度較大，是將其與其他巖性區(qū)分開的最佳因子。從分類結(jié)果中可以看出吉祥峰組、上庫力組與正長花崗巖3種巖性之間有一定的混淆，主要原因是三者的上述因子值都有一定程度的重疊。二長花崗巖的剖面曲率偏小，因此將其與其他巖性區(qū)分開的效果最佳。對于大網(wǎng)子組，其高程和地表切割深度值都偏小，因此這2個因子為將其區(qū)分出來的最佳因子。從分類結(jié)果中還可以看出二長花崗巖與大網(wǎng)子組巖性混淆嚴(yán)重主要是因二者5個地形因子值都有一定程度的重疊。對于花崗閃長巖，地表切割深度因子是將其與其他巖性區(qū)分開的最佳因子，但是可以看出分類圖上花崗閃長巖分類精度很低，主要是其在5個地形因子值上與其他巖性都有一定的重疊。對于堿長花崗巖，地表粗糙度因子將其區(qū)分出來的效果最佳，但是其分類的誤分和漏分情況都比較嚴(yán)重，原因可能是其在研究區(qū)內(nèi)巖性區(qū)范圍較小。

4 結(jié)論

1)以研究區(qū)30 m空間分辨率的DEM數(shù)據(jù)為基礎(chǔ)提取地形因子，通過比較分析地形因子區(qū)分巖性的有效性，得出地形因子中高程、剖面曲率、平面曲率、地表切割深度和地表粗糙度5個因子都有其最適宜識別區(qū)分的巖性，將5種因子用于巖性分類，得到的分類結(jié)果中主要巖性區(qū)均能被識別出來。

2)研究結(jié)果表明巖性與地形因子間存在一定的相關(guān)性，主要的巖性類都有其特有的地形特征都可以通過一定的地形因子組合而被識別出來。因此加入地形因子參數(shù)會明顯提高巖性的識別能力，這也為有植被覆蓋的巖性區(qū)解譯精度的提高提供了一種可能。為了深入研究各類巖性的地形特征從而進(jìn)一步提高巖性的解譯精度，還需挖掘巖性的更多地形信息，因此下一步需要研究更多的地形因子。

參考文獻(xiàn)(References)：

[1] 余海闊,李培軍.運用LANDSAT ETM+和ASTER數(shù)據(jù)進(jìn)行巖性分類[J].巖石學(xué)報,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] 王曉東.水系提取方法研究及其地質(zhì)意義[D].長春:吉林大學(xué),2015.

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

[3] 于亞鳳,楊金中,陳圣波,等.基于光譜指數(shù)的遙感影像巖性分類[J].地球科學(xué),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] 黃穎端,李培軍,李爭曉.基于地統(tǒng)計學(xué)的圖像紋理在巖性分類中的應(yīng)用[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] 曾德耀.基于最佳地形因子組合的地貌形態(tài)類型劃分研究[D].重慶:重慶交通大學(xué),2015.

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

[6] 楊晏立,何政偉,楊 斌,等.最佳因子復(fù)合的四川省地貌類型自動劃分[J].陜西理工學(xué)院學(xué)報(自然科學(xué)版)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圖像和地形因子的巖性單元分類——以新疆木壘地區(qū)為例[J].地球信息科學(xué)學(xué)報,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] 周啟明,劉學(xué)軍.數(shù)字地形分析[M].北京:科學(xué)出版社,2006:52-75.

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

[12] 劉少峰,王 陶,張會平,等.數(shù)字高程模型在地表過程研究中的應(yīng)用[J].地學(xué)前緣,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] 楊 昕,湯國安,劉學(xué)軍,等.數(shù)字地形分析的理論、方法與應(yīng)用[J].地理學(xué)報,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的我國地勢起伏度統(tǒng)計單元研究[J].水利學(xué)報,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].中國水土保持科學(xué),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.