基于遙感生態(tài)指數(shù)的塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量評價

何亮， 杜清， 崔麗洋， 關(guān)力偉， 郭帥， 蔣磊*

(1.中國地質(zhì)調(diào)查局烏魯木齊自然資源綜合調(diào)查中心，新疆 烏魯木齊 830026； 2.新疆維吾爾自治區(qū)林業(yè)規(guī)劃設(shè)計院，新疆 烏魯木齊 830049)

改革開放以來，隨著我國經(jīng)濟飛速發(fā)展，人口數(shù)量也在不斷增加，導(dǎo)致人口、資源、環(huán)境三者之間的矛盾越來越顯著，生態(tài)環(huán)境問題已是各級政府和眾多學(xué)者們關(guān)注的重點。堅持綠水青山就是金山銀山理念，構(gòu)建人與自然和諧共存的生態(tài)文明體系已是我國“十四五”期間的發(fā)展戰(zhàn)略。因此，及時有效地監(jiān)測生態(tài)環(huán)境狀況，評價各類因子對生態(tài)環(huán)境質(zhì)量的影響，對決策者準(zhǔn)確把握生態(tài)環(huán)境現(xiàn)狀并采取針對性的治理措施有著重要意義[1]。

遙感技術(shù)在獲取大區(qū)域監(jiān)測數(shù)據(jù)方面具有速度快、成本低、效率高的優(yōu)勢，利用遙感技術(shù)來監(jiān)測林地[2]、濕地[3]、草地[4]、耕地[5]及城市[6]等生態(tài)系統(tǒng)已得到學(xué)術(shù)界廣泛認(rèn)可。目前，政府主要以生態(tài)環(huán)境狀況指數(shù)[7]來評價某區(qū)域的生態(tài)環(huán)境質(zhì)量，該方法獲得的EI指數(shù)介于0～100之間，只能用單一數(shù)值來表示區(qū)域生態(tài)環(huán)境狀況的好壞，無法表示區(qū)內(nèi)不同生態(tài)環(huán)境狀況的空間分布情況[8]。徐涵秋[8]提出以與人類生存息息相關(guān)且最能反映生態(tài)質(zhì)量，同時也方便獲取的綠度指標(biāo)、濕度指標(biāo)、干度指標(biāo)和熱度指標(biāo)為基礎(chǔ)，可構(gòu)建新型遙感生態(tài)指數(shù)(RSEI)，該指數(shù)不但與EI具有相同的量化評分，而且能表達區(qū)域生態(tài)環(huán)境質(zhì)量在時間、空間上的變化情況。自遙感生態(tài)指數(shù)提出后，諸多研究表明[9-13]，該方法能夠快速、有效地評價高植被覆蓋區(qū)的生態(tài)環(huán)境質(zhì)量，但未見該方法運用在低植被覆蓋區(qū)及干旱區(qū)的生態(tài)環(huán)境質(zhì)量評價中。

第二師塔里木墾區(qū)作為塔里木河下游綠色長廊的保衛(wèi)者，自流域水利管理體制改革后，塔里木河下游生態(tài)輸水量得到保障，塔里木墾區(qū)生態(tài)環(huán)境的變化更加受到當(dāng)?shù)卣蛯W(xué)者們的關(guān)注。從塔里木河管理局公布數(shù)據(jù)發(fā)現(xiàn)，塔里木河生態(tài)輸水量從2010年開始下降，2015年后逐年提高。因此，選擇2010年、2015年生態(tài)輸水改善前，以及2020年生態(tài)輸水改善后的三期遙感數(shù)據(jù)，利用RSEI模型評價塔里木墾區(qū)近十年生態(tài)環(huán)境的變化情況，以期為塔里木墾區(qū)農(nóng)業(yè)發(fā)展和生態(tài)格局規(guī)劃提供參考，同時對遙感生態(tài)指數(shù)在干旱地區(qū)的推廣應(yīng)用提供依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

