基于數字影像的北京松山森林物候模擬及其與氣象因子的關系*

李潤東 田文東 于海群 李鑫豪 靳 川 劉 鵬 查天山 田 赟

(1.北京林業大學水土保持學院 北京 100083；2.上海勘測設計研究院有限公司 上海 200335；3.北京林業大學水土保持國家林業和草原局重點實驗室 北京 100083；4.北京市房山區林果科技服務中心 北京 102400；5.北京市園林綠化規劃和資源監測中心(北京市林業碳匯與國際合作事務中心) 北京 100013)

植被物候一般指植被生長發育過程中，如葉片展開、開花和葉片衰老等周期性事件的發生時間和持續時間(Zhouetal.，2013)，不同空間尺度上的植被物候可觀測提供有關生物圈與大氣之間能量交換、碳水通量動態變化的重要信息(張學霞等，2003)。植被物候變化是通過光、溫度、土壤濕度和養分有效性等的微妙變化來協調的，對生態系統生產力具有控制作用，高精度模擬氣候變化條件下的植被物候變化至關重要(代武君等，2020；王連喜等，2010)。目前，在不同時空尺度上，植被物候與氣候變化之間的關系仍然存在不確定性，傳統植被物候監測依賴于人類對離散物候事件或物候期的直接觀察，雖然比較客觀準確，但費時費力，成本較高(高琪等，2019)；遙感觀測具有連續觀測、覆蓋度高等優點，但獲得的物候指標值與地面實際觀測值存在偏差，同時空間分辨率也較低。隨著圖像處理技術逐步成熟，數字相機影像作為物候監測的新方法，彌補了人工觀測和遙感觀測的不足，通過連續拍攝生態系統高時空分辨率影像，已經應用于全球多個生態系統中(Campilloetal.，2008；Browningetal.，2017)。

與此同時，渦度相關技術可以直接測定生態系統-大氣間的碳交換，具有準確、連續、非破壞性等優點，監測數據可用于檢驗各種基于碳過程的模型，大力推動了陸地生態系統碳循環研究的發展(Hinko-Najeraetal.，2016；紀小芳等，2019；徐麗君等，2011)。近年來，植被物候與生態系統碳循環之間的關系成為生態學家們關注的熱點，物候與碳循環的季節變化密切相關，被認為是氣候變化的一項重要生物學指標，人們逐漸利用數字相機技術進行物候監測和生態系統碳循環研究(Keenanetal.，2014；Richardsonetal.，2013；周磊等，2012a)。以往研究表明，數字重復攝影能夠量化落葉闊葉林、草地、農作物的光合作用期持續時間和生態系統總初級生產力(gross primary productivity，GPP)(Zhangetal.，2019)，數字相機衍生的相對綠度指數(relative green chromatic coordinate，Gcc)數據與碳通量數據高度相關(Kangetal.，2016)。對于落葉林，Gcc對生長季開始(start of growthing season，SOS)時間的估計值與野外人工觀測和衛星衍生SOS的估計值非常匹配，與森林生態系統的衛星衍生數據相比，數字相機綠度指數衍生的物候數據與野外人工觀測更相關(Kurcetal.，2010)。通過數字相機獲得的綠度指數，能夠準確描述森林物候年際變化，并與GPP有較強相關性，這說明數字相機不僅可提供在生態系統尺度上較為準確的植被物候信息，而且還可用來測算碳通量精度，同時也為研究植被物候的環境控制因子提供了新方法(Ahrendsetal.，2009；Richardsonetal.，2009)。有學者將數字相機數據與氣象數據相結合研究植被物候與氣象因子之間的關系，結果發現植被物候與環境因子關系密切，溫度、輻射等均對植被物候具有重要影響，降雨對植被物候具有“觸發”作用(周惠慧等，2016；周磊等，2012b)。本研究以北京松山落葉闊葉林生態系統為例，基于2019年數字相機時間序列數據，提取森林植被生長期中的關鍵指標，并將其與通量觀測數據擬合估算出的物候指標進行對比分析，以檢驗物候模擬精度；同時結合觀測站點環境因子數據，分析林地物候與環境因子間的關系，比較不同環境因子對植被物候的影響程度，以期提高植被固碳模型和區域碳固定模擬的準確性。

1 研究區概況與研究方法

1.1 研究區概況

