張玉紅,閆 浩
(1.哈爾濱師范大學地理科學學院,黑龍江 哈爾濱 150025;2.自然資源部第二大地測量隊,黑龍江 哈爾濱 150025)
森林是地球生態系統的主體,在陸地生態系統中擁有最大的面積和最豐富的物質資源,約占陸地總面積的30%[1]。凈初級生產力(Net primary productivity,簡稱NPP)是綠色植物在單位時間和面積內累積的有機物[2],是表征地球表層植被覆蓋生物量的重要指標,它反映了森林固定大氣CO2的能力,對全球氣候研究中起到重要的參考價值[3]。隨著遙感技術(RS)與地理信息系統技術(GIS)的應用與發展,衛星遙感數據被大量應用到林業科學中來。通過遙感與地理信息系統技術,各類觀測數據更容易獲取,空間分析與模擬更為便捷準確,同時估測范圍更加廣泛。例如利用不同遙感平臺的衛星影像進行森林火災的監測與管理[4-6],或者利用地理信息系統強大的空間分析功能進行各類森林火災的模擬研究[7]。作為森林的監測與管理的重要指標,NPP的計算和模擬可以通過遙感技術解決傳統定點估算NPP數值的局限性和長期性等問題,使NPP的估算獲得更多的思路和方法[8]。同時通過遙感技術與不同的模型方法相結合,可以對森林植被的特征參數進行精確反演,獲取更多的信息數據[9]。主要的NPP 估算模型分為統計模型、過程模型和參數模型三類[10]。統計模型主要考慮氣候因子的影響,估算植被的NPP 主要考慮降水、氣溫、太陽輻射等因子與植被NPP的統計回歸與相關性。例如Miami 模型[11]、Thornthwaite Memorial 模型[12]和Chikugo 模型[13],它們輸入參數簡單,基礎數據較容易取得,能夠方便研究,可以對陸地生態系統凈初級生產力進行估算和預測[14]。但相對的能夠作為依據的生態理論缺少嚴密性,資料的準確性也有待提高。過程模型又被稱機理模型,代表模型如TEM 模型[15]、CERTURY 模型[16]、BEPS 模型[17]。這些模型從一個整體系統的角度,對物質和能量的交換進行模擬并建立了相應的模型或模型庫,從機理上模擬植被的光合作用和蒸騰蒸發等過程,從而推進了大空間尺度的植被NPP 研究。但是為了保證嚴謹的特點,過程模型的構造和機理相當繁瑣,需要的參數如植物生態特征、土壤植被、氣象變化數據等,數據需求大并且難以一一獲取,在空間上較難定量化。應用較為廣泛的則是參數模型,即光能利用率模型,其理論基礎是植被光能利用率,將植被累積的生物量認定為是植被將太陽入射輻射吸收和轉化等一系列光合作用過程得到的。這類模型的代表如CASA模型[18]、GLOPEM 模型[19],C-FIX 模型[20]。與統計模型和參數模型相比,這類模型相對簡單便利,它是通過使用遙感數據來計算光合有效輻射。因為容易獲得周期性強和范圍廣泛的遙感數據,使得這類模型的時空分辨率較好,能夠進行大至全球尺度范圍的估測與觀測。這些優勢促進了這類模型的發展研究,成為現階段NPP研究中的主要研究模型。
作為全球變化的敏感區域之一,大興安嶺地區位于北半球的中高緯度,受到了世界各國學者的關注。大興安嶺地區是重要的林業基地,也是我國唯一擁有大面積興安落葉松(Larix gmelinii)的原始林區。該區高達70%以上的森林覆蓋度和巨大的森林蓄積量,使大興安嶺森林地區有著重要的生態作用和社會效益[21-22]。1987 年5 月6 日,大興安嶺地區發生了嚴重的森林火災。這場特大火災使得森林資源遭受到了嚴重損失,植被覆蓋以及森林功能都需要較長時間的恢復。文中基于遙感數據并結合森林調查數據對大興安嶺地區圖強林業局的多年森林NPP的數量和空間分布進行估測,并分析影響森林NPP的主要自然因素以及對森林火災之后的NPP 分布及恢復情況進行遙感監測,可以為林區管理提供基礎信息,為林業規劃和持續發展提供科學的依據。
黑龍江省大興安嶺地區圖強林業局位于122°18'E-123°28'E,52°15'N-53°33'N 之間(圖1),屬于漠河縣管轄范圍,是我國最北的木材生產基地。整個林業局的總面積約為5 000 km2,其中林業用地約為4 970 km2,而森林覆蓋度達到87%。該研究區屬于寒溫帶大陸性氣候,空氣濕度小,氣溫年較差與月較差均較大。適應這種氣候的植被種類較少,大約70%以上的樹木都為興安落葉松(Larix gmelinii),其余部分又以白樺(Betula platyphylla)和樟子松(Pinus sylvestris var.mongolica)為主,分別占整個林業局森林蓄積量的21%和7.5%,其他樹種總共不足0.5%。該區地勢較為平坦,但也有低矮丘陵,多處于東部,整體地勢成南高北低的分布。1987 年5 月份大興安嶺地區發生了特大火災,總過火面積超過10 000 km2,有林面積達到了70%左右。其中,圖強林業局的森林資源受到了嚴重損失,過火區域超過一半以上。

