姜田亮， 張恒嘉， 馬國軍， 石媛媛， 王雅云

(甘肅農(nóng)業(yè)大學(xué) 水利水電工程學(xué)院， 甘肅 蘭州 730070)

1 研究背景

植被恢復(fù)是退化綠洲生態(tài)系統(tǒng)重建的關(guān)鍵，對于水資源供需矛盾突出，生態(tài)環(huán)境相對脆弱的干旱地區(qū)，植被生態(tài)需水量的準確計算對于指導(dǎo)該區(qū)域植被生態(tài)恢復(fù)與生態(tài)建設(shè)具有重要意義[1]。目前植被生態(tài)需水量計算的方法主要有直接計算法、間接計算法、水量平衡方程法以及基于遙感技術(shù)的計算方法[2]。對地下水依賴較強的干旱地區(qū)，天然植被的生態(tài)需水量可以通過潛水蒸發(fā)來估算，即間接計算法[3-5]，計算公式為：

(1)

式中:W為天然植被生態(tài)需水量，m3；Ei為第i種植被的潛水蒸發(fā)強度，mm；Ki是第i種植被影響系數(shù)；Si為第i種植被所占面積，hm2。對上述參數(shù)進行合理計算是準確量化天然植被生態(tài)需水量的關(guān)鍵。

潛水蒸發(fā)強度E的大小決定著影響植物生長的土壤水分狀況，土壤穩(wěn)定蒸發(fā)時，地表蒸發(fā)強度與土壤含水量都保持穩(wěn)定，三者在數(shù)值上相等。目前主要有兩種計算方法：一種是根據(jù)降雨量、溫度、輻射量等氣象因素建立模型，如Penman-Monteith模型、Hargreaves模型等[6]，但這類模型的計算需要長期觀測和積累大量氣象數(shù)據(jù)，很難應(yīng)用于部分工作基礎(chǔ)較差(難以通過試驗獲取所需數(shù)據(jù))且模型參數(shù)獲取困難的干旱地區(qū)；另一種是建立潛水蒸發(fā)量與埋深的關(guān)系模型，常見的關(guān)系模型大體可分為3類：第1類是潛水蒸發(fā)系數(shù)與潛水埋深呈單一的相關(guān)關(guān)系，如葉水庭指數(shù)型公式[7]，由于忽視了E與蒸發(fā)能力相關(guān)的規(guī)律，故有一定的局限性。第2類公式對此進行了改進，同時考慮了E與埋深及蒸發(fā)能力的關(guān)系，如沈立昌公式[8]。第3類是潛水蒸發(fā)系數(shù)與土壤輸水特性及土表蒸發(fā)具有相互關(guān)系，如雷志棟公式[9]。已有研究表明[10]，后兩類公式擬合埋深大于1 m的實測資料較好。由于3類計算公式具有不同的適用條件，本文基于MATLAB軟件利用去極值法和0.618法對7個常用公式所得結(jié)果進行優(yōu)選，并提出一種無實測數(shù)據(jù)情況下的潛水蒸發(fā)量估算方法。

植被系數(shù)K反映了一定區(qū)域內(nèi)植被對潛水蒸發(fā)量的影響，不同植被適宜生存的地下水埋深不同，利用植被地段的潛水蒸發(fā)量除以無植被地段的潛水蒸發(fā)量求得[11]。目前旱區(qū)植被系數(shù)選取，主要參照河西走廊玉門鎮(zhèn)有關(guān)的試驗結(jié)果[12]，即查表法。為克服查表法信息局限的缺點，本文基于已有試驗數(shù)據(jù)建立了潛水埋深與植被系數(shù)間的關(guān)系函數(shù)，為植被系數(shù)K的選取提供更為科學(xué)有效的方法。

