遲凱歌，龐 博，石樹蘭，崔黎壯

(1. 北京師范大學(xué)水科學(xué)研究院，北京 100875；2. 城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875)

0 引 言

植被作為陸地生態(tài)系統(tǒng)的主要組分，是生態(tài)系統(tǒng)存在的基礎(chǔ)，也是聯(lián)結(jié)土壤、大氣和水分的自然“紐帶”，它在陸地表面的能量交換過程、生物地球化學(xué)循環(huán)過程和水文循環(huán)過程中扮演著重要的角色，在全球變化研究中起著“指示器”的作用[1-3]。NDVI(the Normalized Vegetation Index)是研究植被變化的重要工具，能夠較好地反映出區(qū)域下墊面植被動(dòng)態(tài)變化的規(guī)律。

近年來，國內(nèi)外在NDVI預(yù)測模型方面開展了大量的研究。例如，IWASAKI利用GSMaP(降水?dāng)?shù)據(jù)種類)降水量和JRA-25/JCDAS(溫度數(shù)據(jù)種類)溫度對蒙古草原進(jìn)行NDVI預(yù)測，建立了逐步多元線性回歸分析預(yù)測模型[4]；阿西穆等人基于Markov模型對渭干河和庫車河流域的綠洲地區(qū)的植被覆蓋進(jìn)行動(dòng)態(tài)變化預(yù)測研究[5]；林楠等人以吉林東部為研究區(qū)在進(jìn)行植被覆蓋動(dòng)態(tài)變化分析的基礎(chǔ)上，利用馬爾柯夫模型和灰色系統(tǒng)G(1，1)模型對研究區(qū)植被覆蓋進(jìn)行動(dòng)態(tài)模擬和預(yù)測，分析其演變特點(diǎn)及規(guī)律[6]。WANG等人采用多元線性回歸方法對我國南方丘陵地區(qū)的降水和溫度進(jìn)行NDVI模擬[7]。雅魯藏布江流域地形特征復(fù)雜，氣候變化敏感，生態(tài)系統(tǒng)脆弱，如何綜合考慮氣候要素對流域植被覆蓋的影響，構(gòu)建機(jī)制清晰，精度可靠的植被覆蓋預(yù)測模型成為生態(tài)水文領(lǐng)域的挑戰(zhàn)[8-10]。隨著機(jī)器學(xué)習(xí)技術(shù)的發(fā)展，人工神經(jīng)網(wǎng)絡(luò)由于其強(qiáng)大的非線性模擬能力在預(yù)測模型領(lǐng)域得到了廣泛的應(yīng)用。因此，本文綜合采用主成分分析、偏相關(guān)分析和人工神經(jīng)網(wǎng)絡(luò)模型，嘗試構(gòu)建雅魯藏布江流域的NDVI預(yù)測模型，探索為氣候地形復(fù)雜區(qū)域的植被覆蓋研究提供一種新的方法。

本研究采用雅魯藏布江流域及其周圍地區(qū)30個(gè)地面站點(diǎn)數(shù)據(jù)和MODIS數(shù)據(jù)，計(jì)算得到流域內(nèi)NDVI的數(shù)值變化的時(shí)空規(guī)律，對流域內(nèi)植被動(dòng)態(tài)變化進(jìn)行了分析，并進(jìn)一步構(gòu)建了NDVI預(yù)測模型，對進(jìn)一步探求雅魯藏布江流域下墊面的驅(qū)動(dòng)規(guī)律，積極應(yīng)對氣候變化對下墊面帶來的影響具有十分重要的意義。

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

1.1 研究區(qū)概況

雅魯藏布江流域是青藏高原上最大的河流，也是最重要的國際河流之一，發(fā)源于西藏自治區(qū)南部，喜馬拉雅山北麓的杰瑪央宗冰川[11,12]。雅魯藏布江流域橫貫青藏高原南部，河流總長為2 057 km，流域面積為24.2 萬km2[13,14]。雅魯藏布江流域(圖1)是世界上最高的河流之一，平均海拔在4 500 m以上，流域地勢西高東低，其上較大支流有拉薩河、帕隆藏布、年楚河、多雄藏布和尼洋河等[15]。雅魯藏布江流域的氣候條件主要受高原地理位置和地勢特點(diǎn)支配[16]，上游地區(qū)屬高原寒溫帶氣候，年平均降水量不足300 mm。中游地區(qū)屬高原溫帶氣候，年平均降水量為300～600 mm，河谷地帶年平均氣溫在4.7～8.6 ℃之間[17,18]。下游地區(qū)屬熱帶亞熱帶氣候，下游地區(qū)的巴昔卡附近年均降雨量在4 000 mm以上。