圖1 大興安嶺地區圖強林業局的位置示意圖及過火強度圖Fig.1 Location and fire intensity map of Tuqiang Forestry Bureau,Daxing’anling Region
1.2.1 遙感數據
文中所使用的遙感衛星數據均來自美國Landsat 系列衛星,軌道號為122/23。根據NDVI 估值與植被生長季的關系[23],并考慮到7 月份植被的生長進入穩定階段,因此選擇時相為1987 年到2013 年中7 月份的遙感影像,云量皆為5%以下。經過影像數據的預處理之后,獲取相關植被指數數據,并利用決策樹進行植被信息提取。
1.2.2 氣象數據
氣象數據來源于英國東英吉利大學氣候研究所(CRU)提供的1987年到2013年的逐月網格氣溫和降雨數據,以及附近氣象站點記錄的月平均溫度、月總降水量和太陽輻射數據。
1972 年,Monteith 在光能利用率的原理上提出了利用光合有效輻射的觀測值和光能利用率來估算植被生產力的方法,即建立了Monteith 方程[24]。在此基礎上,Potter 等學者相繼提出了計算全球植被凈初級生產力模型并建立了CASA 模型(Carnegie-Ames-Stanford Approach),完成了基于光能利用率方法的陸地凈初級生產力全球估算,是目前對NPP估算研究中使用較多的模型[25]。該模型的主要優點在于:模型相對簡單,需要的參數不多,因而便于計算和處理;由于部分參數來源于遙感數據,因而不受地面實測站點的條件限制,可以進行大區域上的NPP 估測;同時在遙感數據的基礎上可以對植被進行分類來反映植被的變化狀況,能夠借用多時相遙感數據來對NPP進行動態監測等[26-27]。
CASA模型公式如下:

式中,APAR(x,t)代表的是t月中像元x吸收的光合有效輻射(gC·m-2),ε(x,t)代表的是t月中像元x的實際光能利用率(gC·M-1J)。
2.1.1 APAR的估算
通過遙感數據來估算植物葉片所吸收的光合有效輻射(APAR),其公式為:

式中,SOL(x,t)所代表的是像x處在t月的太陽總輻射量,FPAR(x,t)所代表的是植被層能夠吸收的入射光合有效輻射所占比例,常數0.5所代表的是植被能夠使用的太陽有效輻射所占的比例。在一定范圍內,可以利用FPAR 與植被NDVI來確定FPAR 的最大值和最小值。同時也可以使用FPAR 與比值植被指數之間的線性關系來進行估算。前者估算的數值較高于實測值,后者則低于實測值但誤差相對較小。因此將兩個結果結合起來可以提高FPAR的估算精度。
2.1.2 光能利用率的估算
光能利用率一般是指植被在地表單位面積內,通過光合作用所產生的有機物中的化學潛能,與同一時間內該塊面積上收到的有效太陽輻射之比。氣溫、土壤濕度以及大氣濕度、氣壓等環境因子會通過影響植物的光合能力,進而調節植被的NPP數值。所以,在模型中光能利用率是估算NPP的重要因子。其中,最大光能利用率主要受溫度和水分含量的影響,公式如下:

式中:Tε1(x,t)反映在低溫和高溫時植物內在的生化作用對光合作用的限制,從而降低凈第一性生產力。Tε2(x,t)表示環境溫度從最適溫度Topt(x,t)向高溫或低溫變化時植物光能利用率逐漸變小的趨勢。水分脅迫影響系數Wε(x,t)反映了植物所能利用的有效水分條件對光能利用率的影響。εmax是理想狀態下的最大光能利用率(gC·MJ-1)[28-29]。因此,歸納其模型參數主要如表1所示。

表1 模型參數簡表Table 1 Model parameter profile
CASA 模型參數主要包括固定參數、植被類型、NDVI、月平均溫度、月降水總量以及月太陽總輻射。NDVI通過預處理后的遙感數據在ENVI軟件平臺上計算得到;植被類型圖是基于遙感影像并使用決策樹方法獲得;氣溫和降水數據均來自于CRU數據,通過空間插值處理對氣象數據進行網格內插,得到研究區整體的氣象柵格數據;太陽輻射數據則來自于氣象站點和歷史數據統計得到。通過以上數據的處理獲得1987年至2013年共13期模型參數數據。
3.1.1 圖強林業局森林NPP的數量變化
將以上模型參數數據輸入到模型中,獲得13期7月份森林NPP估算結果。由于文中主要研究森林所帶來的NPP產量,所以將非森林地區的數值刻意減少讓它近于0,但這同時也降低了研究區全部整體的NPP平均值。因此在統計過程中,除了統計全部區域的NPP 平均值外,也統計了森林區域的NPP 平均值。圖強林業局地區的1987-2013年7月的NPP產量統計值如圖2所示。

圖2 7月份森林NPP變化趨勢Fig.2 Trend of forest NPP Change in July
從NPP 數值上看,1987 年7 月的總產值、平均值和最大值由于森林火災都遠遠低于其他年份。各類統計值在2010 年達到最大,2005 年達到最小值。其他年份數據均在207 gC/m2~209 gC/m2范圍之內,數值變動很小,平均值為208.03 gC/m2。從整體趨勢上看,研究區7 月份森林NPP 的產值呈上升趨勢。在火災發生之后的第二年即1988 年7 月的數值就直接提升了近乎一倍。在之后的20 年里,數值增加相對較少。1993、1994、1995 三年7 月的NPP 產量平均值為4.817×1011gC,2000、2001、2002 三年7 月的NPP 產量平均值4.975×1011gC,NPP 產值增加了0.16×1011gC;2005年和2006年7月份的NPP 產值下降幅度很大,到了2008年之后開始回升;2009、2010、2011 三年7 月的NPP 產量平均值為5.711×1011gC,與之前相比NPP 產值增加了0.74×1011gC。2013年各統計值略有下降。
從以上結果可知,火災后森林NPP 恢復很快,尤其是最大值。但是平均值的恢復就需要一定的時間。同時也可以看出,除了突發的火災影響之外,森林NPP的變化還受到其它因素的影響。
3.1.2 圖強林業局森林NPP的空間分布變化
文中選取圖強林業局1987年-2013年之間幾個典型的7月森林NPP空間分布圖,包括1987年、1988年、1995年、2002年、2005年、2006年、2009年和2013年。
(1)火災后恢復初期(1987-1988年)
從1987 年7 月的NPP 產量值分布圖清晰的看出圖強林業局受到森林大火的影響,在森林大火過后,整個研究區北部受到了嚴重的損失,很多地區的森林NPP產量值低于10 gC/m2,輕度過火區域的數值也在平均NPP產量值之下。未受火災的地區仍有很多高NPP產值地區,但總值仍較往年減少很多(圖3(a))。1988年研究區進入了初步恢復階段,過火區域的NPP產量得到了顯著提升,但與未過火區域相比,平均NPP產值仍相對較低。由于部分樹木生長快速或恢復密集,在過火區域也出現了高NPP產值(圖3(b))。

