基于季節性時間序列模型的西藏那曲牧草生長預測研究

汪書樂 雒偉群 趙益民 嘎桑久美

1.西藏民族大學信息工程學院，咸陽 712082；2.西藏民族大學醫學院，咸陽 712082

草原是西藏生態環境大整體中關鍵的一個部分，草原也是牧民生活的重要依靠，目前歸一化植被指數（NDVI）預測已成為一個可行的辦法來評估草原環境保護的方法?！段鞑刈灾螀^國家生態文明高地建設條例》要求我們保護生態環境，保護好我們的綠水青山。充分了解草場的情況，采用合理的放牧可以幫助我們更好地保護生態環境，NDVI 的預測讓我們對草原的植被情況有著更好地了解，從而來更好的判斷草場的生長情況，來預估草原的保護是否有效。

NDVI 是利用遙感技術對地表覆蓋進行定量評估的指標之一，它是通過計算紅外波段和可見光波段的比值來反映地表植被覆蓋狀況的指標。NDVI 在氣候變化、環境監測、農業生產等領域中具有重要的應用價值。因此，如何精確預測NDVI的變化趨向，對于相關領域的研究和實踐具有重要的意義。植被生長情況是一個關鍵的環境指標，它與氣候變化、生態環境和農業生產等密切相關。對植被生長情況進行準確的預測和監測，對于環境保護、生態建設和農業生產都具有關鍵的意義。歸一化植被指數（NDVI）是衡量植被覆蓋程度的重要指標，它通過反映地表反射率中紅光和近紅外光的比值來反映植被的生長狀況。NDVI數據能夠經過遙感技術獲取，對植被的生長狀態進行監測和預測。

過去的研究中，在NDVI的研究方面，國內外學者采用了多種方法來預測NDVI的變化趨勢。例如人工神經網絡在NDVI 時間序列預測中的應用［1］，利用深度學習模型對歸一化植被指數(NDVI)實現精準預測［2］，通過空間分析、統計分析和機器學習等方法，探討了草地NDVI 與降水量、氣溫和日照時數的相關關系［3］，以青海湖流域2000—2014 年歸一化植被指數作為輸入數據(8 天，500 米)，利用時間序列模型生成預測NDVI。采用不同時間序列模型:趨勢移動平均法、二次指數平滑法、差分自回歸滑動平均模型(ARIMA)預測2015 年NDVI 的空間分布并進行精度比較［4］，這些方法并不能很好地解決季節性變化的問題。因此，本文采用季節性自回歸移動平均模型(SARIMA)來預測西藏那曲市NDVI 的變化趨勢[5-15]。本文的主要貢獻在于采用了一種基于季節性時間序列模型的NDVI預測方法，該方法能夠很好地解決西藏那曲NDVI 季節性變化的問題。本文將該方法應用于西藏那曲市的NDVI 預測中，并通過實驗驗證了該方法的有效性和準確性。

1 相關技術

研究區域：那曲市屬于高海拔草原生態系統，這種類型的生態系統的特點是氣候寒冷干燥，由于海拔高而缺氧。白天和黑夜溫度差距大，經常刮大風。每年的5 月至9 月有一個暖和的時期，被認為是草原牧草的重要生長期。在此期間，降雨量最豐富，高山草甸和高山草甸草原類型茁壯成長。該生態系統每年的綠色植物生長期約為100 天，并且集中在這個溫暖的季節。

遙感數據：MOD13Q1 和MYD13Q1 是兩種種植被指數產品，源自Terra 和Aqua 中分辨率成像光譜輻射計（MODIS）傳感器獲取的數據。它提供了全球尺度植被活動的估計值，空間分辨率為250m。MOD13Q1和MYD13Q1 的NDVI 數據產品每16 天生成一次，兩者的數據產品時間之間還會相差8 天，我們采用的數據是陸地植被產品中的NDVI，它提供了植被冠層綠度的一致時空比較，這是葉面積、葉綠素和冠層結構的復合特性。植被指數產品廣泛用于各種應用，如監測植被動態、干旱評估和作物產量估算。

植被指數NDVI：NDVI 用于估計給定區域的綠色植被情況。它基于健康的植被吸收更多紅光并反射更多近紅外光的原理。NDVI 公式將反射的近紅外光和紅光之間的差值除以兩者的總和：NDVI=(NIR-R)/(NIR+R)。其中NIR是植被反射的近紅外光量，R是植被反射的紅光量。NDVI 公式在分析植被覆蓋方面具有重要意義，因為它可以區分茂密植被區域、稀疏植被區域和裸露地面區域。NDVI 的值范疇為-1～+1，其中接近-1 的值表示裸露的土壤或水，接近0 的值表示植被稀疏或沒有綠色植被的區域，接近+1 的值表示茂密的綠色植被。因此，通過使用NDVI，研究人員可以識別具有不同植被覆蓋水平的區域，并監測植被隨時間的變化，這對于研究生態系統動態、土地利用和氣候變化非常有用。綜上所述，NDVI 能夠評估牧草的生長情況，NDVI 公式是通過遙感影像數據分析植被覆蓋的有力工具。它的意義在于它能夠區分不同層次的植被覆蓋，可用于研究生態系統動態、土地利用和氣候變化。

