引黃濟青工程沿線生態因子1990-2020年時空變化特征

馬吉剛， 湯宇婷， 張立民， 吳忠勝， 崔子騰， 竇 智

(1.山東省調水工程運行維護中心 山東 濟南 250000； 2.河海大學 地球科學與工程學院 江蘇 南京 211100； 3.山東省調水工程運行維護中心昌邑管理站 山東 濰坊 261000；4.青島引黃濟青水務有限公司 山東 青島 266000； 5.山東省調水工程運行維護中心壽光管理站 山東 濰坊 261000)

跨流域調水工程區域內的生態環境變化，本質上是生物、土壤、氣候、地形等生態因子的演變[1]。因此，基于對水生態環境質量及生態因子演變的研究，有利于工程對生態影響的分析，更有助于環境問題的對癥治理。遙感技術作為調水工程水生態環境研究中的重要手段，具有多空間性[2]、多時間性[3]、高效性[1]。不僅能對水生態環境進行有效的實時監測，還便于展開長期動態工作，快速、準確地反映大面積水環境的生態質量現狀及發展趨勢。目前，國內外不少學者對遙感技術在調水工程中的應用作了深入的研究。包洪福[4]以南水北調中線工程丹江口水源地為研究區，基于遙感技術對區域內水生態環境中的動植物變化規律進行了研究，探究其生物多樣性。Yang等[5]使用無人機技術結合遙感數據對南水北調工程地下水位上升引起的生物量改變和緩沖帶寬度進行了安全評估。Gao等[6]結合了geographic information system (GIS)和remote sensing (RS)技術，以生態系統服務價值法為手段，對江蘇省東部干線水源地的生態系統服務價值進行了計算并對該地的生態補償情況展開了進一步研究。朱長明等[7]綜合地表水、地下水等數據和多源遙感長時間序列數據，對塔里木河下游生態輸水前后區域水環境變化和生態響應進行監測和分析。

引黃濟青調水工程通水以來，有效緩解了青島及沿線城市水資源供需不均的矛盾[8]。但在其建設和運行過程中，也對沿線的水生態環境產生了一定影響[9]。本文以工程沿線為研究區，基于遙感技術對土地利用/覆蓋、植被指數及覆蓋度、地表溫度、土壤干濕度等生態因子的時空變化規律展開討論，進一步探究遠距離調水工程沿岸生態環境的演變規律，為沿線生態的保護和治理提供了數據支撐和科學依據。

1 研究區與數據

1.1 研究區概況

引黃濟青工程是中國“七五”期間重要的水利工程，位于118°07′—120°12′E，36°22′—37°16′N的山東半島[10]，由山東省濱州市打漁張閘引黃河水到青島市嶗山區白沙河凈水廠。該地區地勢平坦，受海洋影響，屬暖溫帶季風氣候區，具有大陸性氣候的明顯特征，年平均氣溫12.0～12.6 ℃，年平均降雨量在600～770 mm，常年平均蒸發量在1 000 mm以上，且起伏變化較大。地帶性土壤為典型棕色森林土，分布在緩坡地和排水良好的平地。天然植被類型為暖溫帶落葉闊葉林，主要樹種有櫟類，以麻櫟居多。

調水工程線路全長291.14 km，其中棘洪灘水庫以上為水源工程，長253.14 km，以明渠為主，水庫以下為供水工程，長38 km，以管道為主。本文以水庫以上水源工程輸水段為研究對象，綜合考慮研究結果的代表性、準確性，建立了干渠向外左右各2 km的緩沖區，以此展開討論。

1.2 數據來源

研究區范圍廣闊、地形復雜，適合衛星遙感進行大面積、多角度、高時效的生態環境監測。考慮到研究周期和影像精度的綜合要求，最終選用Landsat系列影像作為初始數據。

工程通水至今已有30 a余，為研究其沿程生態因子的歷史演變規律，選擇以4個時期作為時間研究對象，分別是1990，2000，2010，2020年。此外，考慮到春夏兩季植被生長旺盛，識別度高，而夏季云層覆蓋度廣，可見度低等多種因素，最終選擇了5,6月日期較接近的八景影像進行數據的處理和分析，具體影像資料見表1。

表1 研究區遙感影像基本信息

為解決光學遙感數據獲取時的誤差問題，需先對下載好的影像進行預處理，包括了輻射定標和大氣校正過程。其中輻射定標是為了消除傳感器本身誤差，大氣校正是為了消除大氣散射、吸收、反射引起的誤差[11]。校正后的影像再根據研究區的范圍進行裁剪和鑲嵌處理。

2 研究方法

2.1 技術路線