塔里木墾區(qū)位于塔里木盆地東部，塔克拉瑪干沙漠東北部(86°34′—88°00′E,40°30′—41°05′N)，北倚孔雀河，南瀕塔里木河，地勢平坦，海拔大部分在855 m以下，西高東低，南高北低，主要由河流沖積平原以及大片風(fēng)蝕地組成，荒漠地帶多沙丘和古河道分布。墾區(qū)屬大陸性荒漠氣候，日照豐富，年平均氣溫10.8 ℃，最高氣溫7月，最低氣溫1月，年均降水量34.7 mm，最大降水量為96.8 mm，平均年蒸發(fā)量為2 408.6 mm。當(dāng)?shù)刂脖簧L稀疏，以胡楊、沙棗、梨樹、紅柳、梭梭、黑枸杞、蘆葦、駱駝刺為主。塔里木墾區(qū)隸屬新疆生產(chǎn)建設(shè)兵團第二師管轄，轄區(qū)內(nèi)設(shè)立3個團場，總面積2 207.69 km2，區(qū)內(nèi)建有兩個水庫，恰拉水庫位于墾區(qū)西北部，大西海水庫位于墾區(qū)中南部，水源補給為孔雀河和塔里木河。

1.2 數(shù)據(jù)來源

為構(gòu)建RESI指數(shù)，研究選取三期Landsat遙感數(shù)據(jù)，其中2010年為Landsat5 TM數(shù)據(jù)，2015年、2020年為Landsat8 OLI數(shù)據(jù)。塔里木河生態(tài)輸水大多發(fā)生在每年8月左右，如2010年8月，2015年7月，2020年8月，此月份中塔里木墾區(qū)水域面積較大，植被長勢最好，生態(tài)環(huán)境質(zhì)量達到峰值。三期數(shù)據(jù)獲取時間分別為2010年8月13日、2015年7月26日、2020年8月8日。遙感數(shù)據(jù)來自美國地質(zhì)調(diào)查局官網(wǎng)(http://earthexplorer.usgs.gov/)以及地理空間數(shù)據(jù)云(http://www.gscloud.cn/)，大氣剖面數(shù)據(jù)來自NASA官網(wǎng)(https://atmcorr.gsfc.nasa.gov/)。數(shù)據(jù)預(yù)處理方面，采用 ENVI 5.3 軟件進行輻射定標(biāo)、大氣校正、鑲嵌、裁剪等，得到研究區(qū)三期Landsat反射率數(shù)據(jù)。

1.3 研究方法

徐秋涵[8]提出的RSEI模型包含4個與人類生存環(huán)境密切相關(guān)的指標(biāo)，分別為綠度指標(biāo)、濕度指標(biāo)、干度指標(biāo)和熱度指標(biāo)，為方便數(shù)據(jù)獲取和評價，通常用歸一化植被指數(shù)(NDVI)代表綠度指標(biāo)，纓帽變換后的濕度分量(WET)[14,15]代表濕度指標(biāo)，裸土指數(shù)(SI)[16]和建筑指數(shù)(IBI)[17]共同組成干度指標(biāo)(NDBSI)，熱輻射反演地表溫度(LST)代表熱度指標(biāo)，其表達式為：

RSEI=f(NDVI,WET,NDBSI,LST)

(1)

式中：RSEI表示遙感生態(tài)指數(shù)；NDVI為歸一化植被指數(shù)；WET為濕度分量；NDBSI為干度指標(biāo)；LST為地表溫度。

1.3.1 綠度指標(biāo)

歸一化植被指數(shù)(NDVI)能夠較好地反映區(qū)域內(nèi)植物生長狀態(tài)及分布規(guī)律，故采用NDVI來代表綠度指標(biāo)，其公式為：

NDVI=(ρnir-ρred)/(ρnir+ρred)

(2)

式中：NDVI為歸一化植被指數(shù)；ρnir為近紅外波段；ρred為紅光波段。

1.3.2 濕度指標(biāo)

以Landsat TM和OLI數(shù)據(jù)為基礎(chǔ)，采用Crist[14]和Baig[15]提出的計算方法，通過纓帽變換得到濕度分量，計算公式如下：

WETTM=0.0315ρ1+0.2021ρ2+0.3102ρ3+

0.1594ρ4-0.6806ρ5-0.6109ρ7

(3)

WETOLI=0.1511ρ2+0.1973ρ3+0.3283ρ4+

0.3407ρ5-0.7117ρ6-0.4559ρ7

(4)

式中：ρi(i=1,…,6,7)為TM和OLI數(shù)據(jù)各波段的反射率。

1.3.3 干度指標(biāo)

地表干度表現(xiàn)在兩方面，一是裸土，二是人工建筑，因此采用裸土指數(shù)(SI)和建筑指數(shù)(IBI)平均的方法代表干度指標(biāo)(NDBSI)，計算公式為：

(5)

(6)

NDBSI=(SI+IBI)/2

(7)

式中：ρblue為藍光波段；ρgreen為綠光波段；ρred為紅光波段；ρnir為近紅外波段；ρswir1為短波紅外波段。

1.3.4 熱度指標(biāo)

熱度指標(biāo)的直觀表現(xiàn)是地表溫度，本文采用熱輻射傳導(dǎo)方程來反演地表溫度，具體步驟如下：

(1)計算植被覆蓋度，公式為：

Fv=(NDVI-NDVIs)/(NDVIv-NDVIs)

(8)

式中：Fv為植被覆蓋度；NDVI為歸一化植被指數(shù)；NDVIs為無植被覆蓋的NDVI值；NDVIv為完全植被覆蓋的NDVI值。本文根據(jù)研究區(qū)植被特點，取NDVIs=0.05，NDVIv=0.70。

(2)計算地表比輻射率，水體的比輻射率取值0.995 0[18]，植被和建筑的比輻射率計算公式如下：

εplant=0.9625+0.0614Fv-0.0461Fv2

(9)

εbuilding=0.9589+0.0860Fv-0.0671Fv2

(10)

(3)計算黑體輻射亮度值，公式如下：

B(Ts)=[Lλ-L↑-τ(1-ε)L↓]/τε

(11)

式中：B(Ts)為黑體在Ts溫度時的熱輻射亮度值；Lλ為熱紅外波段的輻射亮度值；L↑、L↓分別為大氣向上、向下輻射亮度值；τ為大氣在熱紅外波段的透過率，ε為比輻射率。L↑、L↓、τ數(shù)值取自大氣剖面官網(wǎng)。

(4)根據(jù)普朗克公式計算地表溫度：

LST=K2/ln(K1/B(Ts)+1)

(12)

式中：LST為地表溫度；K1、K2為Landsat衛(wèi)星定標(biāo)參數(shù)；對于TM：K1=607.76 W/(m2·um·sr)，K2=1 260.56 K；對于OLI：K1=774.89 W/(m2·um·sr)，K2=1 321.08 K。

1.3.5 遙感生態(tài)指數(shù)構(gòu)建

對4個指標(biāo)進行主成分分析，可有效避免主觀意識對指標(biāo)權(quán)重的干擾，使結(jié)果精確可靠[9]。為避免水庫、湖泊等大面積水域?qū)χ鞒煞州d荷分布的影響，利用歸一化水體指數(shù)[19]將水域信息從4個指標(biāo)中去除。由上文公式計算得到的4個指標(biāo)量綱不同，需將其正規(guī)化，公式如下：