圖3 圖強林業局7月份NPP分布Fig.3 NPP distribution map of Tuqiang Forestry Bureau in July
(2)火災后恢復中期(1988-2009年)
1995 年后數值呈緩慢上升趨勢,森林的NPP 產值進入了平穩恢復階段。這一期間,森林的恢復不像1988年一樣迅速,仍可以從圖中看出火燒痕跡,北部有較多地區數值偏低。但對比1988年的數據可以看出,與未過火區域相連的區域或輕度過火區域恢復進展良好,已經出現大面積的NPP高產值地區。以1988年最先恢復的地區為中心向外的區域的數值也達到了標準值以上(圖3(c))。到2002年,7月份森林NPP 的產值相比7 年前有了整體的上升,數值較低的區域也減少很多,但仍能看出北部數值較南部數值偏低,重度過火區域仍然達不到森林植被NPP 產值的平均值。相比之下最北部的區域已經基本恢復,高數值區域也較1995 年有所增多(圖3(d))。2005 年森林NPP 各統計值突然降低,無論南北均出現大量的低數值區域。2006年較2005年相比,雖然有部分區域正常,但也可以清晰地看出多處數值偏低區域。從歷史資料來看,這兩年并未發生重大事件。從分布圖上看,區域整體數值變低,并非是局部區域減少。結合現有資料,初步判斷是這兩年7月的氣候較往年不同而影響了森林植被的NPP產值。從圖中還可以看出,雖然整體數值偏低,但是南北的差異已經并不明顯(圖3(e),圖3(f))。
(3)火災后恢復后期(2009-2013年)
到了2009 年,森林的NPP 產值進入到穩定階段,之后的3 年時間里7 月的森林NPP 產值上下浮動不大,且整體數值都明顯高于10 年前。從分布上很難看到曾經的過火區域,南北分布趨勢也呈現一致,高產值區域也大面積出現(圖3(g))。因此判斷從2009 年開始,森林已經近乎恢復到原有的生態環境。2013 年的數據由于基礎遙感影像有云覆蓋的原因,導致南部底端部分森林沒有數值,大大降低了NPP的數值,但當月的NPP最大值浮動并不大,總產值與前幾年相比略有下降,但仍然保持在平均值以上(圖3(h))。
從整個研究時段來看,研究區北部地區受到大火影響,多年總體數值一直較低,且變化明顯呈增長趨勢,到2009 年后則趨近平穩;南部未過火區域沒有明顯變化,多年來數值多由當月氣候決定。整體上來看,非森林地區的周邊數值變動較為明顯。結合2011 年和2013 年的分布圖進行對比,可以看出非森林周邊多處的NPP 數值明顯降低,距離遠處數值則相對偏高。結合當地考察與專家咨詢得知,圖強林業局作為重要的木材儲蓄地,要對森林木材進行采集,道路兩邊的樹木方便采集,這會間接導致NPP 產值降低;采集過后會撒下樹種,重新出現新生樹木,而新生樹木地區的森林覆蓋度相對較低,使得該區域的NPP數值偏低。
森林NPP 估算所利用的CASA 模型中考慮到的因子主要有氣溫、降水、太陽輻射、樹種。文中選取的數據估算的均為7月份的森林NPP值,考慮到當地每年7月的太陽總輻射和溫度變化都不大,因此文中重點分析降水和不同樹種對森林NPP的影響。
3.2.1 降水對森林NPP的影響
圖4為1987年至2013年期間研究區的溫度和降水變化,結合前面的森林NPP估算結果可以發現幾個特殊年份。如1993 年和2002 年的降水量遠高于鄰近年份,但NPP 值卻十分相似。相反2005 年,2006 年,2008年降水相對較少,NPP 數值也較往年減少很多。為了更清楚的研究降水對NPP 生產量的影響,本文制作控制實驗數據,即保留1987年、1993年、2002年、2005年、2006年和2008年的地面遙感數據、溫度、太陽總輻射的數據,使用2009年的氣象降雨數據,再根據模型估算出在此情況下的研究區NPP的生產量情況,與實際的估測值進行對比。