研究地點位于北京松山國家級自然保護區內(115°47′11″E，40°30′48″N)，占地面積4 671 hm2，海拔1 165 m，距北京市區約104 km。屬大陸性季風氣候，是暖溫帶與中溫帶、半干旱與半濕潤的過渡地帶，年均氣溫8 ℃，年降水量424.6 mm，主要集中在7、8月，占年降水量的62%。年日照時數2 726 h，日照率62%，植物生長期160天左右，生長季葉面積指數3.65。土壤類型為棕壤，pH 6.34，土壤有機質含量149.75 g·kg-1。林分平均冠層高度4 m，結構為復層結構，喬木層包括核桃楸(Juglansmandshurica)、大果榆(Ulmusmacrocarpa)、華北五角楓(Acertruncatum)、大葉白蠟(Fraxinusrhynchophylla)、蒙古櫟(Quercusmongolica)等；灌木層包括繡線菊(Spiraeasalicifolia)、鐵線蓮(Clematisflorida)等；草本層包括異穗薹草(Carexheterostachya)、等齒委陵菜(Potentillasimulatrix)、牛扁(Aconitumbarbatumvar.puberulum)等(李潤東等，2020)。

1.2 數據獲取

本研究主要數據來源于近地表相機記錄的Gcc、渦度相關技術觀測的通量數據、歸一化植被指數(normalized difference vegetation index,NDVI)數據以及微氣象觀測系統觀測的空氣溫度(air temperature，Ta)、土壤溫度(soil temperature，Ts)、光合有效輻射(photosynthetically active radiation，PAR)、飽和水氣壓差(vapor pressure deficit，VPD)和降水量(precipitation，P)。

數字相機安裝在觀測塔上高度20 m處，圖像采集于2019年，從當地時間早上9點到下午4點每30 min采集1次圖像，并自動存儲在數據采集器(CR-1000X)中，格式為JPG。松山觀測站拍攝圖像樣例如圖1所示。Gcc數據依據所采集圖像提取，計算過程詳見數據處理部分，最終獲取Gcc數據樣本365個。NDVI數據由SRS光譜反射傳感器(安裝在觀測塔上高度6 m處)測得的入射光和反射光計算得到。渦度相關通量數據采用2019年GPP數據，由SmartFlux在線計算模塊計算每30 min數據的平均值。微氣象觀測數據包括空氣濕度和溫度、降雨、土壤含水量(soil volume water content，SWC)、土壤濕度和溫度、光合有效輻射、風向、風速等。各種氣象因子傳感器每分鐘采集1次數據，由采集器(CR-1000X)記錄，各種數據每周下載。最終計算空氣溫度、土壤含水量、光合有效輻射等日均值以及降雨量日累計值，得到各環境因子數據樣本365個。

圖1 北京松山落葉闊葉林樣品攝影(a)及森林生長動態(b)Fig.1 Sample photography(a)and broad-leaved forest growing dynamic(b)in Songshan,Beijing紅色多邊形表示感興趣區域，用于計算相對綠度指數。Red polygon represents the ROI (region of interesting)used to calculate relative green chromatic coordinate (Gcc).

1.3 數據處理

1.3.1 相對綠度指數(Gcc)計算 首先篩選采集的圖像，剔除因降雨等原因造成的質量較差的圖像，采用相鄰日期圖像進行插補。構建數字圖像時間序列后，選取一張質量較高的圖像作為模板，在模板圖像上定義合適的感興趣區域(region of interesting，ROI)。利用ROI對所有圖像進行掩蔽，計算每個像素的RGB顏色通道信息(Zhou，2019)。本研究計算白天9:00 AM—4:00 PM每張圖像內ROI的Gcc，公式如下：

(1)

式中：Gdn、Rdn和Bdn分別為ROI中像素點R、G、B光譜波段的平均值。

為了降低高頻Gcc數據噪聲并檢測季節周期，采用3天移動窗口法對Gcc時間序列進行平滑處理，稱為第90百分位法(Richardsonetal.，2018)。對每30 min的Gcc結果取平均值，得到每日的值。

1.3.2 時間序列數據處理 日時間序列環境因子數據包括Ta、Ts、SWC和PAR等，均采用9天窗寬的移動平均法進行過濾(Zhou，2019)。為了確定Gcc與環境因子間的時間序列關系，Gcc的每日值也采用相同移動平均法進行平滑處理。通過繪制生長曲線和檢驗Pearson相關系數分析Gcc與環境因子間的相關關系。