受特殊地質(zhì)地貌和水熱因素的影響，流域自然景觀多樣。流域下游水汽條件良好,主要植被類型有熱帶低山半常綠雨林、亞熱帶常綠闊葉林、亞熱帶山地半常綠闊葉林和常綠針葉林；中游地區(qū)植被類型為山地與河谷灌叢草原，普遍分布著灌叢草原植被,草本植物為中溫型禾草；上游地區(qū)主要植被類型為高寒草原、高寒草甸、高寒灌叢以及高寒墊狀植物。

1.2 數(shù)據(jù)來源與預(yù)處理

本文利用2000-2015年雅魯藏布江的NDVI資料，研究時(shí)段內(nèi)的雅魯藏布江的植被覆蓋變化。在神經(jīng)網(wǎng)絡(luò)的模擬模式中，輸入為實(shí)測的降雨和氣溫，輸出為模擬的NDVI指數(shù)，NDVI數(shù)據(jù)來自MODIS遙感衛(wèi)星數(shù)據(jù)，已經(jīng)過云處理。雅魯藏布江降雨和氣溫?cái)?shù)據(jù)收集自國家地球系統(tǒng)科學(xué)數(shù)據(jù)平臺(tái)，研究數(shù)據(jù)跨度時(shí)間是2000年2月至2015年12月。其中NDVI數(shù)據(jù)和氣象數(shù)據(jù)已經(jīng)分別應(yīng)用于雅魯藏布江的植被變化研究[19]和極端氣候事件研究[20]，數(shù)據(jù)在下載前通過了均一性檢驗(yàn)以及質(zhì)量控制，具有較高的可靠性。本文用到的雨量站總共有30個(gè)，并且根據(jù)雅魯藏布江地區(qū)自然水系將區(qū)域分為5個(gè)子區(qū)域，一區(qū)的雨量站有普蘭、改則、聶拉木、定日、拉孜和申扎；二區(qū)的雨量站有日喀則、帕里、江孜、南木林、尼木和浪卡子；三區(qū)的雨量站有當(dāng)雄、拉薩、墨竹工卡和那曲；四區(qū)的雨量站有嘉黎、波密、八宿、洛隆和比如；五區(qū)的雨量站有林芝、加查、米林、澤當(dāng)、貢嘎。子流域分區(qū)情況見表1。

表1 雅魯藏布江子流域分區(qū)表Tab.1 Yarlung Zangbo River Basin Division

2 研究方法

2.1 主成分分析

主成分分析(Principal Component Analysis, PCA)最早是由Pearson于1901年提出的一種通過對協(xié)方差矩陣進(jìn)行分析，主要目的是在降低數(shù)據(jù)維數(shù)的條件下保持?jǐn)?shù)據(jù)集對方差貢獻(xiàn)最大的一種統(tǒng)計(jì)分析方法。其原理為通過降維的方法將原來的多個(gè)指標(biāo)轉(zhuǎn)化為一個(gè)或幾個(gè)綜合指標(biāo)，即主成分(PCS)，而各主成分之間互不相關(guān)。即

(1)

式中：x1，x2，…，xp，為原始變量；Z1，Z2，…，Zm(m≤p)為新變量指標(biāo)；lmp為先變量權(quán)重。

2.2 偏相關(guān)分析

偏相關(guān)分析(Partial Correlation Analysis，PAR)是在多要素構(gòu)成的自然系統(tǒng)中，對其他變量的影響進(jìn)行控制的條件下，衡量多個(gè)變量中某兩個(gè)變量之間的線性相關(guān)關(guān)系的密切程度的分析方法，以偏相關(guān)系數(shù)來度量偏相關(guān)程度。偏相關(guān)計(jì)算公式的通式為(r為兩個(gè)變量之間的相關(guān)系數(shù))：

