趙健赟，彭軍還，宋芊，張波

(1.青海大學(xué) 地質(zhì)工程系，西寧 810016；2.中國地質(zhì)大學(xué)(北京) 土地科學(xué)技術(shù)學(xué)院，北京 100083；3.北京麥格天寶科技股份有限公司，北京 100043)

0 引言

全球變化及其帶來的各種生態(tài)環(huán)境問題受到各國政府和國際組織越來越多的關(guān)注，已經(jīng)成為地球科學(xué)領(lǐng)域的研究熱點之一[1-5]。植被作為陸地生態(tài)系統(tǒng)的重要組成部分，是連接大氣和土壤的重要紐帶，植被覆蓋是區(qū)域環(huán)境、地理條件和氣候變化等多種因素作用的結(jié)果，植被覆蓋變化是全球變化的重要指示器[6-7]。

歸一化植被指數(shù)(normalized difference vegetation index，NDVI)與植被綠度、物候、光合作用緊密聯(lián)系，是反映全球或區(qū)域植被環(huán)境變化的重要指標(biāo)。目前，MODIS、AVHRR、SPOT、TM/ETM等衛(wèi)星傳感器均提供不同時空分辨率的NDVI時間序列產(chǎn)品，已有很多學(xué)者利用NDVI產(chǎn)品對全球或區(qū)域的植被覆蓋時空變化進行了研究[8-13]，如鄧興耀等[12]利用GIMMS NDVI對中亞干旱區(qū)植被覆蓋的時空變化進行了研究，申麗娜等[10]利用GIMMS NDVI3g對三北防護林的植被覆蓋變化圖譜進行了分析，郭繼凱等[11]利用MODIS NDVI對塔里木河流域植被覆蓋變化及其驅(qū)動因素進行了研究等。

但是，目前大多數(shù)研究利用NDVI時間序列數(shù)據(jù)進行植被覆蓋變化分析時，僅對其進行線性建模，掩蓋了地理現(xiàn)象和過程的真實特征，未考慮NDVI時間序列的非平穩(wěn)性[14-16]。處理非平穩(wěn)序列的傳統(tǒng)方法是先利用對數(shù)、差分等算法將其平穩(wěn)化后，再利用自回歸滑動平均(autoregressive and moving average，ARMA)進行建模處理[17-21]。隨著機器學(xué)習(xí)及人工智能技術(shù)的迅速發(fā)展，人工神經(jīng)網(wǎng)絡(luò)、支持向量機、遺傳算法等為時空非平穩(wěn)序列建模提供了新的方法[22-26]。本文基于人工神經(jīng)網(wǎng)絡(luò)(artificial neural network，ANN)模型，研究顧及NDVI數(shù)據(jù)非平穩(wěn)特性下的時間序列建模方法，利用GIMMS NDVI3g數(shù)據(jù)對青海省1982—2015年的植被覆蓋變化進行建模，并對其未來變化進行預(yù)測。

1 數(shù)據(jù)和研究方法

1.1 數(shù)據(jù)資料

選用美國國家航空航天局提供的1982—2015年15天合成8 km分辨率AVHRR GIMMS NDVI3g數(shù)據(jù)集。該數(shù)據(jù)集融合了NOAA-7/9/11/14/16/17/18/19的AVHRR/2和AVHRR/3兩種不同傳感器數(shù)據(jù)，并且進行了軌道偏移、大氣水汽、輻射、氣溶膠、除云、幾何畸變等校正與處理[27]。

選用中國科學(xué)院計算機網(wǎng)絡(luò)信息中心國際科學(xué)數(shù)據(jù)共享平臺的青海省行政區(qū)劃數(shù)據(jù)，對NDVI3g數(shù)據(jù)集進行裁切獲得青海省NDVI數(shù)據(jù)集，并利用SG濾波算法對數(shù)據(jù)集進行降噪處理[28]。

選用美國航天航空局 SRTM 90 m DEM Version 4共 8景DEM 數(shù)據(jù)，對其進行拼接、裁切等處理后形成青海省DEM 數(shù)據(jù)。

