周柯妙，林 輝，宋仁飛，蔣 仟，陳忠明，杜 凱

（1.中南林業(yè)科技大學 林業(yè)遙感信息工程研究中心，湖南 長沙 410004； 2.林業(yè)遙感大數(shù)據(jù)與生態(tài)安全湖南省重點實驗室，湖南 長沙 410004； 3.南方森林資源經(jīng)營與監(jiān)測國家林業(yè)與草原局重點實驗室，湖南 長沙 410004）

隨著遙感技術的不斷發(fā)展，尤其是近些年來高光譜遙感技術的大力發(fā)展，使快速獲取植物葉片的的精細光譜信息成為了可能。對銀杏Ginkgo bilobaL.葉片的落葉期光譜特征參量進行以宏觀上的時間為自變量的曲線擬合，有助于了解銀杏落葉期葉片光譜特征參量的變化趨勢，從而了解銀杏的生長情況，對其進行決策和控制，并為植被的大尺度遙感動態(tài)監(jiān)測提供了方法參考。

國內(nèi)外學者利用遙感技術對植被的NDVI 進行時間序列監(jiān)測取得了大量的成果。在空間大尺度遙感數(shù)據(jù)方面，夏照華、張克斌等人采用多月多分辨率遙感數(shù)據(jù)，合成NDVI 數(shù)據(jù)，在加上同時期的降水、溫度數(shù)據(jù)，從多角度分析了近些年NDVI 大尺度的時空分布規(guī)律[1]，吳文斌等利用NOAA/AVHRR NDVI 時間序列數(shù)據(jù)，利用兩種比較方法對Savitzky-Golay 濾波法和非對稱性高斯函數(shù)擬合法進行了比較分析[2]；蔣雪冰等采用MODIS NDVI 影像和野外實測的反射率數(shù)據(jù)，利用線性內(nèi)插和擴展卡爾曼濾波組合重構方法進行重構處理[3]；宋盼盼等基于最小二乘原理，采用多項式擬合方法對中稻、晚稻NIR、R、NDVI 3 種植被指數(shù)進行時間序列特征曲線擬合[4]； 在地物高光譜遙感方面，盧霞等利用ASD 光譜儀，對滸苔水體的反射率數(shù)據(jù)進行連續(xù)測定，再提取計算滸苔水體的光譜特征參量，進而對其進行時間序列擬合研究[5]；同時，盧霞等還對ASD 光譜儀測定獲取滸苔水體的反射率數(shù)據(jù)，對滸苔藻體時間序列高光譜特征參量進行曲線擬合[6]。但在關于地表植被的NDVI 時間序列模擬的研究中，大部分都是使用的大尺度多光譜遙感數(shù)據(jù)，少有用地面高光譜數(shù)據(jù)作分析；而在用地面高光譜數(shù)據(jù)來分析特征參量的數(shù)據(jù)模擬研究中，又鮮有對喬木樹種的研究。

以銀杏為研究對象[7]，利用SVC HR-1024I 全波段地物光譜儀，對銀杏落葉期葉片數(shù)據(jù)進行周期為15 d 的高光譜觀測，對獲取的高光譜數(shù)據(jù)進行去重處理，提取特征參數(shù)，并分別進行二次函數(shù)擬合及ARIMA 時間序列擬合，得到最優(yōu)擬合方法。

1 材料與方法

1.1 研究區(qū)概況及位置圖

研究區(qū)位于湖南省長沙市天心區(qū)中南林業(yè)科技大學植物園內(nèi)，112°59′31″E，28°8′11″N，為亞熱帶季風氣候，四季分明，春溫多變，夏秋多晴，嚴冬期短，暑熱期長的特點，適合植物生長。植物園內(nèi)植物種類繁多，銀杏作為主要行道樹之一，廣泛分布在植物園內(nèi)。

1.2 實驗數(shù)據(jù)獲取及數(shù)據(jù)預處理

1.2.1 高光譜數(shù)據(jù)獲取

高光譜觀測采用的儀器為SVC HR-1024I 全波段地物光譜儀，波長范圍350 ～2 500 nm，包含 1 024 個波段，視場有4°標準和25°光纖兩種，本次實驗選用的是4°標準視場角。

