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

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

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

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

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

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

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

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

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

研究區(qū)范圍廣闊、地形復(fù)雜，適合衛(wèi)星遙感進行大面積、多角度、高時效的生態(tài)環(huán)境監(jiān)測?？紤]到研究周期和影像精度的綜合要求，最終選用Landsat系列影像作為初始數(shù)據(jù)。

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

表1 研究區(qū)遙感影像基本信息

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

2 研究方法

2.1 技術(shù)路線

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

2.2 土地利用/覆蓋的量化

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

2.3 植被的量化

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

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

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

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

2.4 地表溫度的量化

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

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

(3)

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

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

ε=0.004PV+0.986

(5)

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

式中：K1，K2為定標(biāo)系數(shù)，可通過影像的元數(shù)據(jù)獲取。

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

2.5 土壤干濕度的量化

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

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

Tsmax=a1+b1×NDVI

(8)

Tsmin=a2+b2×NDVI

(9)

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

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

3 結(jié)果與分析

3.1 土地利用/覆蓋

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

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

表2 研究區(qū)土地利用分類面積及比例統(tǒng)計

由圖1可知，在時間維度上，1990年研究區(qū)內(nèi)土地結(jié)構(gòu)以綠地和裸地為主，建筑用地較少，水域地區(qū)更為稀缺，說明工程沿線土地資源豐富，但城鎮(zhèn)化水平較低，城市基礎(chǔ)設(shè)施建設(shè)較落后。1990—2010年期間，綠地面積變化幅度不大，相比之下建筑用地面積增勢明顯，同時裸地區(qū)域也在以肉眼可見的速度減少。2020年區(qū)域內(nèi)地類主要以綠地和建筑用地為主，裸地區(qū)域被大面積開墾利用，可見當(dāng)?shù)爻鞘谢ㄔO(shè)發(fā)展迅速。

在空間維度上，按研究區(qū)內(nèi)4個地級市區(qū)段劃分并對各區(qū)段地類進行統(tǒng)計，進一步分析得到濱州和青島區(qū)段的水域面積較其他區(qū)段相對較大?？赡苁芷涞乩砦恢玫挠绊懀瑸I州地處引黃濟青工程黃河入水口，而青島東、南瀕臨黃海，帶來的河湖水源較充足。同時，水源和地理的優(yōu)勢也使得青島區(qū)段的耕地、林地、草地發(fā)展的較好，其綠地面積比例在4個區(qū)段中一直都居于前列。建筑用地區(qū)域在4個區(qū)段中區(qū)別并不大，說明研究區(qū)內(nèi)各地區(qū)城鎮(zhèn)化發(fā)展相對較為均衡。

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

3.2 植被指數(shù)及覆蓋度

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

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

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

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

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

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

表3 研究區(qū)各級植被覆蓋度比例統(tǒng)計 %

據(jù)表3可知，1990年引黃濟青工程沿線面積的45.57%為中低植被覆蓋程度區(qū)域，低植被覆蓋程度區(qū)域次之，比例約36.91%，高植被覆蓋程度區(qū)域比例最小，僅有2.08%，總體植被覆蓋程度以低、中低覆蓋為主。2000年，研究區(qū)內(nèi)低、中低覆蓋區(qū)域比例較1990年下降了約20%，中高植被覆蓋區(qū)域相對增多，但總體覆蓋程度仍以低、中低植被覆蓋度為主。2000—2010年間，區(qū)域內(nèi)的植被發(fā)育較好，覆蓋度等級也發(fā)生了變化，由原先的以中低植被覆蓋為主轉(zhuǎn)變?yōu)橐灾懈咧脖桓采w為主，其中高植被區(qū)域比例約31.01%，處于遙遙領(lǐng)先地位，相比之下其他覆蓋等級相對均衡。現(xiàn)階段，研究區(qū)已完全以高植被覆蓋區(qū)域為主，且所占面積近總區(qū)域的一半。與前文中土地利用分類和植被指數(shù)的解譯結(jié)果一致，說明研究區(qū)內(nèi)的植被不僅區(qū)域在增加，長勢也越來越茂盛。

3.3 地表溫度

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

圖4 研究區(qū)地表溫度空間分布

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

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

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

3.4 土壤干濕度

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

圖6 研究區(qū)TVDI空間分布

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

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

現(xiàn)階段，研究區(qū)內(nèi)的土壤干濕程度已相對均衡，各地級市區(qū)段的平均TVDI指數(shù)值保持在0.39～0.45的區(qū)間內(nèi)。

4 結(jié) 論

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

(1) 研究區(qū)生態(tài)環(huán)境質(zhì)量改善明顯，城鎮(zhèn)化水平發(fā)展迅速。30 a間，水域面積比例增長1.61%，綠地面積比例增長6.93%，裸地面積比例由29.52%縮減至2.61%，多用以建設(shè)基礎(chǔ)設(shè)施和發(fā)展生態(tài)、農(nóng)業(yè)等。此外，各地級市區(qū)段土地結(jié)構(gòu)變化與研究區(qū)總趨勢大致相同。

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

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

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