將研究區按地級市范圍劃分為4個區段，分別為濱州市、東營市、濰坊市及青島市區段，基于生態遙感監測技術，選擇了土地利用/覆蓋、植被指數及覆蓋度、地表溫度、土壤干濕度等生態因子，探究其時空演變規律，討論工程運行后沿線生態環境的變化。

2.2 土地利用/覆蓋的量化

土地利用/覆蓋中綠地、水域、建筑用地、裸地面積的提取，是在遙感影像經過預處理的基礎上，利用ENVI軟件中監督分類的方式進行劃分的[12]，結合影像和地圖對分類樣本進行多次調整直至樣本間具有較好的可分離性。分離器的選擇也是同時對比了多種分離試驗結果，最終選擇了支持向量機法進行分類。分類結果再以高分辨率的影像作為參照，對混淆部分進行人工目視解譯，解決同物異譜或異物同譜帶來的誤差問題[13]，最后利用ENVI軟件中的混淆矩陣精度評價法進行驗證，并對比“地理國情監測云平臺”(http:∥www.dsac.cn/)相關數據產品，分類結果的精度達到80%以上。

2.3 植被的量化

植被的量化采用的是NDVI指數法[14]，通常反映植被覆蓋、生長等信息，其計算公式為：

式中：NIR為近紅外波段的反射率；R為紅波段的反射率。

在現有NDVI值的基礎上，進一步分析植被的覆蓋情況，這里采用了像元二分法[15]計算得到植被覆蓋度PV[16]，其具體計算模型如下：

式中：PV為植被覆蓋度； NDVI為植被指數； NDVISoil為完全是裸土或無植被覆蓋區域的NDVI值； NDVIVeg為完全被植被所覆蓋的像元的NDVI值。

2.4 地表溫度的量化

本文采用大氣校正法[17](輻射傳輸方程)進行地表溫度的反演，其本質是通過消除大氣對地表熱輻射的影響誤差，將得到的地表熱輻射強度轉化為相應的地表溫度。衛星傳感器接收到的熱紅外輻射亮度值Lλ由3部分組成： ①大氣向上輻射亮度； ②地面的真實輻射亮度經過大氣層后到達衛星傳感器的能量； ③大氣向下輻射到達地面反射的能量，其計算公式如下：

Lλ=〔εB(Ts)+(1-ε)L↓〕τ+L↑

(3)

式中：ε為地表比輻射率；Ts為地表真實溫度(K)；B(Ts)為黑體熱輻射亮度；τ為大氣在熱紅外波段的透過率；L↓為大氣下行輻射亮度；L↑為大氣上行輻射亮度。則溫度為T的黑體的輻射亮度B(Ts)可表示為：

式中：τ，L↑，L↓參數可在NASA提供的網站上(http:∥atmcorr.gsfc.nasa.gov/)通過提供時間、經緯度等數據信息進行獲取。而地表比輻射率ε則采用NDVI閾值法計算，公式如下：

ε=0.004PV+0.986

(5)

獲得所需數據參數后，利用普朗克公式計算得到地表真實溫度Ts，具體公式如下：

式中：K1，K2為定標系數，可通過影像的元數據獲取。

結合沿線濱州、東營、濰坊、青島市氣象站點收集到的實測均溫數據，對地表溫度的反演結果進行驗證，結果顯示兩者波動變化一致，說明方法具有可行性。

2.5 土壤干濕度的量化

土壤干濕度的量化方法選用的是NDVI即溫度植被干旱指數法[18]，其原理是在光學與熱紅外遙感通道數據的基礎上反演植被覆蓋土層表面的水分情況，本質是結合歸一化植被指數(NDVI)和地表溫度(Ts)兩個參數，構建Ts-NDVI特征空間，其公式表示為：

式中：Ts為地表溫度值；Tsmin為NDVI相同值所對應的最低地表溫度值，為Ts-NDVI特征空間的濕邊；Tsmax為NDVI相同值所對應的最高地表溫度，為特征空間的干邊。根據像元構造的特征空間，同時對干邊和濕邊進行線性擬合，其方程為：

Tsmax=a1+b1×NDVI

(8)

Tsmin=a2+b2×NDVI

(9)

式中：a1，b1為干邊擬合方程的系數；a2，b2為濕邊擬合方程的系數。如此TVDI指數的計算方法可以進一步轉化為：

TVDI的值在0-1之間，干邊對應的TVDI值為1，濕邊對應的TVDI值為0，TVDI的值越趨向于0，土壤的濕潤度越高，反之，土壤的濕潤度越低。

3 結果與分析

3.1 土地利用/覆蓋

基于監督分類和目視解譯結合的方法將研究區內土地分為了：綠地、水域、建筑用地、裸地4大類。其中，綠地包括了林地、草地、耕地等；水域涵蓋了江、河、湖、渠等；建筑用地指城鄉住宅、公共用地等；裸地一般為無植被生長的未利用土地，總體空間分布見圖1。為探究各地類的變化規律，進一步統計了1990，2000，2010，2020年各類土地利用/覆蓋信息的面積及比例情況如表2。