(2)

一階偏相關(guān)系數(shù)計(jì)算公式為：

(3)

即消除第3變量影響，分析兩變量的相關(guān)關(guān)系，二階偏相關(guān)系數(shù)計(jì)算公式為：

(4)

即消除第 3變量和第 4 變量的影響，分析兩變量的相關(guān)關(guān)系。其他多階偏相關(guān)系數(shù)計(jì)算公式都可以通過N階偏相關(guān)系數(shù)計(jì)算通式推導(dǎo)出來，對求得的相關(guān)系數(shù)，可以用顯著性檢驗(yàn)方法進(jìn)行顯著性檢驗(yàn)。對相關(guān)系數(shù)顯著性檢驗(yàn)方法多選取F檢驗(yàn)進(jìn)行檢測，F(xiàn)檢驗(yàn)又稱方差齊性檢驗(yàn)，主要通過對比2個(gè)變量數(shù)據(jù)的方差，以確定其吻合度是否有顯著性差異，當(dāng)數(shù)據(jù)變量個(gè)數(shù)充分大時(shí)，統(tǒng)計(jì)量近似服從F分布。

2.3 人工神經(jīng)網(wǎng)絡(luò)模型

人工神經(jīng)網(wǎng)絡(luò)是借鑒大腦和神經(jīng)系統(tǒng)存儲(chǔ)和處理信息的特征抽象出來的一種數(shù)學(xué)模型。前饋BP(Back Propagation )神經(jīng)網(wǎng)絡(luò)是應(yīng)用廣泛的靜態(tài)網(wǎng)絡(luò)。BP網(wǎng)絡(luò)以其良好的泛函逼近能力以及自組織、自適應(yīng)性能和容錯(cuò)性等優(yōu)點(diǎn)為水文系統(tǒng)的建模、識(shí)別和預(yù)測提供了一種有效的研究手段[21]。BP 神經(jīng)網(wǎng)絡(luò)是一種包含輸入層、隱含層和輸出層的前向網(wǎng)絡(luò)結(jié)構(gòu)，它的基本特點(diǎn)是信號(hào)經(jīng)輸入層、隱含層逐步處理后向前傳遞直達(dá)輸出層，神經(jīng)網(wǎng)絡(luò)每一層神經(jīng)元的狀態(tài)只由上一層決定，如果輸出層得到的輸出結(jié)果與期望輸出不在誤差允許范圍內(nèi)，則信號(hào)逆向反饋，調(diào)整網(wǎng)絡(luò)的權(quán)值和各神經(jīng)元的閾值，使得網(wǎng)絡(luò)輸出與期望的輸出無限靠近，直至誤差允許的范圍內(nèi)[21]，具有以非線性適應(yīng)性為特征的信息處理能力、十分優(yōu)良的非線性逼近能力以及較好的容錯(cuò)性，因而能夠?qū)Υ笠?guī)模數(shù)據(jù)進(jìn)行并行處理，且自組織，自學(xué)習(xí)與自適應(yīng)能力突出，被廣泛應(yīng)用于預(yù)測、分類、模式識(shí)別和聚類等領(lǐng)域[22]。

BP模型分為訓(xùn)練與檢驗(yàn)兩個(gè)部分，模型由率定數(shù)據(jù)集進(jìn)行訓(xùn)練，當(dāng)達(dá)到最小誤差限制或者最大迭代次數(shù)時(shí)，模型訓(xùn)練結(jié)束。模型的數(shù)學(xué)表達(dá)由下式所示：

(5)

式中：xi為節(jié)點(diǎn)i的輸入值；Ok為節(jié)點(diǎn)k的輸出值；g1為隱含層的激活函數(shù)；g2為輸出層的激活函數(shù)；M和N分別是輸入層和隱含層中神經(jīng)元的個(gè)數(shù)；wj0為隱含層中第j個(gè)神經(jīng)元的偏差；wk0為輸出層中第k個(gè)神經(jīng)元的偏差；wji為輸入節(jié)點(diǎn)i與隱含節(jié)點(diǎn)j之間的權(quán)重；wkj為隱含節(jié)點(diǎn)j和輸出節(jié)點(diǎn)k之間的權(quán)重。