植被覆蓋面積S是一個地區(qū)生態(tài)環(huán)境的綜合反映，現(xiàn)階段，植被覆蓋面積一般通過遙感影像人工解譯土地利用覆蓋圖，在此基礎(chǔ)上運用ARCGIS9.0軟件統(tǒng)計出各種天然植被的覆蓋情況[13]。研究植被面積的變化趨勢的常用方法有時間序列法、灰色模型法、生長曲線法、生產(chǎn)函數(shù)法、回歸分析法和馬爾科夫模型等[14]，其中灰色模型GM(1,1)不需要大量樣本就能得到較為準確的預(yù)測結(jié)果，可應(yīng)用于干旱地區(qū)植被面積的預(yù)測分析。

通過對公式中參數(shù)值的優(yōu)選及算法改進，有助于植被生態(tài)需水量的準確量化和間接計算法的推廣，為旱區(qū)水資源的分配及生態(tài)恢復(fù)政策的制定提供合理依據(jù)。

2 試驗地點及數(shù)據(jù)來源

2.1 民勤縣概況

民勤綠洲位于38°03′～39°28′N、101°49′～104°12′E，沿石羊河兩岸呈西南-東北向帶狀分布，西、北、東三面被巴丹吉林沙漠和騰格里沙漠包圍。氣候干燥，降水稀少，蒸發(fā)強烈，溫差較大，屬于典型的干旱地區(qū)，多年平均蒸發(fā)量2 644 mm[15]。歷史記載中民勤水草豐盛，草甸及沼澤植被是過去綠洲主要的原生植被[16]?，F(xiàn)如今地下水的嚴重超采致使地下水埋深急劇加大，土壤鹽漬化、沙化，水質(zhì)惡化日益嚴峻，進而導(dǎo)致植被覆蓋率降低,生態(tài)環(huán)境遭到嚴重破壞，給經(jīng)濟發(fā)展和區(qū)域內(nèi)人民的生活帶來了嚴重的威脅[17-18]。因此，民勤地區(qū)植被生態(tài)需水量準確計算對指導(dǎo)該地區(qū)水資源的優(yōu)化配置與生態(tài)建設(shè)具有重要的意義。

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

民勤地區(qū)各類植被面積數(shù)據(jù)以《民勤縣國民經(jīng)濟和社會發(fā)展統(tǒng)計資料匯編》[19]為準；對照組使用張奎俊[16]對石羊河流域下游進行研究所得潛水埋深-蒸發(fā)埋深數(shù)據(jù)；阿維里昂諾夫參數(shù)的選取依據(jù)中科院寒區(qū)旱區(qū)研究所和甘肅省水文二隊在典型干旱區(qū)玉門、張掖試驗所標定的值[20-21]；各植被最適地下水埋深數(shù)據(jù)以栗曉玲[22]和楊秀英等[23]的研究為依據(jù)。

3 研究方法

3.1 潛水蒸發(fā)量E的優(yōu)選

由于現(xiàn)有潛水蒸發(fā)模型適用條件不同，若利用常規(guī)方法對潛水蒸發(fā)量計算，則對不同埋深下的潛水蒸發(fā)量反映不夠理想。因此，利用MATLAB將各類公式(表1中公式1～7)編入同一個程序，分別對不同埋深的潛水蒸發(fā)點進行擬合，排除不滿足相關(guān)度的模型公式后，結(jié)合去極值法和0.618法進行結(jié)果優(yōu)選，編程思路如下：

由于民勤縣極限地下水位埋深hmax為4.5m[18]，因此h以[1,4.5]為區(qū)間、0.01為步長增長，代入各公式得到350組不同的潛水蒸發(fā)量，構(gòu)成矩陣E：

去掉每行的最大值、最小值得到新的矩陣E′：

用下式求得優(yōu)化結(jié)果E優(yōu)化:

(2)

(i=1,2,…,350)

表1 潛水蒸發(fā)公式

注：a,b,c為經(jīng)驗參數(shù)，即擬合系數(shù)；E0為大氣蒸發(fā)量,mm/d；h為地下水埋深,m。