以上數(shù)據(jù)均采用Albers Conical Equal Area投影(南標(biāo)準(zhǔn)緯線：25°S；北標(biāo)準(zhǔn)緯線：47°N；中央經(jīng)線：105°E)。

1.2 時間序列數(shù)據(jù)的非平穩(wěn)性

平穩(wěn)時間序列數(shù)據(jù)一般指序列數(shù)據(jù)中不包含趨勢、周期等變化特征，即期望或均值為常數(shù)，其方差、協(xié)方差等統(tǒng)計量不隨時間和空間的變化而發(fā)生改變[29-34]。設(shè)t時刻i位置的觀測值均值為：

μzi(t)=E[zi(t)]

(1)

若時間序列平穩(wěn)，則在(t+k)時刻(i+h)位置滿足：

(2)

若不滿足式(1)、式(2)，則該數(shù)據(jù)為非平穩(wěn)序列。

1.3 非平穩(wěn)NDVI時間序列模型方法

時間序列過程一般由總體確定性變量和局部隨機變量構(gòu)成，總體變異即為變化趨勢，而局部變異是在時間序列過程中分離出總體趨勢后的隨機變異量，局部變異量一般為平穩(wěn)過程。因此，一個NDVI序列過程可用下式表示[35-36]：

(3)

式中：NDVIi(t)為t時刻i位置的NDVI觀測值；μi(t)為數(shù)學(xué)期望總體趨勢；f(i，t)為非線性函數(shù)；ei(t)為隨機變異量。本文提出基于ANN與ARMA組合的NDVI時間序列模型，流程如圖1所示。

圖1 基于ANN與ARMA的NDVI序列建模流程

考慮到ANN模型具有自組織、自學(xué)習(xí)能力和較強的非線性和容錯性，比較適合于解決復(fù)雜的非線性動態(tài)問題[37]，可以直接建立基于ANN非線性NDVI時間序列模型，其中包含一個隱含層、一個輸出層的ANN函數(shù)模型如下[38-40]：

(4)

非平穩(wěn)NDVI時間序列數(shù)據(jù)建模，涉及到確定ARMA模型和非線性ANN模型參數(shù)的問題，選用最小信息準(zhǔn)則(akaike information criterion，AIC)估計最佳模型參數(shù)：

(5)

2 NDVI時序數(shù)據(jù)建模結(jié)果與分析

2.1 NDVI時間序列數(shù)據(jù)的非平穩(wěn)性檢驗

基于青海省GIMMS NDVI3g數(shù)據(jù)，在像元、區(qū)域級別分別構(gòu)建15天、年尺度NDVI時間序列。為描述構(gòu)建的NDVI序列數(shù)據(jù)特性，繪制對應(yīng)的QQ圖(正態(tài)概率圖)[31-32，36，41]，結(jié)果如圖2所示。利用小波時序分析方法[42-45]，對像元級別的時間序列進行周期特征和趨勢特征進行分解，結(jié)果如圖3所示。

圖2 NDVI序列正態(tài)性檢驗

圖3 NDVI序列季節(jié)和趨勢特征分解

從圖2、圖3可以看出，時間序列數(shù)據(jù)均不符合正態(tài)分布，顯然不滿足式(1)、式(2)的平穩(wěn)條件，數(shù)據(jù)序列具有明顯的趨勢特征，且像元級別的時間序列數(shù)據(jù)具有顯著的周期變化特征。因此，NDVI3g數(shù)據(jù)在像元、區(qū)域級別上均具有明顯的非平穩(wěn)特性。

2.2 NDVI時間序列數(shù)據(jù)的非平穩(wěn)建模

利用ANN與ARMA組合模型、ANN非線性模型對像元年際序列數(shù)據(jù)進行分析，獲得2種模型參數(shù)的AIC指標(biāo)如表1所示。

表1 像元年際尺度的2種模型AIC指標(biāo)

按照最小信息準(zhǔn)則，確定組合模型的參數(shù)為(6，6)，非線性ANN模型的參數(shù)為(4，1，1)。利用上述參數(shù)分別建立ANN-ARMA組合模型和非線性ANN模型，結(jié)果如圖4所示。