銀杏為銀杏科、銀杏屬落葉喬木，4月初發(fā)芽，每年9月底至12月為落葉期，葉片逐漸由綠變黃并掉落，至1月初完全落葉。實驗選在銀杏葉片落葉期（9月中旬至12月中旬）進行，選取3 株健康、生長環(huán)境相同的中齡銀杏為葉片采集對象，以15 d為觀測周期，利用高枝剪對其冠層向陽處葉片進行采集，將采集到的葉片置于保鮮袋，迅速將葉片帶到室驗室內(nèi)，進行高光譜測量，每株樹觀測5 次。將采集到的進行整理與篩選，剔除異常光譜曲線后，選取3 株樹的光譜反射率均值為實驗數(shù)據(jù)，得到銀杏落葉期各期光譜值（圖1）。

圖1 銀杏落葉期各時期光譜曲線圖Fig.1 Spectral curves of Ginkgo leaves in various periods

1.2.2 高光譜數(shù)據(jù)重疊校正

由于儀器在測量時有的波段進行了重復測量，光譜曲線在波長967 ～1 013 nm 部分有一段重疊曲線，因此需要對原始數(shù)據(jù)進行重疊校正，使光譜曲線變成一條平滑且無重復點的曲線，校正后結果如圖2 所示。

圖2 銀杏落葉期光譜曲線重疊校正結果Fig.2 Overlapping correction results of spectral curve of Ginkgo leaves in various periods

由圖2 可以看出光譜曲線重疊部分已經(jīng)被去除，光譜波段也由原始的1 024 個減到了960 個。

1.2.3 銀杏葉片NDVI 值計算

歸 一 化 植 被 指 數(shù)（Normalized difference vegetation index，NDVI）是反映作物長勢和營養(yǎng)信息的重要參數(shù)之一[8]，植物生長狀態(tài)健康、長勢旺盛的時候，其NDVI 值比較高；但當植物處于不健康狀態(tài)，或者趨于衰老的時候，其NDVI值會隨之降低，NDVI 計算公式為：

其中，NIR—近紅外波段的反射率，R—紅光波段的反射率。

不同樹種由于反射能力不同，而NIR 與R 所對應的反射率也不一樣。對于SVC 全地物光譜儀獲取的光譜值數(shù)據(jù)，目前普遍采用波長719 nm 處的反射率為NIR 處反射率，波長680 nm 的反射率為R 的反射率[5]，據(jù)公式（1）計算銀杏落葉期各期NDVI 值（表1）。

1.3 研究方法

1.3.1 二次函數(shù)擬合法

最小二乘法中的二次函數(shù)擬合法，原理是所選擇的回歸模型應該使所有觀察值的殘差平方和達到最小[9]，其模型為：對給定的數(shù)據(jù)列(xi,yi)，i=0，1…，m，滿 足y=a0x2+a1x+a2，a0，a1，a2的計算如公式（2）（3）（4）。

表1 銀杏落葉期各時期NDVI 值Table 1 NDVI in the deciduous period of Ginkgo leaves

1.3.2 ARIMA 時間序列擬合法

差分自回歸移動平均模型（Autoregressive integrated moving average model, ARIMA）擬合是根據(jù)實測的數(shù)據(jù)[10]，通過對原始不平穩(wěn)序列進行差分使序列平穩(wěn)，再通過曲線擬合對系統(tǒng)進行客觀的描述，可對時間序列進行分析、預測、決策和控制等[12]。一個ARIMA 過程包括自回歸移動平均過程（ARMA）以及ARIMA 過程，其中ARMA 模型有以下三種形式[11]：

一是自回歸模型（AR：Auto-regressive），其基本形式為：

式中εt為白噪聲序列，此時稱時間序xt為服從p階的AR 模型；

二是移動平均模型（MA：Moving-average）， 其基本形式為：

式中εt為白噪聲序列，此時稱時間序列xt為服從q 階的MA 模型；

三是自回歸滑動平均模型（ARIMA），其初始序列一般不具有平穩(wěn)性，對序列進行d 差分后使序列平穩(wěn)，其基本形式為：

式中εt為白噪聲序列，此時稱xt為服從(p,d,q)階的ARIMA 模型。

1.3.3 零均值化處理

零均值化處理能去除直流分量的影響，避免在后面的處理中引入不必要的分量從而對結果造成影響。零均值化處理計算公式如公式（5）。

—零均值化處理后數(shù)據(jù)，xi—原始數(shù)據(jù)，—原始數(shù)據(jù)平均值。