1.3.3 時間序列擬合與物候指標計算 基于平滑插值的時間序列數據，采用Klosterman方法對NDVI時間序列數據、Gcc數據以及GPP數據進行平滑和插值，以擬合森林生長季節軌跡(Klostermanetal.，2014)。地面通量實測數據定義當生長初期第1次出現連續3天總GPP小于夏季最大碳固定的5%時，為生長季的開始日期；而當生長末期第1次出現連續3天GPP小于夏季最大碳固定的5%時，為生長季的結束日期(Zhaetal.，2009；Xieetal.，2016)。從GPP中推導出的樣本生長曲線和提取的主要物候指標如圖2。

圖2 從GPP中推導出的樣本生長曲線和提取的主要物候指標Fig.2 A sample growth curve and key phenological indicators derived from GPPSOS：生長季開始Start of growing season；POP：峰值位置Position of peak value (maximum)；MSP：春季平均值Mean spring value；MGS：生長季平均值Mean growing season value；RSP：春季生長速率Rate of spring green up；RAU：秋季衰老速率Rate of autumn senescence；EOS：生長季結束End of growing season；MAU：秋季平均值Mean autumn value；LOS：生長季長度Length of growing season ；PEAK：峰值Peak value (maximum).下同The same below.

1.3.4 數字相機數據與地面通量數據和NDVI數據的差異 本研究認為，實測GPP時間序列數據提取的物候指標為真實物候指標，對Gcc、NDVI和GPP時間序列數據采用Klosterman方法擬合(周玉科，2018)，并計算各自生長季物候指標，進而計算所提取物候指標的相對差異，用于評價基于數字相機和SRS-NDVI數據的物候指標精度，公式如下：

(2)

式中：P_VIi表示Gcc與NDVI時間序列數據提取的第i個物候指標的值；P_GPPi表示GPP時間序列數據提取的第i個物候指標的值。

1.3.5 統計分析 植被指數能夠反映植被生長狀況，而植被生長發育受多種環境因子影響，這些環境因子對植被生長的影響程度需要準確估計。采用Pearson相關分析方法判定Gcc與各環境因子間的關系，看變量間有無一定程度的相關性。將松山觀測站環境因子數據利用主成分分析法(principal component analysis，PCA)進行分析。主成分分析是一種多元統計技術，其采用正交變換法將一組相關變量轉換成一組正交的、不相關的軸，具有降低數據空間維數的作用，是識別高維變量線性組合最常用的技術(Yaoetal.，2012)。為了定量分析Gcc與各環境因子間的關系，采用全子集回歸分析法篩選Gcc與環境因子間的變量方程，利用貝葉斯信息準則(Bayesian information criterion)和調整決定系數(adjust coefficient of determination)對全子集篩選效果進行評價，最終建立Gcc與環境因子間的回歸方程(譚丞軒等，2020)。

2 結果與分析

2.1 相對綠度指數與環境因子的關系

Gcc的變化趨勢與Ta、Ts的變化趨勢相似，Pearson相關系數分別為0.88和0.86(P<0.01)(圖3a、b，表1)。Gcc與SWC相關性較低(圖3c，表1)(R=0.29,P<0.01)，可能是降雨突發性影響土壤水分變化所致。對比降雨量(圖3f)與土壤水分發現，土壤水分對降雨量有正響應關系，如7—8月時SWC隨降雨增加而增加。PAR的季節變化也表現出與Gcc的一致性(圖3d，表1)，Pearson相關系數R=0.45(P<0.01)。VPD與Gcc的季節變化趨勢一致(圖3e，表1)(R=0.65，P<0.01)。松山森林生態系統的Gcc與降雨雖然相關性較低(圖3f，表1)(R=0.25，P<0.01)，但是降雨對森林變綠期具有“觸發”作用，即降雨能夠使土壤中微生物呼吸加強，促使土壤有機物分解，釋放出植被生長發育所需養分，進而促進植被生長。Gcc隨降雨變化呈現波動狀態，這表明基于數字相機的Gcc具有可實時捕捉到森林生長動態變化的能力。

表1 相對綠度指數與環境因子的相關性①Tab.1 Correlation relationships between Gcc and environmental factors