時間序列模型：時間序列模型是功能強大的工具，通過可識別表現出隨時間變化的數據中的模式和趨勢，以及預測未來值。時間序列模型被寬泛作用于經濟、氣象、環境等領域的數據分析和預測中。在時間序列模型中，ARIMA 模型對于非季節性時間序列數據建模特別有用，而SARIMA 模型能夠更好地處理表現出隨著時間和季節性變化而變化數據。ARIMA模型通過對時間序列數據進行差分來消除非平穩性，進而通過對差分后的序列進行模型選擇、參數估計和模型診斷等步驟來預測未來的數據變化趨勢。ARIMA 模型，需要選擇自回歸項的階數（AR 部分的參數p）、差分的次數（d 參數）和移動平均項的階數（MA 部分的參數q）。這些參數的選擇可以通過觀察自相關性和偏自相關性函數來進行。參數估計能夠運用最大似然估計或最小二乘法來實現。完成參數估計后，還需要進行模型診斷，以確定模型是否滿足假設和是否需要進一步調整。SARIMA 模型的顯示為SARIMA(p,d,q)(P,D,Q)s，其中，p：非季節性自回歸項的階數（即AR 部分的參數）；d：進行非季節性差分的次數；q：非季節性移動平均項的階數（即MA 部分的參數）；P：季節性自回歸項的階數；D：進行季節性差分的次數；Q：季節性移動平均項的階數；s：季節性的周期長度，要將非平穩時間序列（例如表現出趨勢和季節性模式的NDVI）轉換為平穩時間序列，必須使用微分方法。平穩時間序列的特點是具有恒定的平均值且沒有趨勢變化，而非平穩時間序列則表現出向上或向下的趨勢。確保時間序列的穩定性至關重要，應用適當的變換對于實現這一目標至關重要。季節性時間序列通常表現出循環模式，需要通過微分來消除周期性變化。為此，使用了可以處理季節性時間序列數據的SARIMA 模型。在季節性時間序列模型中，時間序列數據可以被看作由一個長期趨勢、一個季節性周期和一個隨機波動組成的線性組合。長期趨勢是時間序列數據隨著時間的變化而整體呈現出的變化趨勢，季節性周期是時間序列數據在特定時間范圍內的重復變化趨勢。

精度檢驗方法：

（1）MAPE:平均絕對百分比誤差（Mean Absolute Percentage Error）

平均絕對百分比誤差（MAPE）是用于評估預測模型的準確性或測量實際值和預測值之間的誤差的指標。它表示為百分比，并計算為實際值和預測值之間的絕對百分比差異的平均值。

（2）RMSE：均方根誤差（Root Mean Square Error）

均方根誤差（RMSE）是用于評估回歸模型或預測模型性能的廣泛使用的指標。它以與原始數據相同的單位測量預測值與實際值的平均偏差。

真實值：y={yi,y2,…yn}，n是模型中的參數數量。

2 時間序列模型的實現

本部分主要介紹實驗設計和結果分析。

季節性時間序列模型的構建主要分為以下幾個步驟：數據收集、數據預處理、模型選擇、模型訓練與評估、結果分析等步驟，如圖1所示。

圖1 時間序列模型構建步驟流程圖Figure 1 Flow chart of time series model construction steps

數據收集：從NASA官網收集2012年到2020年遙感數據產品(如MOD13Q1 和MYD13Q1)中NDVI（歸一化差異植被指數）。數據包含2012 至2020 共9 個年份，每年12×4=48 項數據，數據集由9×48=432 項數據組成。

數據預處理：對收集到的NDVI 時間序列數據進行預處理，包括缺失值處理，異常值，數據平滑等。應用ArcGIS、ENVI 對其中的MOD13Q1 和MYD13Q1 指數產品處理，處理完成后將16天最大值為輸入數據。

缺失值處理：ArcGIS 最大值合成法（MVC）方法用于用相鄰單元的最大值替換柵格數據集中的缺失值。使用ArcGIS中的“像元統計數據”工具來應用此方法。選擇兩張相鄰8 天的MOD13Q1 和MYD13Q1 數據作為輸入柵格數據，配置疊加統計為MAXIMUM，使用最大合成法補全缺失值。得到周期為16 天的NDVI數據集。

異常值處理：利用ENVI 中band math 的工具進行篩選出大于-1 小于1，輸入公式：(b1 lt 0)×0+(b1 ge 0 and b1 le 10000)×b1×0.0001，式中b1 為NDVI 波段，lt意思為小于，ge 意思為大于等于，le 意思為小于等于，整體上式子為小于0 的都乘0，大于等于0、小于等于10000 的都保留本身并乘以0.0001，因此，該公式通過將NDVI 值固定在［-1,1］范圍內，并在有效范圍內重新縮放值，有效地處理異常值。