(13)

式中：NIi為正規(guī)化后的某指標(biāo)；Ii為指標(biāo)在像元i的值；Imax為該指標(biāo)的最大值；Imin為該指標(biāo)的最小值。

正規(guī)化后的4個指標(biāo)通過ENVI波段融合形成一幅含有4個波段信息的遙感數(shù)據(jù)，對其進行主成分分析，得到第一主成分(PC1)。當(dāng)植被覆蓋度高、生態(tài)質(zhì)量好的區(qū)域PC1為負(fù)值時，為使PC1值越高代表生態(tài)質(zhì)量越好，需用1減去PC1得到初始RSEI0，公式如下：

RSEI0=1-PC1[f(NDVI,WET,NDBSI,LST)]

(14)

在自然環(huán)境中，水域面積是衡量一個區(qū)域生態(tài)質(zhì)量好壞的重要指標(biāo)。由上式得到RSEI0缺少水域信息，為更完整地評估生態(tài)環(huán)境質(zhì)量，將水域的MNDWI進行正規(guī)化處理，壓縮到0.8～1的區(qū)間內(nèi)，通過ArcGIS中的鑲嵌功能填補RSEI0中缺失的水域信息。計算公式如下：

RSEIwater=MNDWI(1-0.8)+0.8

(15)

