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

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

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

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

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

1 材料與方法

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

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

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

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

高光譜觀測(cè)采用的儀器為SVC HR-1024I 全波段地物光譜儀，波長(zhǎng)范圍350 ～2 500 nm，包含 1 024 個(gè)波段，視場(chǎng)有4°標(biāo)準(zhǔn)和25°光纖兩種，本次實(shí)驗(yàn)選用的是4°標(biāo)準(zhǔn)視場(chǎng)角。

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

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

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

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

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

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

1.2.3 銀杏葉片NDVI 值計(jì)算

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

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

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

1.3 研究方法

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

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

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

1.3.2 ARIMA 時(shí)間序列擬合法

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

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

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

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

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

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

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

1.3.3 零均值化處理

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

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

1.3.4 AIC 準(zhǔn)則

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

式中n—樣本數(shù)量，ln—對(duì)數(shù)運(yùn)算符，A—?dú)埐钇椒胶?，p—回歸方程中自變量的個(gè)數(shù)。

2 結(jié)果與分析

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

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

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

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

式（7）中：y—NDVI，x—時(shí)間編號(hào)。

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

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

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

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

2.2 ARIMA 時(shí)間序列擬合結(jié)果

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

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

表4 單樣本Kolmogorov-Smirnov 檢驗(yàn)Tab.4 Single sample Kolmogorov-Smirnov test

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

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

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

表5 ARMA 模型相關(guān)性特征Table 5 ARMA model correlation characteristics

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

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

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

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

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

采用Matlab 軟件對(duì)模型進(jìn)行ARIMA 時(shí)間序列曲線擬合得到曲線，如圖6 所示。

表6 AIC 計(jì)算結(jié)果Table 6 AIC criteria to calculate the results

圖6 銀杏落葉期NDVI 值A(chǔ)RIMA 時(shí)間序列擬合結(jié)果Fig.6 Deciduous period of Ginkgo leaves NDVI values of ARIMA time series fitting results

采用Eviews 軟件對(duì)模型進(jìn)精度評(píng)價(jià)，R2為0.811，擬合曲線與實(shí)際曲線走勢(shì)一致，擬合效果較好。

3 結(jié)論與討論

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

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

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

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

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