圖1 研究區1990，2000，2010，2020年土地利用分類

表2 研究區土地利用分類面積及比例統計

由圖1可知，在時間維度上，1990年研究區內土地結構以綠地和裸地為主，建筑用地較少，水域地區更為稀缺，說明工程沿線土地資源豐富，但城鎮化水平較低，城市基礎設施建設較落后。1990—2010年期間，綠地面積變化幅度不大，相比之下建筑用地面積增勢明顯，同時裸地區域也在以肉眼可見的速度減少。2020年區域內地類主要以綠地和建筑用地為主，裸地區域被大面積開墾利用，可見當地城市化建設發展迅速。

在空間維度上，按研究區內4個地級市區段劃分并對各區段地類進行統計，進一步分析得到濱州和青島區段的水域面積較其他區段相對較大。可能受其地理位置的影響，濱州地處引黃濟青工程黃河入水口，而青島東、南瀕臨黃海，帶來的河湖水源較充足。同時，水源和地理的優勢也使得青島區段的耕地、林地、草地發展的較好，其綠地面積比例在4個區段中一直都居于前列。建筑用地區域在4個區段中區別并不大，說明研究區內各地區城鎮化發展相對較為均衡。

從表2數據分析，1990—2020年，綠地面積比例共增長6.93%，長期占據研究區內一半以上土地；建筑用地面積比例增加18.37%，尤其以1990—2000年增勢明顯，增幅近15%。由此可見，這10 a的城鎮化發展效果顯著。相比之下水域面積增加并不明顯(30 a間面積增長17.18 km2)，增幅1.61%，大部分河流湖泊位于濱州和青島地段，主要受黃河和黃海影響。當然調水工程作為沿線重要的生態紐帶也會直接促進水資源面積的增長。此外，城市的發展離不開城鄉建設，工程運行以來，大面積的未利用土地被開墾翻新以開展房屋、基礎設施建設，相對的裸土地區域大幅縮減，占總面積比例降低近27%，現階段研究區內僅有2.61%的區域為裸土區。

3.2 植被指數及覆蓋度

3.2.1 植被指數 植被指數的值一般在[-1,1]區間，通常植被長勢越好、覆蓋程度越高對應的指數越高。整個研究區的NDVI空間分布見圖2。

圖2 研究區1990，2000，2010，2020年NDVI空間分布特征

利用均值法統計了不同時間研究區各地市級區段的平均NDVI指數(如圖3所示)。引黃濟青工程沿線各個時期的NDVI均值分別為0.296 2,0.335 7,0.430 1,0.487 4。從整體變化趨勢來看，指數呈現不斷上升趨勢。在1990年植被指數處于較低的狀態，說明當時的研究區內植被生長并不茂盛，植被的覆蓋率也比較低。經過30a多的發展，植被指數均值增加了近0.2，增勢明顯，可見當地的生態發展進步顯著。

圖3 總研究區及各地級市區段NDVI均值變化趨勢

通過對圖3的進一步分析，發現各地市級區段的NDVI均值變化趨勢與沿線總趨勢基本保持一致。其中，青島區段的均值變化較為顯著，大體上都高于其他地區，前文曾提及青島區段中的水域面積也居于前列，而水資源對于植被作物的生長不可或缺。加之青島市地理位置優越，氣候適宜，對草木的生長和農作物的培育更是十分有益。此外，可以看到在1990年，濱州、濰坊、青島3區段的NDVI均值幾乎接近，而東營區段的植被指數卻遠低于3者，說明當時東營區段內雖然耕地、林地、草地的占地面積不少，但植被的茂密程度相對于其他區域較弱。到了2020年，4地市級區段的植被指數均值都相接近，可見當下研究區各地的生態發展均衡且逐漸趨于穩定。

3.2.2 植被覆蓋度 植被覆蓋度相較于植被指數，是更能衡量地表植被覆蓋狀況的重要指標，實際研究中常將兩者結合來反映植被的生長情況。植被覆蓋度的計算值在0～1之間，根據閾值范圍和研究區植被覆蓋特點將其分為5個等級，具體分級和各覆蓋度面積比例統計情況見表3。

表3 研究區各級植被覆蓋度比例統計 %