圖4 1987-2013年7月氣溫降雨變化圖Fig.4 The variation of temperature and total precipitation in July,1987-2013
結果表明,1993年和2002年結果與控制實驗結果相比,雖然控制實驗結果各項數值都有所減少,但變化不大,總產值、最大值與平均值均減少的數值都不到0.8 gC/m2。因而證明并不是降水的數值越多NPP 的數值就會越大。若水分過多超出植被和土壤能夠吸收利用的數值,還有可能抑制植被NPP 的產值。1987 年、2005 年、2006 年、2008 年的控制實驗結果數據則表明:隨著月總降雨量的增加,NPP 的產值得到了大量的提升。
結合2005 年實際估測值與控制數據的區域分布圖(圖5)來看,其實在2005 年,就已經看不到曾經的過火區域,南北分布趨勢也呈現一致,與2009年以后的分布圖相似。判斷影響2006年和2008年NPP產值的降低的主要原因為水分過低,沒能達到較適宜的降水量。

圖5 2005年實際模擬分布圖與控制降水數據分布圖對比Fig.5 Comparison of actual simulated distribution and control data in 2005
3.2.2 不同樹種對NPP 產值的影響
植物的特性影響其光能利用率,如樹葉的形狀與大小、樹干的高度、對陽光的喜好、水分的利用吸收等。一般情況下,闊葉林的光能利用率要高于針葉林。因此在模型中,不同的樹種有著不同數值的最大光能利用率,光能利用率較高的樹種的NPP產值會相對較高。文中研究區林種主要有白樺、落葉松和樟子松。表2為不同樹種的NPP產值統計結果。

表2 研究區不同樹種NPP產值數據Table 2 NPP value of different tree species in the study area
白樺、落葉松和樟子松的NPP數值增長趨勢與森林整體的增長趨勢保持一致。同樣都是在1988年大幅度增長后逐漸緩慢,在2005年和2006年數值減少,到2008年數值增加并再度進入平穩。
白樺是喜光的樹種,因為它生命力強、生長周期快等特征,使它在森林火燒后能夠最先生長起來,形成大片的新生白樺林。白樺的產值為圖強林業局總產值的一半,最多比例為51.41%,最少為49.17%。7 月產值最多的在2010 年為2.863×1011gC,平均值從1988 年的144.16 gC/m2提升到了2010 年的166.97 gC/m2,比落葉松的同月的NPP 平均值高出36.92 gC/m2,其主要原因是白樺樹屬于闊葉林,落葉松和樟子松屬于針葉林,相比這2 種樹種,白樺樹本身的葉面積指數,植被冠層覆蓋度都相對較高,白樺樹葉內吸收的入射光合有效輻射以及光能利用率也要高于落葉松和樟子松,因此單位面積的NPP產量要高于落葉松和樟子松。落葉松也是喜光的樹種,適應性強、耐寒,它的樹干與樹高一般高于白樺樹,且全年生長。落葉松的NPP 多年均值大約占總產值的34.38%左右,落葉松最多比例為35.83%。樟子松也是適應力強的樹種,但由于它樹冠稀疏,針葉稀少,再加上樟子松在圖強林業局的占地面積較少僅有7.5%,因此樟子松的總產值與其他兩類樹種相比數值相差很多,NPP 的平均值也比其它兩類樹種低很多,多年均值只有95.58 gC/m2,只在2009 年、2010年、2011 年3 年超過了100 gC/m2。從NPP 的總產值比例來看,相比于其它年份的樟子松比例大多數都低于10%,1987 年的15.33%就顯得格外突出,這是由于1987 年落葉松大面積被燒毀,從而導致樟子松的NPP 產值比例相對上升。
可見,不同的樹種對NPP 的數值有很大影響。在圖強林業局,7 月的主要NPP 產量來源為白樺樹,約占全部產值的50.13%,落葉松和樟子松約則分別約占34.38%和10.27%。
研究區在1987 年的特大火災中過火面積超過一半,其中絕大部分是重度過火區域,且均位于圖強林業局總局較遠的北部地區,能夠投入的人力物力與南部地區相比較少,無法在種植樹苗之后進行長期的人工撫育措施,栽種的新生樹苗的生長環境得不到人為的改善,導致過火后北部的森林NPP 產值沒能得到顯著提升。
從森林NPP 總產值上來看,1987 年7 月的過火地區NPP 產值僅占當月總產值的3.97%,第二年7 月的產值大幅度提升,占當月總產值的21.92%,其主要原因是在火災過后,林業局采取了人工造林的恢復措施,在過火地區大面積種植落葉松的幼苗,并采取人工更新等技術加快森林的恢復。1993 年到1995 年期間,過火區域的總產值并沒有較大的提升,所占比例也略微有所下降。這是由于圖強林業局只有6、7、8 三個月適宜植被的生長,其余月份溫度都較低并且幾乎無降雨,在這樣的氣候環境下,幼苗的成活率會大大降低,生長速度也會變得較為緩慢。2000年開始,過火地區的NPP產值上升趨勢就變得較為明顯。白樺和落葉松本來就是適應當地氣候的樹種,存活下來的幼樹會逐漸進入高生長旺盛期,此時即使沒有人工的撫育樹木也能夠正常的生長,因此過火地區的森林開始的自然恢復程度較好,產值也逐年提升。2005 年由于降雨量過低的原因,過火地區的總產值略微有所下降,但是占總產值的比例仍然提升到了25.02%,此時過火區域與非過火區域的差異已經逐漸減少。2009年以后占總產值的比例一直在26%以上。
結合研究區7 月份的過火森林區域NPP 平均值和全部森林區域NPP 平均值來看,可以明顯看出其動態變化(圖6)。1987 年二者相差較大,這一差距在1988 年大幅度縮減,到1995 年開始差距逐漸縮小,到2009年以后就趨近一致,說明過火區域的森林植被基本恢復。