為了方便不同時期遙感生態(tài)指數(shù)的量化對比，對RSEI0進行正規(guī)化處理：

(16)

2 結(jié)果與分析

2.1 各指標(biāo)結(jié)果分析

根據(jù)PCA統(tǒng)計，結(jié)果如表1。三期生態(tài)指標(biāo)中，PC1的特征貢獻率均超過93%，說明PC1包含了4個指標(biāo)的主體信息，其中NDVI和WET為負(fù)值，NDBSI和LST為正值，說明研究區(qū)PC1中以干度、熱度指標(biāo)為主，且干度指標(biāo)對PC1的影響大于熱度指標(biāo)，原因是塔里木墾區(qū)超過一半的面積為干旱炎熱的荒漠地。綠度、濕度指標(biāo)與PC1呈負(fù)相關(guān)，說明綠度指標(biāo)、濕度指標(biāo)不是研究區(qū)的主體，且綠度指標(biāo)對PC1的影響大于濕度指標(biāo)，原因是塔里木墾區(qū)植被覆蓋度較低，濕地面積小。

表1 不同年份4個生態(tài)指標(biāo)PCA結(jié)果Table 1 PCA analysis of four ecological indexes in different years

計算不同年份各指標(biāo)的均值，結(jié)果如表2。綠度和濕度指標(biāo)的均值逐年增長，2020年達到最高，分別為0.395和0.635，說明塔里木墾區(qū)植被覆蓋面積和土壤濕度逐年增加。干度指標(biāo)呈現(xiàn)逐年降低的趨勢，這與兵團開荒種田、滴灌技術(shù)運用密切相關(guān)。熱度指標(biāo)變化幅度不大，保持在0.673左右。從遙感生態(tài)指數(shù)來看，無水體信息的RSEI指數(shù)均值分別為0.287、0.277和0.303，通過ArcGIS鑲嵌水體信息后的RSEI指數(shù)分別為0.301、0.292和0.315，表明塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量處于較差等級。RSEI數(shù)值呈現(xiàn)先降低后升高、總體升高的趨勢，說明塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量略有改善。對比兩種方法，鑲嵌水體信息的RSEI均值都高于無水體信息的均值，說明將MNDWI壓縮、鑲嵌至原始RSEI，該方法對研究區(qū)生態(tài)環(huán)境質(zhì)量評價具有積極影響。

表2 不同年份各指標(biāo)歸一化均值結(jié)果Table 2 Results of normalized mean and standard deviation of each index in different years

2.2 生態(tài)環(huán)境分區(qū)評價

根據(jù)徐涵秋[8]提出的等級劃分標(biāo)準(zhǔn)，將RSEI分為極差(0～0.2)、較差(0.2～0.4)、一般(0.4～0.6)、良好(0.6～0.8)、優(yōu)秀(0.8～1)這5個等級，結(jié)果如圖1。塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量總體較差，極差和較差的RSEI分布面積廣，主要地類為沙漠和荒地。西部生態(tài)環(huán)境質(zhì)量明顯高于東部，因為西部恰拉水庫常年蓄水，人口密集，耕地面積大，農(nóng)業(yè)發(fā)展條件優(yōu)越，且近年來在政府的大力支持下，梨、棗果業(yè)發(fā)展迅猛，不僅促進了當(dāng)?shù)亟?jīng)濟發(fā)展，還提供了防風(fēng)固沙生態(tài)功能，改善了當(dāng)?shù)厣鷳B(tài)環(huán)境質(zhì)量。東部34團場由于沙漠占地面積大，風(fēng)沙侵蝕嚴(yán)重，人口較少，且大西海水庫蓄水能力不足，導(dǎo)致該區(qū)域植被長勢較差，農(nóng)業(yè)發(fā)展受限。