針對部分研究地區(qū)試驗條件有限，缺乏實測潛水蒸發(fā)量-埋深數(shù)據(jù)的現(xiàn)狀，可以根據(jù)研究區(qū)的土壤和氣候特點，使用相似地區(qū)阿維里昂諾夫經(jīng)驗參數(shù)模擬出數(shù)據(jù)點進行計算。阿維里昂諾夫模型形式簡單(公式8)，易于推廣，現(xiàn)已被廣泛應(yīng)用于實際生產(chǎn)，積累了大量的研究資料。選取3組阿維里昂諾夫經(jīng)驗參數(shù)，以3組不同的潛水埋深(1、2、3、4 m);(1.35、2.65、3.35、4.35 m);(1.75、2.25、3.75、3.85 m)生成模擬蒸發(fā)點進行擬合。

3.2 植被影響系數(shù)K的計算

K值的常規(guī)解法是查表法[22-23]，如表2所示。但當潛水埋深取值不在表中時，K只能通過專業(yè)人員估值得到，降低了結(jié)果的準確度。利用MATLAB中的二階Gaussian函數(shù)對數(shù)據(jù)進行準確擬合，結(jié)果如公式(3)和圖1。

(3)

式中:h為潛水埋深,m。經(jīng)檢驗，R為0.9957，擬合度較高。

圖1 Ki值擬合曲線

表2 干旱區(qū)潛水埋深與植被影響系數(shù)表

3.3 植被面積S的研究預(yù)測

利用灰度模型GM(1,1)分別對喬木林、灌木林、疏林地、高蓋度草地、中蓋度草地2005-2014年的面積數(shù)據(jù)進行率定，通過后驗差檢驗，得到2015年各植被面積，與記載資料[19]進行對比驗證，進而對2020年植被面積進行預(yù)測。灰度模型具體思路如下[29]：

(1)各植被面積原始序列數(shù)據(jù)為x(0)(k)={x(0)(1),x(0)(2),…,x(0)(n)}，并對其作累加處理得：x(1)(k)={x(1)(1),x(1)(2),…,x(1)(n)}，k=1,2,…,n。

(2)灰度模型GM(1,1)的建立。即：

(4)

式中:μ為內(nèi)生控制灰數(shù);α為發(fā)展灰數(shù)。

(3)向量Y和矩陣B的構(gòu)建：

(5)

Y=[x(0)(2),x(0)(3),…,x(0)(n)]T

(6)

(4)利用最小二乘法估計模型參數(shù)可得：

a=BT(BTB)-1Y

(7)

(5)預(yù)測模型求解：

(8)

k=0，1，2，…，n

(6)利用累減還原法得到灰色預(yù)測模型：

(9)

(7)后驗差法計算殘差，得：

(10)

k=1,2,…,n

(11)

(12)

后驗差比為:

C=S2/S1

(13)

精度檢驗等級如表3所示。

表3 精度檢驗等級參照表

4 結(jié)果分析

4.1 潛水蒸發(fā)量

在充分考慮不同植被蓋度的蒸騰量、不同水質(zhì)以及地下水埋深與植被蒸騰量的相關(guān)關(guān)系的基礎(chǔ)上，對表4中阿維里昂諾夫初始參數(shù)生成的模擬點進行擬合，結(jié)果如圖2～5。各公式擬合度均大于R0.05，線性關(guān)系合理，與3種經(jīng)典公式相比，優(yōu)選曲線擬合程度更好。與張奎俊[16]在民勤所做潛水蒸發(fā)試驗結(jié)果進行對比，相關(guān)系數(shù)R分別為0.918，0.917，0.907，0.926，P值分別為0.912，0.849，0.900，0.911，0.934，其中優(yōu)化結(jié)果的相關(guān)系數(shù)最大，且無差異的顯著性水平最高，如表5所列，說明了優(yōu)化的合理性。

表4 阿維里昂諾夫初始參數(shù)

圖2 沈立昌公式擬合結(jié)果 圖3 葉水庭公式擬合結(jié)果

圖4 雷志棟公式擬合結(jié)果 圖5 優(yōu)選結(jié)果

表5 各公式擬合結(jié)果與張奎俊模擬點相關(guān)性分析及差異性檢驗

4.2 植被系數(shù)