模型計(jì)算在MATLAB(2016)環(huán)境下進(jìn)行。ANN模型采用MATLAB中BP神經(jīng)網(wǎng)絡(luò)工具箱，隱含層和輸出層的傳遞函數(shù)分別使用Sigmoid函數(shù)和線性函數(shù)，并且為了避免過度擬合采用了早停函數(shù)earlystopping有效避免過度迭代。其中，2000-2009年數(shù)據(jù)進(jìn)行模型率定，2010-2015年數(shù)據(jù)為模型驗(yàn)證。

2）修改的后64卦卦序（圖8），揭示了64卦序數(shù)／陰陽對稱規(guī)律，這也許是一個(gè)史無前例的、劃時(shí)代的發(fā)現(xiàn)！呈現(xiàn)出既陰陽對稱又?jǐn)?shù)量平衡，陰中有陽，陽中有陰，陰陽相抱，自然天成。這豈不是精美的“陰陽／數(shù)量太極圖”嗎？64卦若采用比干八宮卦序，那么乾為首卦就要改變了，即坤為1，乾為2。

本研究的主要計(jì)算過程如下：首先利用2000-2015年的MODIS數(shù)據(jù)，利用歸一化植被指數(shù)的方法得到16年間NDVI逐月的數(shù)值大小，以及流域NDVI的多年平均值。其次，使用線性擬合的方法得到流域內(nèi)5個(gè)分區(qū)NDVI的增長速率，并計(jì)算NDVI在各個(gè)分區(qū)上的Hurst指數(shù)，以分析NDVI變化的持續(xù)性。最后，分別采用主成分分析法和偏相關(guān)分析法提取得到影響NDVI的主要影響因子，作為輸入因子輸入人工神經(jīng)網(wǎng)絡(luò)模型。模型基于MATLAB計(jì)算平臺(tái)，其中未采用因子篩選的模型簡稱為ANN，采用主成分分析進(jìn)行因子篩選的模型簡稱為ANN-PCA, 采用偏相關(guān)分析進(jìn)行因子篩選的模型簡稱為ANN-PAR，最后得到NDVI的預(yù)測模型。

3 模型評(píng)價(jià)指標(biāo)

本文中選取4種模型評(píng)價(jià)指標(biāo)來評(píng)價(jià)各個(gè)模型的表現(xiàn)能力：納西效率系數(shù)(NASH)、均方根誤差(RMSE)、平均絕對誤差(MAE)和相關(guān)系數(shù)(R)，其計(jì)算公式如下所示：

(6)

(7)

(8)

(9)

4 研究結(jié)果與分析

4.1 NDVI的時(shí)空變化趨勢

雅魯藏布江流域2000-2015的NDVI年均空間分布特征如圖2所示，如圖2所示，流域NDVI呈從西向東，從北向南的增加趨勢。受印度洋暖濕氣流影響，區(qū)域5的NDVI顯著高于其他區(qū)域。

其中，區(qū)域1到區(qū)域5的NDVI多年平均值分別為0.13、0.18、0.24、0.23、0.31。5個(gè)區(qū)域的植被狀況穩(wěn)定，呈現(xiàn)不顯著的增長趨勢。以區(qū)域5為例，對子流域的NDVI的多年平均數(shù)值進(jìn)行線性擬合，從圖中可以看到，在2000-2015年期間，NDVI區(qū)域5上呈現(xiàn)出增加的趨勢，其傾向率為0.002/(10 a)，變化幅度較小，這與全球氣候變化背景下，中高緯度降水增加的變化趨勢相一致。

本文對雅魯藏布江地區(qū)NDVI的變化的持續(xù)性也做出了研究。各個(gè)區(qū)域Hurst指數(shù)計(jì)算結(jié)果如表2所示。如圖，區(qū)域1到區(qū)域5的Hurst指數(shù)呈現(xiàn)上升趨勢。其中區(qū)域1呈弱反持續(xù)性，當(dāng)前該區(qū)域的趨勢不是很明顯比較穩(wěn)定，在未來植被有可能會(huì)稍微減少。區(qū)域2到5都呈現(xiàn)持續(xù)性，表明未來NDVI將持續(xù)增加，而且從上游到下游的持續(xù)性越來越強(qiáng)。