圖1 塔里木墾區(qū)各期RSEI等級分布圖Figure 1 Distribution of RSEI grades in different periods in Tarim Reclamation Area

計算RSEI各等級所占面積，結(jié)果如表3。2010—2020年，塔里木墾區(qū)不同生態(tài)環(huán)境等級所占面積變化差異明顯。極差等級所占面積為50%左右，且呈現(xiàn)逐年增加的趨勢，由2010年的48.69%增長到2020年的58.90%，增幅達到20.96%，原因是塔里木墾區(qū)位于兩大沙漠之間，氣候干旱，水源不足，導(dǎo)致該地區(qū)土地沙化嚴(yán)重。較差、一般、良好等級占地面積呈減少趨勢，其中較差、良好等級占地面積減幅逐漸減小，一般等級占地面積減幅保持穩(wěn)定。優(yōu)秀等級占地面積從2010年的5.73%增長到2020年的14.62%，一方面是塔里木河生態(tài)輸水量增加，從2010年的3.6×108m3增加到2020年的4.6×108m3，墾區(qū)水庫面積從60.45 km2增加到89.74 km2；另一方面是開荒及荒地復(fù)墾增加了耕地面積，耕地面積從415.50 km2增加到439.10 km2。

表3 塔里木墾區(qū)不同等級RSEI面積變化表Table 3 Area changeTable of different grade RSEI in Tarim Reclamation Area

2.3 生態(tài)環(huán)境時空變化分析

將不同時期的RSEI做空間差值，以數(shù)值表示生態(tài)環(huán)境變化情況，將其分為5類：惡化(<-0.6)、變差(-0.6～-0.2)、不變(-0.2～0.2)、改善(0.2～0.6)、優(yōu)化(>0.6)，結(jié)果如圖2。2010—2015年間，生態(tài)環(huán)境改善的區(qū)域分布在團場內(nèi)部農(nóng)作區(qū)，優(yōu)化的區(qū)域集中在大西海水庫南側(cè)；生態(tài)環(huán)境變差和惡化的區(qū)域較小，變差的區(qū)域呈星點狀分布在團場耕作區(qū)及恰拉水庫一帶，惡化的區(qū)域集中在大西海水庫西側(cè)，環(huán)境惡化的原因是2015年塔里木河流量減小，水庫蓄水量低。2015—2020年間，生態(tài)環(huán)境變化劇烈的區(qū)域集中在水庫周邊，恰拉水庫大面積生態(tài)環(huán)境改善，大西海水庫北側(cè)生態(tài)環(huán)境改善明顯，而南側(cè)發(fā)生大面積惡化，原因是2020年生態(tài)輸水效用顯著，恰拉水庫和大西海水庫蓄水量增加，由于大西海水庫南側(cè)維修加固，使水庫干涸。十年間，塔里木墾區(qū)生態(tài)環(huán)境優(yōu)化的區(qū)域集中在大西海水庫中部，改善和變差的區(qū)域星點狀分布在團場農(nóng)作區(qū)，惡化的區(qū)域集中在大西海水庫西側(cè)。

圖2 塔里木墾區(qū)2010—2020年RSEI變化監(jiān)測圖Figure 2 Variation of RSEI in Tarim Reclamation Area from 2010—2020