據(jù)粟曉玲[22]和楊秀英等[23]的研究，民勤地區(qū)喬木林的最適地下水埋深為1～4.5 m，灌木林3～4 m，疏林地3.5～4.5 m，高覆蓋度草3～4.5 m，中覆蓋度草4～4.5 m。選取中間值作為計算埋深代入公式(2)，得到各類植被影響系數(shù)及潛水蒸發(fā)量，結(jié)果如表6所列。

表6 各類植被不同地下水埋深植被影響系數(shù)及潛水蒸發(fā)量

4.3 植被面積

利用民勤縣2005年至2014年的數(shù)據(jù)對模型進行率定，依據(jù)公式(4)和(7)計算得出 ，代入式(8)則分別得到喬木林、灌木林、疏林地、高覆蓋度草及中覆蓋度草的面積累計預(yù)測公式，得到5個單一灰色GM(1,1)模型，結(jié)果皆通過后驗差法進行檢驗?；贛atlab編程并代入面積數(shù)據(jù)運行，結(jié)果如表7。

在此基礎(chǔ)上計算，到2015年，喬木林、灌木林、疏林地面積分別為1 475、7 181、13 025 hm2，高覆蓋度草地和中覆蓋度草地面積分別為1 571、6 873 hm2?！睹袂诳h2015年國民經(jīng)濟和社會發(fā)展統(tǒng)計資料匯編》統(tǒng)計，喬木林、灌木林、疏林地面積分別為1 498、7 943、12 901 hm2，高覆蓋度草地、中覆蓋度草地面積分別為1 600、6 954 hm2。兩組數(shù)據(jù)p值為0.9611，無差異。2015年結(jié)果如圖6～10所示。代入公式 (1)，得到2015年民勤縣生態(tài)需水量為3 726.8×104m3，其中喬木林需水量為299.0×104m3，灌木林1 208.9×104m3，疏林地需水量為1 310.2×104m3，高覆蓋度草為217.3×104m3，中覆蓋度草地為691.4×104m3。

表7 各植被類型面積計算公式及后驗差檢驗

注：k=2～10，代表2006-2014年的時間序列，x(0)(1)為2005年的植被面積。

圖6 喬木林面積增長趨勢 圖7 灌木林面積增長趨勢 圖8 疏林地面積增長趨勢

因此可以預(yù)測，到2020年，民勤縣喬木林面積將達到1 623 hm2，灌木林7 899 hm2,疏林地14 328 hm2,高覆蓋度草1 728 hm2,中覆蓋度草7 141 hm2。屆時民勤縣天然植被生態(tài)需水量為4 060.8×104m3，其中喬木林329.0×104m3，灌木林地1 329.8×104m3，疏林地1 441.3×104m3，高覆蓋度草地239.0×104m3，中覆蓋度草地721.7×104m3。

圖9 高覆蓋度草面積增長趨勢 圖10 中覆蓋度草面積增長趨勢

5 討 論

本研究探求潛水蒸發(fā)量E是對眾多潛水蒸發(fā)公式所得結(jié)果的優(yōu)選，無具體的模型公式，所得潛水蒸發(fā)量-埋深曲線是由大量優(yōu)選點組成的折線，需要借助計算軟件(如MATLAB)進行編程求解。由圖2～5可見，當埋深較淺時，對潛水蒸發(fā)量取值影響較大的是第一類公式；當埋深大于3.5 m，由更適宜計算該埋深條件的二、三類公式?jīng)Q定潛水蒸發(fā)量降低速率，因此圖5中優(yōu)化的潛水蒸發(fā)量-埋深曲線更具有代表性。