1.3.4 AIC 準則

最小信息量準則（An information criterion，AIC 準則）是一種考評綜合最優(yōu)配置的指標，也是擬合精度和參數(shù)位置個數(shù)的加權函數(shù)，在模型定階過程中，使AIC 值達到最小值的值被認為是最優(yōu)模型階數(shù)。在現(xiàn)實生活中，要求模型越簡單越好，通常階數(shù)≤2。AIC 計算公式如公式（6）。

式中n—樣本數(shù)量，ln—對數(shù)運算符，A—殘差平方和，p—回歸方程中自變量的個數(shù)。

2 結果與分析

2.1 二次函數(shù)擬合結果

銀杏落葉期葉片NDVI 值擬合過程采用Matlab 與Eviews。在Eviews 中 對 銀 杏 落 葉 期NDVI 值（N）與其所對應的時間編號（T）進行方程二次函數(shù)擬合，根據(jù)公式（2）（3）（4）得到的參數(shù)結果如表2。

表2 銀杏落葉期葉片NDVI 二次函數(shù)擬合分析結果表Table 2 Analysis of NDVI quadratic function fitting of leaves of Ginkgo leaves

由表2可知模型二次項系數(shù)為-0.022 1，一次項 系數(shù)為0.054 7，常數(shù)項為0.710 6，即擬合方程為：

式（7）中：y—NDVI，x—時間編號。

在Matlab 中對運用最小二乘法對銀杏落葉期葉片的NDVI 值進行曲線估計得到曲線，如圖3所示：

圖3 銀杏落葉期NDVI 值二次函數(shù)擬合結果Fig.3 Fitting results of quadratic function to NDVI values in deciduous period of Ginkgo leaves

如圖4 中擬合曲線是一條平滑的曲線，在時間編號范圍1 ～3 之間時，預測值與實際值大體相符；在時間編號范圍3 ～4.5 之間時，預測值小于實際值，且實際值不是平緩下降，而是在編號4之后才開始有明顯的下降趨勢；在編號范圍4.5 ～7之間，下降的速度慢慢變緩，而在編號4.5 ～6.5之間，預測值略高于實際值。

從擬合的結果看，R2為0.926，說明擬合效果較好；根據(jù)兩個系數(shù)項和常數(shù)項的t 統(tǒng)計量取值情況，研究中所用到的NDVI 值個數(shù)為7 個，所以這些t 統(tǒng)計量要服從自由度為5 的t 分布，根據(jù)查得的t 分布表[11]，當t=5 時，檢驗水平為0.1、0.05和0.01 時所對應的臨界值分別為2.015、2.571 和4.032，本次擬合中3 個t 統(tǒng)計量只有常數(shù)項的t 統(tǒng)計量顯著，2 個系數(shù)項的t 統(tǒng)計量均遠小于檢驗水平的臨界值，且常數(shù)項的P值小于0.05，而兩個系數(shù)P值遠大于0.05，這說明所構建模型擬合結果具有隨機性[13]。

2.2 ARIMA 時間序列擬合結果

在Matlab 中利用公式（5）對原始NDVI 數(shù)據(jù)進行零均值化處理，得到一組新的數(shù)據(jù)NDVI。在SPSS 軟件中對得到的新的數(shù)據(jù)組進行單樣本K-S檢驗，看數(shù)據(jù)組是否服從正態(tài)分布，檢驗結果如表3 和表4 所示。

表3 K-S 檢驗描述性統(tǒng)計量表Table 3 K-S test descriptive statistical scale

表4 單樣本Kolmogorov-Smirnov 檢驗Tab.4 Single sample Kolmogorov-Smirnov test

由表3 表4 可知，數(shù)組雙側漸進顯著性均大于0.1，數(shù)組服從正態(tài)分布。

對數(shù)列進行模型識別，即計算時間序列其自相關函數(shù)和偏自相關函數(shù)的拖尾及截尾情況[14]。自相關系數(shù)（ACF）指的是同一事件在兩個不同時期之間的相關程度，目的是度量過去的事件對現(xiàn)在的時間所造成的的影響；偏自相關系（PACF）指的是單純考慮單獨事件某一時期對另一時期的影響或相關程度，只考慮事件兩個時期之間的相關性程度。ARMA 模型三種形式的相關性特征及模型選擇，如表5 所示。

