基于遙感技術的南水北調水源區(qū)(河南段)石漠化遙感調查與評價

許軍強，張斌，袁晶，盧意愷

1.河南省航空物探遙感中心，鄭州450053；2.河南省地質科學研究所，鄭州450001；3.河南省地礦局第五地質勘查院，鄭州450001

0 引言

目前，遙感技術在巖溶石漠化信息提取、空間分析等領域得到了廣泛應用，并取得了一系列顯著成果。如陳希等[1]利用SAR數(shù)據(jù)結合光學影像融合對喀斯特地區(qū)石漠化等級劃分，提取石漠化信息，結果表明，多源遙感數(shù)據(jù)支持下的石漠化信息提取方法比傳統(tǒng)提取方法更加準確和有效。周迪[2]等基于影像紋理和地形數(shù)據(jù)輔助面向對象方法進行喀斯特地區(qū)石漠化信息的提取，可有效地減少因復雜地形導致石漠化信息提取結果“椒鹽化”現(xiàn)象，提取精度明顯優(yōu)于基于像素的監(jiān)督分類和非監(jiān)督分類方法。夏學齊等[3]提出“植被線”概念，并定義“石漠化幾何指數(shù)”(GRI)，客觀穩(wěn)定地提取石漠化程度。涂杰楠等[4-5]通過歸一化植被指數(shù)與植被覆蓋度進行巖溶石漠化遙感信息提取，得到兩期石漠化空間分布特征，并對兩期植被覆蓋度數(shù)據(jù)進行疊加，分析生態(tài)治理區(qū)石漠化的演變趨勢。劉芳等[6]基于通過對喀斯特地區(qū)巖石熱特性探究，提出基于熱紅外遙感的石漠化監(jiān)測方法，進行石漠化等級劃分并分析石漠化變化。甘敬民等[7]運用多因子綜合分析、專家打分和層次分析等方法，建立多因素石漠化強度等級評價及石漠化演化敏感性分析地理模型，得出石漠化等級及石漠化演化敏感性空間分布結果。熊康寧等[8]基于GIS/RS調查了貴州省巖溶石漠化的分布現(xiàn)狀和特征。陳飛等[9]利用CA-Markov預測方法，探尋研究區(qū)1990—2016年并預測了2016—2021年石漠化時空演變規(guī)律。郭麗琴等[10]選取TM和高分一號、資源三號和遙感二號2期衛(wèi)星遙感數(shù)據(jù)為數(shù)據(jù)源，采用RS和GIS技術，獲取了石漠化面積變化情況，探索石漠化時空演變特征及成因分析。李松等[11]采用遙感變化檢測的方法對TM/ETM+數(shù)據(jù)進行石漠化信息提取和精度檢驗，取得了較好的效果。李建存等[12]利用相隔10年的3期影像數(shù)據(jù)，進行碳酸鹽巖遙感解譯和信息提取及野外核查，生成碳酸鹽巖分布圖、石漠化分級分布圖以及石漠化趨勢圖等，為石漠化的防治提供依據(jù)。

1 生態(tài)環(huán)境背景

綜合丹江口水庫(河南轄區(qū))匯水區(qū)域和碳酸鹽巖的分布選取南水北調水源區(qū)(河南段)的范圍。南水北調水源區(qū)(河南段)主要位于河南省南陽市西南部，南陽盆地西緣，包括淅川南部、內鄉(xiāng)和鄧州西南部一隅。東鄰內鄉(xiāng)縣、鄧州市，北依西峽縣，西與陜西省交界，南與湖北省相接，地處32°35′～33°08′N， 111°08′～111°44′E之間，總面積約2 584 km2(圖1)。研究區(qū)屬北亞熱帶向暖溫帶過渡的季風性氣候區(qū)，年平均氣溫15.8℃。年降雨量在391.3～1 423.7 mm之間，多年平均降雨量為804.3 mm，多年月平均降雨量68.2 mm，區(qū)內雨量充沛。

圖1 研究區(qū)交通位置圖Fig.1 Traffic location map of research area

研究區(qū)地形為西北高、東南低，呈馬蹄狀。西北部地貌為侵蝕中山與低山區(qū)，中部為丘陵區(qū)，東南部為崗地與河流階地緩沖區(qū)。區(qū)內出露地層從老到新依次是震旦系、寒武系、奧陶系及石炭系，碳酸鹽巖分布廣泛，有灰?guī)r、白云質灰?guī)r、白云巖及大理巖等。廣泛分布的各類碳酸鹽巖為巖溶作用提供了物質條件，成為石漠化的主要物質基礎。

2 研究方法

2.1 波譜特征分析

圖2 碳酸鹽巖地區(qū)主要地物波譜曲線Fig.2 Spectrum curves of main objects in area composed of limestone