數據平滑：利用Savitzky-Golay（S-G）濾波算法于NDVI（歸一化植被指數）數據進行平滑處理，簡單介紹S-G 濾波算法的示例表達式：smoothed_nevi=savitzky_golay(b1,window_size,polynomial_order) 其 中smoothed_nevi是平滑后的NDVI數據。b1是原始NDVI數據。window_size 是S-G 濾波的窗口大小，它用于確定在每個時間點應用的數據點數量。polynomial_order是多項式階數，它決定了S-G濾波所使用的多項式的次數。圖2 為展示批量經過MVC 合成為30 天最大值處理，異常值處理后，2012年那曲市各月NDVI數據。

圖2 2012年那曲NDVI 時間序列圖（250m，30天）Figure 2 Naqu NDVI time series chart in 2012(250m,30days)

模型選擇：初步選擇為差分自回歸模型，主要是看加了季節性因素后，兩個模型是否為合適的ARIMA模型和SARIMA 模型。利用赤池信息量準則(AIC)或貝葉斯方法(BIC)對所選的ARIMA 模型和SARIMA 模型進行參數估計，較小的AIC 和BIC 能夠說明模型具有較好的擬合效果，經過比較AIC和BIC值比較之后，我們選擇了SARIMA 模型，因為它在這兩個信息準則下都具有較小的數值。這意味著SARIMA 模型更能有效地描述和預測我們的數據，從而成為最終的模型選擇。AIC的計算公式：AIC=-2×ln(L)+2×k，BIC的計算公式：BIC=-2×ln(L)+k×ln(n)，其中，L是似然函數的最大值。k是模型中的參數數量。n是數據點的數量。

模型訓練與評估：將預先處理后的NDVI 時間序列數據歸類為訓練集和測試集，采用訓練集對選定的模型進行訓練，并應用測試集進行模型的分析。使用評估指標如平均絕對百分比誤差（MAPE），均方根誤差（RMSE）來分析和評估模型的預測性能。預測NDVI的變化趨勢，將模型預測結果與實際觀測值進行比較，制作預測曲線和實際觀測值曲線，查看它們之間的擬合水平圖3、圖4所示。再經過計算RMSE、MAPE來評測所選季節性時間序列模型的預測精度，表1 所示。較低的RMSE和MAPE值說明模型具備較高的預測準確性。

表1 ARIMA模型和SARIMA模型預測精度評價表Table 1 ARIMA model and SARIMA model prediction accuracy evaluation table

圖3 基于ARIMA模型的預測圖Figure 3 Forecast chart based on ARIMA model

圖4 基于SARIMA模型的預測圖Figure 4 Prediction chart based on SARIMA model

結果分析：SARIMA 模型是適合擁有季節性波動的數據來分析其時間序列的建模辦法。首先，確定了數據集是由16天為一個周期的數據組成，一年便是時間周期為s=24，并提取了非季節性成分來調整每個季節周期內部值的變化。建立模型流程如下：第一步，取2012—2020 年NDVI 數據為數據集，對數據集進行差分，平穩時間序列；第二步，利用ADF測試確定序列是否平穩，并通過PACF方法確定了SARIMA模型的p和q值（p=2，q=1）；第三步，利用赤池信息量準則（AIC）或貝葉斯方法(BIC)標準來估算所選差分自回歸平均模型（ARIMA）和季節性自回歸移動平均模型（SARIMA）的參數，并在選擇最佳模型時權衡了模型的復雜度和擬合數據的好壞；第四步，進行殘差測試，利用RMSE、MAPE等指標來評估模型的準確性，ARIMA 模型中MAPE=0.045，SARIMA 模型中MAPE=0.043，ARIMA 模型中RMSE=0.0115，SARIMA 模型中RMSE=0.0108，結果表明SARIMA模型的殘差誤差比ARIMA模型更小，RMSE和MAPE的值較低說明SARIMA模型具有更高的預測準確性。

3 總結與展望

綜上所述，實驗設計和結果分析應包括數據收集、數據預處理、模型選擇、模型訓練與評估、結果分析等步驟，以全面評估基于季節性時間序列模型的NDVI 預測研究的有效性和可靠性。本文探索了兩種常見的時間序列模型來探索NDVI 的預測精度，選擇ARIMA 模型和SARIMA 模型這兩種模型進行理論分析，通過NDVI 的訓練集進行對比實驗，實驗證明SARIMA 模型具備較高的擬合精度，也因為加入了季節性參數，SARIMA 模型在具有季節性數據預測方面，精度更高。這為草原牧草生長和環境監測等方面提供了重要的參考和決策依據。

基于季節性時間序列模型的NDVI預測研究仍有很大的發展空間。一方面，可以進一步優化和改進季節性時間序列模型，如采用多元季節性時間序列模型、基于深度學習的季節性時間序列模型，基于神經網絡相關的季節性時間序列模型研究等以提高不同下預測精度和預測效率。另一方面，可以將季節性時間序列模型與其他遙感數據進行融合，如降雨數據、溫度數據等，用來取得更為精確和充分的結果?？傊?，基于季節性時間序列模型的NDVI 預測研究有著廣闊的應用前景和研究價值，有望將來實現更多的功能和應用。