表2 各個(gè)子流域區(qū)域Hurst指數(shù)計(jì)算結(jié)果Tab.2 Results of Hurst index of each sub-watershed area

各個(gè)子流域的NDVI的MK突變檢驗(yàn)結(jié)果如圖4所示。如圖所示，區(qū)域突變不明顯，NDVI變化相對穩(wěn)定。

4.2 影響因子分析

分別采用主成分分析法和偏相關(guān)分析法對影響NDVI變化的主要?dú)庀笠蜃舆M(jìn)行分析。其中，主成分分析方法本文選取累計(jì)方差貢獻(xiàn)率高于90%的主成分作為主要影響因子，結(jié)果如表3所示。偏相關(guān)分析的結(jié)果如表3所示，且本文選擇偏相關(guān)系數(shù)大于0.3的變量作為主要影響因子(結(jié)果加粗的為偏相關(guān)分析中篩選出的影響因子)。

從表3和表4可以看出，氣溫和降水對各區(qū)域NDVI的影響分別隨著時(shí)間推移逐步衰減。主成分結(jié)果統(tǒng)計(jì)表明，各區(qū)域前3個(gè)月的氣溫和降水會(huì)對NDVI有主要影響，而偏相關(guān)分析中影響因子的范圍縮小主要集中在前兩個(gè)月的降雨和氣溫。兩種方法的結(jié)果表明，總體上同時(shí)期氣溫對NDVI的影響要略高于對降水的影響，利用兩種方法對因子進(jìn)行篩選，并將篩選出的因子作為輸入輸進(jìn)人工神經(jīng)網(wǎng)絡(luò)模型進(jìn)行NDVI模擬可以提升精度。偏相關(guān)分析與主成分分析的結(jié)果與相關(guān)研究的結(jié)論較為一致。例如郭兵[23]等提出雅魯藏布江流域植被覆蓋度受降水和氣溫的影響較大，但是作用到植被覆蓋具有一定的滯后性。本文采的影響因子相關(guān)分析結(jié)果與其結(jié)論相似。

表3 主成分累計(jì)貢獻(xiàn)率及所選主成分 %

表4 各個(gè)成分偏相關(guān)計(jì)算結(jié)果Tab.4 Partial correlation calculation results

4.3 模型預(yù)報(bào)結(jié)果

輸入數(shù)據(jù)的選擇是人工神經(jīng)網(wǎng)絡(luò)模型構(gòu)建的關(guān)鍵問題之一[22]。如何有效地利用數(shù)據(jù)信息，同時(shí)避免由于過多的輸入數(shù)據(jù)導(dǎo)致網(wǎng)絡(luò)結(jié)構(gòu)過于復(fù)雜，出現(xiàn)過擬合現(xiàn)象，成為研究的熱點(diǎn)[24]。本文分別使用主成分分析和偏相關(guān)分析進(jìn)行人工神經(jīng)網(wǎng)絡(luò)的輸入因子篩選，通過比選構(gòu)建NDVI預(yù)測模型。其中采用主成分分析進(jìn)行因子篩選的模型簡稱為ANN-PCA, 采用偏相關(guān)分析進(jìn)行因子篩選的模型簡稱為ANN-PAR。為了比較，加入采用全部預(yù)報(bào)因子作為輸入的人工神經(jīng)網(wǎng)絡(luò)模型，簡稱ANN。表5給出了人工神經(jīng)網(wǎng)絡(luò)模型的率定和驗(yàn)證結(jié)果。

表5 人工神經(jīng)網(wǎng)絡(luò)計(jì)算結(jié)果Tab.5 Artificial neural network results

結(jié)果表明，基于人工神經(jīng)網(wǎng)絡(luò)的NDVI預(yù)測模型能夠取得較高的預(yù)報(bào)精度。其中，ANN-PCA模型在各分區(qū)的率定期的NASH系數(shù)分別為0.68、0.74、0.79、0.75、0.78，檢驗(yàn)期分別為0.67、0.73、0.77、0.75、0.72。模型在各區(qū)域的率定驗(yàn)證結(jié)果均達(dá)到了0.7左右，取得了較高的模擬預(yù)測精度。