圖4 像元尺度的組合模型與ANN模型

由建立的NDVI時間序列模型可以看出，2種模型均能較好地模擬和反映NDVI變化情況。組合模型較為光滑，能夠凸顯植被的總體變化趨勢，而非線性ANN模型與觀測值更為吻合，適宜于描述植被覆蓋的細(xì)節(jié)變化。

2.3 分析與預(yù)測

為定量評價NDVI非平穩(wěn)時間序列模型的精度，在像元、區(qū)域級別年尺度上分別計算ANN-ARMA組合模型和非線性ANN模型的殘差，結(jié)果如表2所示。

表2 像元、區(qū)域級別年尺度的2種模型殘差

從表2可以看出，2種模型的估計值與NDVI原始值之間的差異并不顯著。利用ANN-ARMA組合模型和非線性ANN模型，按照以下方法分別預(yù)測2011—2015年區(qū)域級別上的NDVI值：

(6)

表3 ANN與組合模型的精度檢驗

由表可以看出，ANN-ARMA組合模型的誤差絕對值均小于0.001，比ANN模型誤差小接近一個數(shù)量級，且ANN模型誤差值域分布比組合模型更加離散，離差達到0.006 5，而組合模型離差僅有0.000 1。因此，相對于直接利用ANN建立的非線性模型，ANN-ARMA組合模型具有更好的預(yù)測精度，也能更加準(zhǔn)確地描述NDVI3g數(shù)據(jù)的非平穩(wěn)特征，有利于數(shù)據(jù)的建模與分析。

青海省的氣候暖區(qū)主要分布在海拔2 000 m 以下的河湟谷地地區(qū)，次暖區(qū)主要分布在海拔3 000 m左右的柴達木盆地等地區(qū)，冷區(qū)主要分布在4 000 m以上的祁連山區(qū)和青南高原地區(qū)?？紤]到2 000 m以下河湟谷地面積極小，將青海省劃分為3 000 m以下、3 000～4 000 m和4 000 m以上3個高程帶，利用ANN-ARMA組合模型對青海省2016—2025年各高程分區(qū)的植被覆蓋變化進行預(yù)測，結(jié)果如圖5所示。

從預(yù)測結(jié)果來看，2016—2025年青海省NDVI整體有上升趨勢，3 000 m以下地區(qū)的NDVI值有降低趨勢，3 000 m以上地區(qū)有上升趨勢。而從近年來有關(guān)研究區(qū)植被覆蓋變化的研究成果來分析，2000—2016年青藏高原植被覆蓋改善的地區(qū)要比退化的地區(qū)面積大，其NDVI呈增加趨勢[46]，與預(yù)測結(jié)果吻合，表明ANN-ARMA組合模型的預(yù)測結(jié)果可靠。

圖5 各高程分區(qū)組合建模與預(yù)測

3 結(jié)束語

本文在對青海省GIMMS NDVI3g時間序列數(shù)據(jù)進行非平穩(wěn)分析的基礎(chǔ)上，利用人工神經(jīng)網(wǎng)絡(luò)模型研究了NDVI數(shù)據(jù)非平穩(wěn)時間序列建模方法?；?982—2015年的NDVI3g數(shù)據(jù)建立了青海省植被覆蓋變化ANN非線性模型和ANN-ARMA組合模型，分析了模型精度，并預(yù)測了青海省2016—2025年的植被覆蓋變化規(guī)律，主要結(jié)論如下：

①GIMMS NDVI3g時間序列數(shù)據(jù)具有顯著的非平穩(wěn)特性。

②ANN非線性模型和ANN-ARMA組合模型均能較好地模擬植被變化情況。組合模型較為光滑，能夠凸顯植被的總體變化趨勢，而非線性ANN模型與觀測值更為吻合，適宜于描述植被覆蓋的細(xì)節(jié)變化。

③ANN-ARMA組合模型的精度高于ANN非線性模型精度。

④2016—2025年青海省植被覆蓋整體有上升趨勢，3 000 m以下地區(qū)植被覆蓋有降低趨勢，3 000 m以上地區(qū)有上升趨勢。