借鑒碳酸鹽巖地區(qū)主要地物波譜[13]，碳酸鹽巖地區(qū)主要地物類型為3種：①植被；②碳酸鹽巖；③土壤(圖2)。從圖中可以看出，植被的反射率在0.48 μm附近有一個18%的反射波峰，在0.55 μm附近有一個反射率為10%～15%的吸收波谷，即在TM3波段中紅光波段為植被吸收波段，在0.76～0.82 μm區(qū)間反射率陡增，在0.82～1.3 μm區(qū)間植被反射率相對穩(wěn)定在40%～50%的高反射區(qū)間，在1.5 μm附近植被反射率陡降到17%的吸收谷，在1.60 μm升到32%的次反射峰，說明TM4波段中的近紅外波段(0.76～0.90 μm)是植被強反射區(qū)間，可以區(qū)分植被類型和覆蓋率的差異。碳酸鹽巖波譜曲線相對平緩，在各個波段的反射率無明顯差異，TM1、TM2、TM3、TM4、TM5波段在0.40～1.75 μm的反射率基本穩(wěn)定在20%±，其中1.55～1.75 μm是絕大部分造巖礦物波譜響應的強反射波譜，碳酸鹽巖在TM5波段具有較強的識別能力，從1.75～2.35 μm反射率緩慢下降至17%±，TM7波段對造巖礦物的識別能力較弱。土壤的波譜曲線在0.40～2.10 μm區(qū)間，反射率整體呈上升趨勢，但波動較大，在0.40～0.80 μm區(qū)間，土壤反射率呈線性逐漸上升到42%，反射率達到一個波峰，0.80～0.90 μm區(qū)間逐漸下降到36%后又持續(xù)上升至46%，在TM5波段反射率達到最大53%。在TM7波段升降波動明顯。

石漠化在遙感圖像上的光譜響應特征可以通過波段DN值計算獲得，選取良好的石漠化訓練樣本對表征不同石漠化侵蝕程度的光譜特征以及石漠化分類能夠達到良好的效果。為保證樣本選取的準確性與可靠性，以野外踏勘選取的石漠化樣本點為主，投影到影像上選取典型的潛在石漠化、輕度石漠化、中度石漠化和重度石漠化的訓練樣區(qū)?；跇颖緟^(qū)石漠化光譜DN值，生成不同程度的石漠化光譜特征曲線(圖3)。

圖3 石漠化程度光譜特征Fig.3 Spectral characteristics of rocky desertification

圖3中，從均值的可分性看，重度石漠化均值最大，反射率最高，且與中度石漠化、輕度石漠化、潛在石漠化的反射率有顯著差異，中度石漠化均值其次，與輕度石漠化、潛在石漠化的反射率有明顯差異，輕度石漠化與潛在石漠化的反射率接近，可分性最弱。從波段特征上看，第4波段對重度石漠化和中度石漠化有較好的區(qū)分，第1波段對潛在石漠化識別較好，第2、3波段對中度和輕度石漠化有較好的分區(qū)，第5波段對重度石漠化區(qū)分較好。

2.2 石漠化信息提取方法

以TM/ETM+、OLI數(shù)據(jù)為主要數(shù)據(jù)源，對無石漠化和非石漠化做掩膜處理后的石漠化范圍進行石漠化程度信息提取，在不同程度石漠化訓練樣本的基礎上，基于改進增強型植被指數(shù)[14]和監(jiān)督分類的方法提取石漠化侵蝕程度信息。本文采用最大似然法、支持向量機和神經(jīng)網(wǎng)絡分類法，通過改進增強型植被指數(shù)與數(shù)據(jù)源的合并，參與分類，實現(xiàn)石漠化信息提取。

(1)最大似然分類

像元數(shù)據(jù)x，其對于給定的類別k的似然度Lk被定義為x出現(xiàn)時，其歸屬于類別k的條件概率的計算公式如下：

Lk=P(k|x)=P(k)×P(x|k)/∑P(i)×P(x/i)

(1)

式中：x為像元的觀測數(shù)據(jù)，它是由待分類像元n個特征變量的觀測值所組成的列向量。

對于同一分類類別和相同的訓練樣本，上式是個常量，可以簡化為：

Lk=P(k|x)=P(k)×P(x|k)

(2)

式中：先驗概率P(k)已知，求像元x與類別k的似然度，P(x|k)是唯一的未知變量。要計算x與類別k的似然度，必須先給出總體的概率密度函數(shù)P(x|k)。對于未知類別X，它屬于類別A的概率為Pa，它屬于類別B的概率Pb，且Pa>Pb，所以未知類別X屬于類別A。

(2)支持向量機分類