其次，從不同神經(jīng)網(wǎng)絡(luò)的模型的結(jié)果可以看出，ANN-PCA模型和ANN-PAR模型均顯著優(yōu)于ANN模型。在率定期，ANN-PCA模型比ANN模型在各分區(qū)的NASH系數(shù)分別高0.17、0.11、0.15、0.11、0.06，在檢驗(yàn)期，分別高0.11、0.10、0.13、0.10、0.12。兩種篩選方法對比，ANN-PCA模型略優(yōu)于ANN-PAR模型，在率定期，ANN-PCA模型比ANN-PAR模型在各分區(qū)的NASH系數(shù)分別高0.03、0.05、0.01、0.04、0.06，在檢驗(yàn)期，分別高0.04、0.02、0.02、0.08、0.05。結(jié)果體現(xiàn)了NDVI預(yù)測模型構(gòu)建時(shí)因子篩選的重要性。流域分區(qū)在率定期和驗(yàn)證期NASH效率系數(shù)如圖5所示。結(jié)果表明，ANN-PCA和ANN-PAR模型由于經(jīng)過因子篩選，避免了冗余數(shù)據(jù)的輸入，有效提高了人工神經(jīng)網(wǎng)絡(luò)的訓(xùn)練效率。兩個(gè)模型的精度顯著高于未經(jīng)過因子篩選的ANN模型。因此，無效或冗余的輸入數(shù)據(jù)是人工神經(jīng)網(wǎng)絡(luò)誤差產(chǎn)生的重要原因，而因子篩選江有助于人工神經(jīng)網(wǎng)絡(luò)模型精度的提高。

5 結(jié) 論

在全球氣候變化的背景下，雅魯藏布江流域的生態(tài)系統(tǒng)也隨之發(fā)生了變化。本研究基于2000-2015年雅魯藏布江流域的MODIS數(shù)據(jù)和氣象站點(diǎn)數(shù)據(jù)，分析雅魯藏布江流域植被的空間分布特征，并通過主成分分析和偏相關(guān)分析方法進(jìn)行影響因子篩選，構(gòu)建了基于人工神經(jīng)網(wǎng)絡(luò)的NDVI預(yù)測模型。主要結(jié)論如下：

(1)NDVI在雅魯藏布江流域整體上植被覆蓋比較穩(wěn)定，呈現(xiàn)不顯著的增長趨勢。Hurst指數(shù)的分析結(jié)果表明，雅魯藏布江干流上游流域(拉孜站以上)呈反持續(xù)性，未來植被覆蓋將可能減少，下游區(qū)域呈持續(xù)性，未來植被覆蓋將持續(xù)增加，且持續(xù)性從上游到下游逐漸增強(qiáng)。

(2)采用主成分分析和偏相關(guān)分析分別對影響雅魯藏布江各子流域NDVI的因子進(jìn)行了分析，結(jié)果表明對NDVI具有顯著影響的降雨和氣溫主要集中在前3個(gè)月，同期氣溫的影響要高于降水的影響。

(3)本文所構(gòu)建的NDVI預(yù)測模型具有較高的模擬與預(yù)報(bào)精度。其中，影響因子的篩選能夠顯著提高人工神經(jīng)網(wǎng)絡(luò)的預(yù)測精度。本研究中，使用主成分分析篩選得到的模擬結(jié)果相比于偏相關(guān)分析得到的模擬果更接近真實(shí)值。

論文所建模型和分析結(jié)論能夠?yàn)槲鞑厣鷳B(tài)環(huán)境治理提供技術(shù)支持。由于人工神經(jīng)網(wǎng)絡(luò)模型的精度依賴于訓(xùn)練數(shù)據(jù)，因此該方法適用于具備較長觀測資料，氣象數(shù)據(jù)較為完整的區(qū)域。同時(shí)缺乏物理機(jī)制的具體分析是人工神經(jīng)網(wǎng)絡(luò)模型的主要缺點(diǎn)，下一步將結(jié)合雅魯藏布江流域的植被類型，將在拉薩河等子流域進(jìn)行典型研究，深入分析其物理機(jī)制，并考慮更多的局地影響因子(如高程等)進(jìn)行綜合分析評(píng)價(jià)。

□