圖6 過火區和全部森林區域的NPP平均值對比Fig.6 Comparison of NPP mean value between the burned area and total forest area
文中以黑龍江省大興安嶺地區圖強林業局為研究區,充分利用GIS 和RS 技術,結合生態模型對火災后的森林植被恢復情況進行模擬和監測。從森林植被凈初級生產力NPP 入手,對研究區1988 年以來至2013年的森林恢復情況進行遙感估算,并通過控制實驗對比分析影響森林NPP估算模型的主要因素和過火區域的植被恢復情況。研究表明:1987年森林大火造成的研究區森林NPP數值下降嚴重。但是在適當的氣候條件下,森林NPP 逐年恢復。1988 年進行了森林人工種植,之后逐年恢復;1993 年之后過火區域的NPP 產值逐漸增加,到2009年與未過火區域的NPP產值接近,植被覆蓋恢復到平均水平,恢復時間大約在20年左右。從空間分布來看,研究區北部地區受火災影響較大,多年森林NPP總體數值變化明顯呈增長趨勢,到2009年后則趨近平穩;南部未過火區域沒有明顯變化。在降水量適中的情況下,森林NPP產值呈現增加的趨勢,恢復較快;同時不同的樹種的恢復程度有所差別。
可以看出,結合參數生態模型的遙感監測解決了傳統定點估算NPP 數值的局限性與長期性等問題,同時能夠更有效地計算光合有效輻射,使得大尺度的估測與監測得以實現。但是作為一種參數模型,CASA 生態模型存在一定的不足之處,例如缺乏生物學基礎,對植物、大氣與土壤之間的聯系的解釋尚不完整等。但是與傳統的統計模型相比,利用容易獲得的周期性較強且大范圍的遙感數據來計算光合有效輻射,其時空分辨率更好,更容易進行大尺度的NPP 估測與觀測,是目前NPP 研究的主要模型。文中由于遙感數據和其它資料的獲取限制,缺乏研究時段內的逐年和逐月的森林NPP產值估算;同時,研究區在1987年特大火災之后原有的火燒跡地又發生過多起火災,對森林的恢復程度與恢復時間均產生了影響,因此不能全面分析而精準分析出森林的恢復情況;今后應嘗試其他渠道獲取研究區的遙感數據和環境數據,多方面多因素綜合考慮以提高結果的準確性;同時要充分利用實測數據對模型估算的結果進行精度驗證。