支持向量機是以結構風險最小化為原則，在線性分類的基礎上提出并應用到線性分類器中。

對非線性可分的訓練樣本點，用一個非線性函數(shù)ψ將訓練樣本點空間映射到一個高維的特征空間并進行線性分類。若可以找到映射ψ，則內積運算(xi·x)可用(ψ(xi) ·ψ(x))代替。用函數(shù)K(xi，x)=(ψ(xi)·ψ(x))表示內積運算而不直接利用函數(shù)ψ，稱函數(shù)為核函數(shù)。

最優(yōu)分類超平面可表示為：

(3)

在ENVI中的SVM經(jīng)常使用的核函數(shù)為徑向基核函數(shù)：

(4)

在ENVI軟件下，選擇菜單欄中的Classification模塊中的support Vector Machine，核函數(shù)類型選擇徑向基核函數(shù)。

(3)神經(jīng)網(wǎng)絡分類

神經(jīng)網(wǎng)絡是模仿人的大腦神經(jīng)元行為活動，具有分布式并行信息處理特征的一種數(shù)學模型。它是一個復雜而又龐大的系統(tǒng)，系統(tǒng)內存在大量的像人腦神經(jīng)元一樣的節(jié)點，節(jié)點之間相互連接，數(shù)據(jù)信息在神經(jīng)網(wǎng)絡中不斷的輸入和輸出，從而達到數(shù)據(jù)處理的目的。

BP人工神經(jīng)網(wǎng)絡算法的分類步驟描述如下：

① 對權值Wij置初值。對每一層的權值Wij的初始值設為較小的隨機數(shù)，其中Wi，n+1=-θ(θ為閾值)；

② 給定訓練樣本集，輸入訓練樣本和相應的目標輸出，即輸入向量X和期望輸出Y；

(5)

⑤ 更新權值Wij和閾值θ。其計算表達式為：

(6)

式中：η為學習速率，一般取值為0.2～0.6。

當使用動量時，更新權值Wij數(shù)學表達式可采用如下：

(7)

=Wij(t)-Wij(t-1)

(8)

式中：α為動量系數(shù)，0<α<1。

⑥ 在得出各層各個權系數(shù)之后，判斷誤差函數(shù)E是否收斂所給出的學習精度ε(E≤ε)，如果滿足學習精度要求，則算法結束；如果未滿足要求，則返回③繼續(xù)執(zhí)行。選擇ENVI4.8 Classification模塊下的Neural Net做神經(jīng)網(wǎng)絡分類，分類結果如圖4所示。

a.最大似然分類圖;b.支持向量機分類;c.神經(jīng)網(wǎng)絡分類圖。圖4 石漠化分類圖Fig.4 Rocky desertification classification

2.3 分類精度評價

評價石漠化強度分級結果的分類精度時，通過建立的解譯標志選取分類訓練樣本，把選取的樣本隨機分為兩部分，其中70%用于分類，30%用于精度評價，提高精度評價的可靠性，分類精度和Kappa系數(shù)如表1所示。

表1 分類精度評價

3 石漠化遙感監(jiān)測與演變分析

利用ArcGIS10.1空間疊加分析技術對兩期石漠化程度進行疊加分析處理，分別以1996、2006和2006年、2016年石漠化程度信息完成兩期的石漠化變化趨勢信息提取。基本了解1996—2006年和2006—2016年兩個動態(tài)時期內的石漠化侵蝕程度變化(表2、圖5)，掌握各石漠化程度的穩(wěn)定、改善或惡化情況(圖6)。從表2、圖5可以看出，1996—2006年的變化趨勢是以穩(wěn)定狀態(tài)為主，局部石漠化程度處于趨好的狀態(tài)，輕微改善面積為209.47 km2，明顯改善面積22.63 km2；2006—2016年的變化趨勢依然以穩(wěn)定狀態(tài)為主，生態(tài)地質狀況總體良好，經(jīng)野外調查發(fā)現(xiàn)，部分嚴重石漠化區(qū)域經(jīng)治理已有明顯改善情況，但局部石漠化依然有惡化現(xiàn)象。

表2 1996—2006年、2006—2016年二期石漠化面積變化趨勢統(tǒng)計表

Table 2 Statistic table of trend of rocky desertification area in 1996—2006 and 2006—2016

類型1996—2006年/km22006—2016年/km2嚴重加劇22.7375.15加劇48.27144.80穩(wěn)定1 323.621 189.29輕微改善209.47167.40明顯改善22.6330.00

a.1996—2006年；b.2006—2016年。圖5 石漠化變化趨勢柱狀圖Fig.5 Histogram of rocky desertification change trend from 1996—2006 and 2006—2016

a.1996—2006年；b. 2006—2016年。圖6 石漠化面積變化趨勢圖Fig.6 Chart of change trend of rocky desertification area