由于序列本身具有不平穩(wěn)性，先將數(shù)據(jù)進行差分使序列平穩(wěn)，在Matlab 軟件中差分后的樣本序列的自相關系數(shù)和偏自相關系數(shù)，如圖4 和圖5 所示：

表5 ARMA 模型相關性特征Table 5 ARMA model correlation characteristics

圖4 樣本數(shù)據(jù)自相關系數(shù)Fig.4 Sample data autocorrelaion coefficient

圖5 樣本數(shù)據(jù)偏自相關系數(shù)Fig.5 Sample data partial self-correlation coefficient

根據(jù)NDVI 值零均值化后的樣本序列所計算得到的自相關系數(shù)以及偏自相關系數(shù)結果來看，自相關函數(shù)自延遲1 階后，在2 倍標準差的范圍內(nèi)波動，表明數(shù)列短期內(nèi)相關，因此自相關系數(shù)可視為拖尾；偏自相關函數(shù)延遲1 階的偏自相關函數(shù)顯著大于標準差，其他的自相關系數(shù)在2 倍標準差范圍內(nèi)做小值的波動，非零相關系數(shù)衰減為小值的過程十分突然，自延遲3 階后相關系數(shù)衰減為0，即偏自相關系數(shù)視為截尾；根據(jù)ARMA 模型的相關特征，初步將模型視為AR（3）模型，由于序列組數(shù)據(jù)比較少，為了提高擬合的準確性，具體模型的識別及定階通過AIC 準則進行計算確定[15]。

AIC 準則計算根據(jù)公式（6）在Matlab 軟件中進行，計算結果如表6：

由表6 可知，當p=2，q=2 時，AIC 值最小，此時模型為最優(yōu)，即ARIMA 過程中的ARMA 過程為ARMA（2，2）模型，加上一次差分過程，綜合模型為ARIMA（2,1,2）。

采用Matlab 軟件對模型進行ARIMA 時間序列曲線擬合得到曲線，如圖6 所示。

表6 AIC 計算結果Table 6 AIC criteria to calculate the results

圖6 銀杏落葉期NDVI 值ARIMA 時間序列擬合結果Fig.6 Deciduous period of Ginkgo leaves NDVI values of ARIMA time series fitting results

采用Eviews 軟件對模型進精度評價，R2為0.811，擬合曲線與實際曲線走勢一致，擬合效果較好。

3 結論與討論

以中南林業(yè)科技大學校植物園為研究區(qū)，在銀杏的落葉期內(nèi)，以健康生長的中齡銀杏葉片為研究對象，以半個月為周期，利用SVC HR-1024I全波段地物光譜儀對其進行定點的高光譜數(shù)據(jù)采集，對采集的數(shù)據(jù)進行數(shù)據(jù)預處理后，根據(jù)計算觀察到的NDVI 趨勢，選用二次函數(shù)曲線擬合和ARMA 時間序列擬合法對各期的NDVI 值進行擬合，得到銀杏落葉期間NDVI 值的動態(tài)變化規(guī)律，得出的主要結論如下：

1）二次函數(shù)擬合結果為NDVI=-0.022 1T2+ 0.054 7T+0.711，決定系數(shù)R2為0.926，t 值不顯著，樣本結果隨機性大，不具廣泛性；

2）ARIMA 時間序列擬合結果為ARIMA（2,1,2）模型，R2為0.811，樣本結果與實際趨勢大致一致，擬合效果較好；

3）二次函數(shù)擬合和ARIMA 時間序列擬合都能較好地對本次實驗數(shù)據(jù)進行擬合，但由于二次函數(shù)擬合結果具有很大的隨機性雖然二次函數(shù)擬合結果的決定系數(shù)大于ARIMA 時間序列擬合結果，但二次函數(shù)擬合結果具有很大的隨機性，總體上講ARIMA 時間序列擬合的方法要優(yōu)于二次函數(shù)擬合方法。

基于高光譜遙感數(shù)據(jù)的喬木樹種動態(tài)一直是研究的一個難點和重點。本研究采用地面實測高光譜數(shù)據(jù)，通過兩種方法對銀杏落葉期葉片NDVI值進行曲線擬合，為植被大尺度的高光譜空間監(jiān)測提供方法參考，在下一步的研究中，將嘗試用其他高光譜特征參量對銀杏及其他樹種的光譜特征參量進行年度上的微觀變化分析。