圖3 相對綠度指數與環境因子的關系Fig.3 Relationships between Gcc and environmental factors

2.2 相對綠度指數與環境因子的回歸分析

從上述得知，Ta、Ts、PAR、P、SWC等均對植被生長發育存在不同程度影響，但任何一個環境因子都不能單獨對植被生長起決定性作用，各環境因子信息相互疊加，存在相關關系，作為一個整體影響植被生長發育。利用主成分分析法(PCA)分析松山觀測站環境因子之間的相關性(各變量數據樣本365個)，根據特征值大于1的原則選取2個主成分，2個主成分對總變異的解釋貢獻率分別為57.03%和19.81%，累計貢獻率為76.83%。各環境因子在2個主成分中的荷載如表2所示，荷載越高表示貢獻率越高。由表2可知，在第1主成分中各指標載荷均為正值，說明環境因子與第1主成分具有一定相關性。土壤體積含水量在第2主成分中的載荷最大且為正值，說明第2主成分主要代表土壤體積含水量。

表2 環境因子在2個主成分中的荷載Tab.2 Environmental factors’load in the two principle components

采用全子集回歸分析法，選擇Ta、Ts、PAR、P、SWC為自變量，Gcc為因變量進行分析(各變量數據樣本365個)。篩選不同數目自變量的隨機組合方式，根據貝葉斯信息準則和調整決定系數篩選出最優組合，結果如表3所示，當組合變量包括Ta、Ts、VPD、PAR時調整決定系數最大，此時方程為最優模型，最終得到Gcc與環境因子間的回歸方程如下：

表3 全子集篩選最佳組合統計結果Tab.3 Best combination result after total subset selection

Gcc=0.34+(1.68Ta-0.39Ts+

0.015PAR-9.89VPD)×0.001。

(3)

2.3 Gcc、GPP和NDVI數據提取的物候指標對比

物候學是全球生態學和陸地生態系統碳循環研究的新線索，在環境監測和管理中具有重要應用價值(Berraetal.，2021)。本研究比較基于數字相機時間序列數據提取的Gcc數據、地面實測的通量數據和NDVI數據的關鍵物候事件時間序列，關鍵物候指標不僅包括生長季開始和結束時間，還包括春季生長速率和秋季衰老速率。采用Klosterman方法對Gcc時間序列數據、NDVI時間序列數據和GPP數據進行擬合并計算物候指標，結果如圖4所示，可以看出，各時間序列數據提取的物候指標差異明顯。通過物候模型擬合提取物候指標的相對差異對比結果見表4。Gcc與GPP提取物候指標對比，除生長季結束時間和秋季衰老速率的相對差異小于0.1外，其余物候指標的相對差異均大于0.1，說明基于Gcc提取的物候指標與GPP差異較大，其中相對差異最大的是春季生長速率，為12.5。這是因為春季植被生長初期光合作用已經進行，但葉子較小，導致數字相機拍攝圖像綠色波段值較低。隨著植被快速生長發育，綠度值快速增長，基于Gcc提取的生長季開始時間較GPP晚，春季生長速率較大，但是秋季物候指標提取結果十分相近，說明數字相機對植被葉子變色過程捕捉敏感，Gcc提取的生長末期物候指標精度大于生長初期。NDVI與GPP提取物候指標對比，生長季開始時間相差6天，相對差異較小，說明SRS-NDVI 測量儀測定的NDVI數據能夠敏感捕捉到植被生長狀態變化。生長季結束時間與GPP的相對差異為0.021，說明提取的生長末期物候指標精度大于生長季初期，與Gcc提取物候指標結果一致。

表4 Gcc、NDVI和GPP提取物候指標的相對差異Tab.4 Relative difference between key phenological indicators derived from Gcc,NDVI and GPP

圖4 利用Gcc(a)、GPP(b)和NDVI(c)反演森林生長曲線和主要物候指標Fig.4 Gcc (a),GPP (b)and NDVI (c)were used to invert forest growth curve and key phenological indicators

3 討論

3.1 相對綠度指數與環境因子的關系