計算各類變化所占面積，結(jié)果如表4。2010—2015年間，塔里木墾區(qū)生態(tài)環(huán)境改善和優(yōu)化的面積分別占到6.57%和0.20%，生態(tài)環(huán)境變差和惡化的面積達到2.59%和0.09%。2015—2020年，改善和優(yōu)化的面積分別占到4.01%和0.50%，生態(tài)環(huán)境變差和惡化的面積達到1.82%和0.15%。十年間，塔里木墾區(qū)生態(tài)環(huán)境改善和優(yōu)化的面積占比9.10%，變差和惡化的面積為2.32%，說明塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量有所改善，這與塔里木河生態(tài)輸水工程密切相關(guān)。

表4 塔里木墾區(qū)RSEI變化監(jiān)測表Table 4 Monitoring Table of RSEI change in Tarim Reclamation Area

3 討論

(1)對塔里木墾區(qū)4個指標(biāo)進行主成分分析，綠度、濕度指標(biāo)與PC1呈負(fù)相關(guān)，干度、熱度指標(biāo)與PC1呈正相關(guān)，其中干度指標(biāo)的貢獻度最大，其次是熱度、綠度、濕度，說明干度指標(biāo)是影響塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量的主要因子。主成分分析得到的結(jié)果與南方高植被覆蓋區(qū)的研究結(jié)果[9,12-13]相反，原因是塔里木墾區(qū)荒漠占地面積大，植被覆蓋度較低。

(2)通過查閱相關(guān)文獻，基于遙感生態(tài)指數(shù)的評價方法大多運用在高植被覆蓋區(qū)，在干旱區(qū)未見運用。本研究選擇典型干旱區(qū)運用遙感生態(tài)指數(shù)模型進行生態(tài)環(huán)境質(zhì)量評價，結(jié)果顯示塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量總體處于較差等級，RSEI保持在0.3左右，較差和極差區(qū)域占比達到72%，生態(tài)環(huán)境優(yōu)秀、良好和一般的區(qū)域集中在恰拉水庫、大西海水庫和團場農(nóng)業(yè)區(qū)。RSEI呈現(xiàn)先降低后升高的趨勢，生態(tài)環(huán)境改善、優(yōu)化的面積大于變差、惡化的面積，說明塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量有所改善。以上評價結(jié)果均有數(shù)據(jù)支撐，說明遙感生態(tài)指數(shù)模型在干旱區(qū)適用。

(3)基于前人研究，RSEI模型在使用中為避免水體對主成分分析的影響，通常刪除了水體區(qū)域，而水體對生態(tài)環(huán)境質(zhì)量的貢獻卻是不可忽視的。本研究對水域內(nèi)的MNDWI進行正規(guī)劃處理，將其壓縮到0.8～1區(qū)間內(nèi)，通過ArcGIS鑲嵌功能填補RSEI缺失的水域信息。結(jié)果表明，鑲嵌水域信息的RSEI值高于一般方法得到的結(jié)果，且鑲嵌后的RSEI在空間分布上完整，便于多期數(shù)據(jù)空間分析。但是，將正規(guī)化后的MNDWI值壓縮至0.8～1，區(qū)間閾值選擇是否合適，有待進一步研究證明。

4 結(jié)論

選取2010年、2015年和2020年三期Landsat遙感數(shù)據(jù)，通過RSEI模型展現(xiàn)塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量分布及變化情況。結(jié)果表明，塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量總體處于較差等級，較差和極差面積占比達到72%，生態(tài)環(huán)境質(zhì)量優(yōu)秀、良好、一般的區(qū)域集中在恰拉水庫、大西海水庫和團場農(nóng)業(yè)區(qū)，周邊地區(qū)生態(tài)環(huán)境質(zhì)量較差。RSEI呈現(xiàn)先降低后升高的趨勢，生態(tài)環(huán)境改善、優(yōu)化的面積大于變差、惡化的面積，說明塔里木墾區(qū)生態(tài)環(huán)境質(zhì)量得到改善，近年來環(huán)境治理得到一定成效，但不能放松警惕，需積極進行生態(tài)環(huán)境規(guī)劃治理，處理好農(nóng)業(yè)發(fā)展和生態(tài)保護之間的關(guān)系。