4 討論

在GIS平臺下，首先分別基于地層巖性、地形地貌、工程地質和水文地質等因素不同類別的確定性系數(shù)CF函數(shù)值來構建影響因素的柵格圖層，然后利用GIS平臺下空間分析模塊中的柵格計算式對各統(tǒng)計單元進行因子疊加分析。通過定義石漠化評價指數(shù)(stony desertification index，SDI)來進行評價石漠化的易發(fā)性。石漠化評價指數(shù)公式如下：

SDI=∑RiCFi

(9)

式中：SDI為評價單元的石漠化易發(fā)性指數(shù)；CFi為該單元各影響因素的CF值；Ri為根據(jù)層次分析法得到的各影響因素的權重。SDI也是一個處于區(qū)間[-1，1]的指標，SDI值越高，說明石漠化易發(fā)性越高，反之，SDI值越低，石漠化易發(fā)性越低。

分級方法采用GIS空間分析的柵格重分類工具，按照自然間距分類(natural break)方法，將石漠化易發(fā)性結果柵格圖分為4級，即石漠化高易發(fā)區(qū)(0.07～0.63)、石漠化易發(fā)區(qū)(-0.28～0.07)、石漠化低易發(fā)區(qū)(-0.64～-0.28)和石漠化不易發(fā)區(qū)(-0.86～-0.64)(圖7)。通過對研究區(qū)石漠化易發(fā)性綜合分區(qū)評價，石漠化易發(fā)性共分為高易發(fā)區(qū)、易發(fā)區(qū)、低易發(fā)區(qū)和不易發(fā)區(qū)4個分區(qū)，4個分區(qū)的空間分布和地層走向大體一致，說明地層巖性在石漠化易發(fā)性綜合分區(qū)評價中的影響最高，其權重系數(shù)在影響因子中最高也得以證明。

圖7 石漠化易發(fā)性綜合分區(qū)圖Fig.7 Comprehensive regionalization map of rocky desertification susceptibility

高易發(fā)區(qū)主要分布在研究區(qū)的中部和北部，中部的高易發(fā)區(qū)呈片狀特征，面積較大，此區(qū)域以奧陶系地層為主，巖性主要為白云質灰?guī)r和灰?guī)r等，北部的高易發(fā)區(qū)呈條帶狀特征，面積較小，此易發(fā)區(qū)域以寒武系、奧陶系地層為主，巖性主要為灰?guī)r、白云質灰?guī)r和含泥質條帶灰?guī)r等，高易發(fā)區(qū)的整體分布與該區(qū)域的巖性具有高石漠化發(fā)生率相吻合。易發(fā)區(qū)面積較大，主要集中在研究區(qū)北部，以石炭系地層為主，巖性主要為白云質灰?guī)r，同樣具有較高的石漠化發(fā)生率。低易發(fā)區(qū)面積較小，主要集中在研究區(qū)西北部、中部及西南部，分布地層有震旦系、志留系及寒武系地層，主要巖性為白云質大理巖、白云巖、泥灰?guī)r和砂礫巖等，此區(qū)域的石漠化發(fā)生率較低。不易發(fā)區(qū)的分布與第四系地層吻合，處于平原地貌，巖性為砂礫石和黏土等，石漠化發(fā)生率最低。

5 結論

(1)利用改進增強型植被指數(shù)和監(jiān)督分類相結合的方法提取石漠化侵蝕程度信息。通過采用最大似然法、支持向量機和神經(jīng)網(wǎng)絡分類均能有效提取研究區(qū)石漠化程度信息，最大似然法提取石漠化總體精度最高， 達到84.93， Kappa系數(shù)為0.792 3。

(2)利用空間疊加分析技術完成對兩期石漠化變化趨勢監(jiān)測。分別對1996、2006、2016年石漠化程度信息進行石漠化空間疊加分析。1996—2006年輕微改善面積為209.47 km2，明顯改善面積22.63 km2；2006—2016年的變化趨勢依然以穩(wěn)定狀態(tài)為主，生態(tài)地質狀況總體良好。

(3)基于石漠化評價指數(shù)進行石漠化易發(fā)性分區(qū)評價。以地層巖性、地形地貌、工程地質和水文地質等因素不同類別的不確定性系數(shù)及權重構建評價指數(shù)，得出石漠化高易發(fā)區(qū)、易發(fā)區(qū)、低易發(fā)區(qū)和不易發(fā)區(qū)。4個綜合分區(qū)指數(shù)分別為0.07～0.63、-0.28～0.07、-0.64～-0.28和-0.86～-0.64， 4個分區(qū)的空間分布和地層走向大體一致，表明該指數(shù)可定量表達石漠化易發(fā)性評價。