本研究中，Gcc與溫度、光合有效輻射等環境因子間存在較強相關性，一些研究結果也表明，植被生長對溫度、光合有效輻射等環境因子具有正響應關系(劉鑫等，2019；Seyednasrollahetal.，2020)，降雨對Gcc具有“觸發”作用(周磊等，2012b；Berraetal.，2019)。在逐步回歸分析中，溫度是Gcc的主要解釋變量，這可能是因為溫度升高會使生長季開始時間提前(Huangetal.，2018)。SWC與植被生長間存在復雜的相互作用，隨降雨發生出現明顯波動，說明降雨可以有效補充土壤水分。土壤溫度對提高植被生產力具有積極影響(徐滿厚等，2013)，與本研究逐步回歸結果一致，說明土壤溫度在控制植被綠度方面起著關鍵作用。飽和水氣壓差可表征空氣干燥程度，VPD與Gcc的Pearson相關系數為0.65，說明VPD對Gcc變化具有重要作用。在過去半個世紀中，溫帶落葉林物種對氣候變化表現出廣泛的物候響應(Elmoreetal.，2012；Chenetal.，2020)，但未來溫度變化如何影響關鍵的物候期事件，如植被葉子出現和衰老的時間，還存在相當大的不確定性。

3.2 Gcc、NDVI和GPP數據提取物候指標的差異

Gcc與GPP的季節變異模式相似(圖4)，即春季升高，初夏達到峰值，隨后略有下降，秋季急劇下降。Gcc時間序列數據存在不確定性，如光照變化、被分析的ROI有其他物體造成背景效應、圖像分辨率、氣象擾動等，這些不確定性影響植被SOS日期(Ahrendsetal.，2008)。通過Gcc數據確定的其他物候指標日期均與GPP提取的物候指標日期相差不大，說明Gcc能夠敏感捕捉到植被生長動態變化。生長初期與生長末期相比，數字相機對生長末期的物候指標提取比生長初期更加準確。之前研究中，通常利用遙感數據提取的植被指數估算生態系統物候指標日期，雖然遙感數據存在時間間隔長、云層遮擋等不確定性，但基于遙感數據提取的植被指數仍是物候研究中使用最廣泛的指標(Duetal.，2019；Wangetal.，2016)。本研究中，NDVI同樣具有類似的季節變化，與GPP提取的物候指標相比差異較小，生長季開始與結束日期的相對差異均小于0.1，說明NDVI數據能夠有效捕捉植被生長動態變化，與以往利用遙感數據提取的植被指數監測結果一致(Melaasetal.，2016)。Gcc與NDVI相比，季節變化趨勢相同，均為先升高，達到峰值后略有降低，隨后急劇下降。衛星植被指數與數字相機提取的物候指標比較發現，Gcc與衛星植被指數具有顯著相關性，Gcc能夠敏感捕捉生態系統植被生長變化(Wangetal.，2017)，說明SRS-NDVI測量儀測定的NDVI數據用于物候指標估算的可行性較高。

3.3 局限與展望

本研究僅針對一個生態觀測站，且時間序列數據較短，后續應繼續構建補充數據集。物候學研究需要跨植物物種和不同時空尺度的植被狀態長期觀測(Sonnentagetal.，2012)，隨著物候觀測網絡的進一步發展，更多生態觀測站將配備物候相機，因此，更多的環境變量數據和更多的全球觀測點物候圖像會顯著改善物候變化和生態系統碳循環研究(Whartonetal.，2016)。未來將收集更多站點的環境變量數據和物候圖像，并添加遙感圖像(如Phenocams、Modis、Landsat 8和Sentinel 2)的長期觀測數據，以評估自然因素對不同植被類型、不同生態系統物候轉換日期的影響，提高物候模擬精度。

4 結論

本研究以北京松山落葉闊葉林生態系統為例，利用數字相機拍攝圖像提取物候信息，結合當地微氣象觀測系統的環境因子數據，分析二者間的相關性；將數字圖像提取的Gcc、NDVI與GPP數據進行對比，分析其差異性，得到以下結論：1)在松山落葉闊葉林生態系統中，基于數字相機提取的Gcc與Ta、Ts、PAR、VPD等環境因子具有很強相關性，溫度是影響植被生長變化的關鍵因子，降雨對Gcc具有“觸發”作用；2)與GPP數據提取的物候指標相比，基于Gcc數據能夠敏感捕捉到植被生長動態變化，且生長末期物候指標捕捉精確度高于生長初期；NDVI數據也可以敏感捕捉植被生長動態變化，與GPP提取物候指標的相對差異較小。