據表3可知，1990年引黃濟青工程沿線面積的45.57%為中低植被覆蓋程度區域，低植被覆蓋程度區域次之，比例約36.91%，高植被覆蓋程度區域比例最小，僅有2.08%，總體植被覆蓋程度以低、中低覆蓋為主。2000年，研究區內低、中低覆蓋區域比例較1990年下降了約20%，中高植被覆蓋區域相對增多，但總體覆蓋程度仍以低、中低植被覆蓋度為主。2000—2010年間，區域內的植被發育較好，覆蓋度等級也發生了變化，由原先的以中低植被覆蓋為主轉變為以中高植被覆蓋為主，其中高植被區域比例約31.01%，處于遙遙領先地位，相比之下其他覆蓋等級相對均衡。現階段，研究區已完全以高植被覆蓋區域為主，且所占面積近總區域的一半。與前文中土地利用分類和植被指數的解譯結果一致，說明研究區內的植被不僅區域在增加，長勢也越來越茂盛。

3.3 地表溫度

利用大氣校正法得到的1990—2020年研究區內總體溫度分布情況(如圖4所示)。為了更清晰的研究其空間分布和變化趨勢，最終根據其值域范圍將地表溫度值劃分為了6個區間，分別為低于20，20～25，25～30，30～35，36～40 ℃，高于40 ℃。

圖4 研究區地表溫度空間分布

對各區域的地表溫度均值進行統計(見圖5)。從圖5地溫均值變化趨勢可知，4個時期研究區地表溫度均值分別為30.263，33.825，31.909，34.865 ℃，整體上呈現波動式變化。從濱州、東營、濰坊、青島4市區段的地溫變化來看，與總體的波動變化趨勢大體一致。其中，平均溫度最低的地區通常為青島區段。原因在于，青島臨靠黃海，受季風性和海洋性氣候共同作用。夏季溫濕多雨，但無酷暑，較其他地區溫度相對低；冬季風大溫和，且持續時間較長。相比之下，濰坊市、東營市、濱州市區段地理較為內陸，受海洋性氣候影響小，因此地表溫度表現的比青島區段高。

圖5 總研究區及各地級市區段地表溫度均值變化趨勢

結合土地分類結果，發現地表溫度值在空間分布上也有一定的差異性，水域地區溫度最低，大多處于10～25 ℃區間內。植被的蒸騰對溫度也有一定的降低作用，使得研究區內農田、林地、草地的地溫都在25～35 ℃之間。而其余的建筑用地和裸地區域受太陽直射，較為炎熱干燥，因此溫度大都在35 ℃以上。

3.4 土壤干濕度

利用溫度植被指數公式可以得到1990—2020年研究區總體土壤干濕度(即TVDI指數)的空間分布(圖6)。如圖6所示，1990，2000，2010，2020年研究區內的TVDI均值分別表現為0.758，0.694，0.568，0.433，由于TVDI指數的值與土壤濕度呈負相關，因此區域內的土壤呈現越來越濕潤的狀態。

圖6 研究區TVDI空間分布

圖7統計了1990—2020年研究區和各地級市區段土壤干濕度均值的變化趨勢，分析可知，濱州、東營、濰坊、青島市區段的TVDI均值變化呈現和總區域相同的下降趨勢，即土壤濕潤度越來越高。其中，青島區段的土壤濕潤度最高，濱州次之，濰坊、東營區段相對于沿海的青海和與黃河接壤的濱州，土壤都較為干旱。1990年，東營區段的TVDI指數高達0.866，土壤濕潤程度遠低于其他3區段，可能受到了水域、溫度、降雨等多方面因素的影響，由于土壤的干濕性對于植被的生長有很大的作用，因此這也在另一方面解釋了同年該區段內NDVI指數值不高的原因。

圖7 總研究區及各地級市區段TVDI均值變化趨勢

現階段，研究區內的土壤干濕程度已相對均衡，各地級市區段的平均TVDI指數值保持在0.39～0.45的區間內。

4 結 論

本文以引黃濟青工程輸水段沿線為研究區，采用遙感技術對生態因子的時空分布特征進行探究，得出結論如下。

(1) 研究區生態環境質量改善明顯，城鎮化水平發展迅速。30 a間，水域面積比例增長1.61%，綠地面積比例增長6.93%，裸地面積比例由29.52%縮減至2.61%，多用以建設基礎設施和發展生態、農業等。此外，各地級市區段土地結構變化與研究區總趨勢大致相同。

(2) 研究區內土壤濕潤度及植被覆蓋度顯著提升，且現階段發展趨于均衡。工程運行期間，土壤含水量增加，直接促進著區域內植被的繁茂生長及覆蓋程度的大幅提高。其中，青島區段植被生長狀況最好，土壤也最濕潤。

(3) 各生態因子間彼此不獨立。溫度的波動性、水域面積的增加、土壤的濕潤程度必然都影響著植被的生長，而植被的生長又促使著綠地面積的增加和裸土區域的利用，說明因子間存在一定程度上的相關性。

(4) 利用遙感技術研究較大區域的生態變化規律是一種可行、實用、便捷的方法，可以推廣應用于類似研究和實踐。