由圖1見，植被系數(shù)K的減小速率隨埋深增先逐漸變小，后趨于穩(wěn)定，最后逐漸變大。已有文獻[30-31]發(fā)現(xiàn)，不同植被影響系數(shù)的大小是一個隨潛水埋深變化的函數(shù)，其值隨地下水埋深的增大而減少。當?shù)叵滤裆钚∮? m時，減小速率隨埋深的加大逐漸變小，當埋深在2～3 m時，K值趨于穩(wěn)定；當埋深大于3 m時，減小速率逐漸上升,可以用二階Gaussian函數(shù)精確的擬合。朱艷紅[32]研究表明,干旱地區(qū)地下水埋深小于2 m時，土壤含水率減小速率隨埋深加大逐漸變小，當潛水埋深大于2 m時，土壤含水率趨于穩(wěn)定。上述變化規(guī)律與圖1中K值變化趨勢一致。因此可得出，在潛水埋深小于3 m的范圍內(nèi)，土壤水蒸發(fā)是影響旱區(qū)植被生態(tài)需水量的主要因素，此時K值的變化速率與土壤含水率變化速率一致。當潛水埋深大于3 m時，土壤水蒸發(fā)對植被生態(tài)需水量的影響逐漸變小，而地表植被由蒸騰作用更為強烈的喬木、灌木變?yōu)檎趄v作用較弱草地，導(dǎo)致K值迅速變小。

Penman-Monteith模型是一種根據(jù)降雨量、溫度、輻射量等氣象因素建立的潛水蒸發(fā)模型，郝博等[33]利用該模型計算得出2015年民勤縣總生態(tài)需水量4 118×104m3，其中喬木林278×104m3，灌木林地1 194×104m3，疏林地1 431×104m3，高覆蓋度草地203×104m3，中覆蓋度草地702×104m3。與本研究所得結(jié)果相比，各類植被生態(tài)需水量相關(guān)度為0.9739，大于R0.05，p值為0.962，兩者無差異。總生態(tài)需水量較本研究結(jié)果偏高，究其原因是該方法計算結(jié)果中的天然林地包括除以上3種林地外的其他林地。

6 結(jié) 論

(1)通過優(yōu)化植被生態(tài)需水量間接法計算公式中每個參數(shù)，更加準確的計算和預(yù)測了民勤縣的天然植被生態(tài)需水量，該方法也為我國部分環(huán)境惡劣、試驗條件艱苦而難以取得準確觀測數(shù)據(jù)的干旱地區(qū)生態(tài)需水量的估值提供了科學(xué)的參照。其中潛水蒸發(fā)量的算法以去極值法和0.618優(yōu)選法為基礎(chǔ)，植被影響系數(shù)的核心是一元函數(shù)的擬合，植被面積的計算及預(yù)測運用灰色模型理論，整個方法所使用的理論簡單有效、可操作性強，是一種容易推廣的計算方法。

(2)預(yù)測得到2020年民勤縣的生態(tài)需水量將達到4 060.7×104m3，即2015年至2020年以每年66.78×104m3的速度增加。為了恢復(fù)綠洲植被面積，實現(xiàn)民勤綠洲的社會經(jīng)濟可持續(xù)發(fā)展的戰(zhàn)略，必須先滿足綠洲天然植被的生態(tài)需水。地下水的持續(xù)超采是民勤綠洲生態(tài)惡化的主要原因，因此石羊河的地表徑流是民勤最主要的供水源。所以，在合理規(guī)劃生活、生產(chǎn)用水的前提下，增加進入民勤綠洲的地表水量，是解決水資源危機、緩解地下水持續(xù)超采、保證民勤綠洲安全的最重要對策，亦是不容延誤的當務(wù)之急。

(3)本研究對民勤縣潛水蒸發(fā)量進行估算時，由于未使用實測潛水蒸發(fā)量-埋深點，因此估算結(jié)果很大程度上取決于阿維里昂諾夫公式初始參數(shù)的選擇，其中初始參數(shù)a與優(yōu)化結(jié)果呈正相關(guān)，b呈負相關(guān)。為了更加準確的估算天然植被生態(tài)需水量，消除估算所得潛水蒸發(fā)量與實際值間偏差，今后應(yīng)考慮觀測和確定不同植被類型在不同地下水位所對應(yīng)的潛水蒸發(fā)參數(shù)，準確掌握植被蒸騰與潛水位之間的關(guān)系，以提高生態(tài)需水的計算精度。