非平穩NDVI時序數據人工神經網絡建模方法

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

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

0 引言

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

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

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

1 數據和研究方法

1.1 數據資料

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

選用中國科學院計算機網絡信息中心國際科學數據共享平臺的青海省行政區劃數據，對NDVI3g數據集進行裁切獲得青海省NDVI數據集，并利用SG濾波算法對數據集進行降噪處理[28]。

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

以上數據均采用Albers Conical Equal Area投影(南標準緯線：25°S；北標準緯線：47°N；中央經線：105°E)。

1.2 時間序列數據的非平穩性

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

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

(1)

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

(2)

若不滿足式(1)、式(2)，則該數據為非平穩序列。

1.3 非平穩NDVI時間序列模型方法

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

(3)

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

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

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

(4)

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

(5)

2 NDVI時序數據建模結果與分析

2.1 NDVI時間序列數據的非平穩性檢驗

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

圖2 NDVI序列正態性檢驗

圖3 NDVI序列季節和趨勢特征分解

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

2.2 NDVI時間序列數據的非平穩建模

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

表1 像元年際尺度的2種模型AIC指標

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

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

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

2.3 分析與預測

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

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

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

(6)

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

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

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

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

圖5 各高程分區組合建模與預測

3 結束語

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

①GIMMS NDVI3g時間序列數據具有顯著的非平穩特性。

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

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

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