陳佩佩, 張利華, 崔越, 陳俊宏

(中國(guó)地質(zhì)大學(xué)(武漢)地理與信息工程學(xué)院,武漢 430074)

0 引言

植被凈初級(jí)生產(chǎn)力(net primary productivity, NPP)是指在單位時(shí)間、單位面積內(nèi),植被通過(guò)光合作用得到的全部有機(jī)物去除植物呼吸消耗能量后的剩余部分[1],NPP不僅在碳循環(huán)中發(fā)揮著無(wú)可代替的作用,同時(shí)在生態(tài)系統(tǒng)中也扮演著不可或缺的角色[2],很多國(guó)際研究項(xiàng)目,如“生物圈計(jì)劃”、“世界氣候研究計(jì)劃”都把它列為研究的中心內(nèi)容[3]。目前,在環(huán)境問(wèn)題以及碳循環(huán)備受關(guān)注的同時(shí),研究植被凈初級(jí)生產(chǎn)力以及其影響因子對(duì)減緩全球變暖以及分析區(qū)域生態(tài)系統(tǒng)循環(huán)過(guò)程等研究具有重要意義。

獲取NPP的方法有很多: 站點(diǎn)實(shí)測(cè)法、遙感數(shù)據(jù)產(chǎn)品以及模型模擬等。建立實(shí)測(cè)站點(diǎn)受研究區(qū)域條件影響較大并且只適用于小范圍研究區(qū)域,人力耗費(fèi)較大。MODIS 的NPP產(chǎn)品(MOD17A3)目前是廣泛被認(rèn)可和應(yīng)用的數(shù)據(jù)產(chǎn)品,但是此產(chǎn)品時(shí)間分辨率為1 a,不適合研究較小時(shí)間尺度上的NPP情況,且空間分辨率為500 m×500 m,數(shù)據(jù)缺失像元較多,對(duì)于小流域尺度的NPP空間異質(zhì)性分析適用性不強(qiáng)。隨著遙感技術(shù)的不斷發(fā)展,遙感資料覆蓋范圍廣,更新速度快,遙感數(shù)據(jù)結(jié)合生態(tài)模型進(jìn)行NPP的估算成為當(dāng)前的研究熱點(diǎn)之一,例如BEPS-Terrainlab模型[4]、InIEC模型[5]、BIOME-BGC[6]模型等。其中BEPS-Terrainlab模型是BEPS生態(tài)模型[7]和Terrainlab分布式水文模型[8]的耦合,是遙感驅(qū)動(dòng)的生態(tài)過(guò)程模型,在實(shí)踐運(yùn)用的過(guò)程中形成了不同的版本[9]。BEPS-Terrainlab V2.0在模擬水文循環(huán)過(guò)程時(shí)不僅考慮了地表覆蓋、植被結(jié)構(gòu)、土壤特性等地理因素,同時(shí)還引入太陽(yáng)輻射、降水等氣候因素對(duì)結(jié)果的影響; 另外考慮到地形因素以及陰陽(yáng)葉對(duì)生態(tài)水文循環(huán)的影響,輸出數(shù)據(jù)的空間分辨率為30 m×30 m,以日為時(shí)間步長(zhǎng),對(duì)于小流域尺度生態(tài)水文參數(shù)(如NPP,ET等)的時(shí)空演變分析更具有適用性。

犟河流域是南水北調(diào)水源保護(hù)的核心區(qū)域,本文在進(jìn)行模型模擬精度分析的基礎(chǔ)上,基于BEPS-Terrainlab V2.0模型對(duì)流域1986—2017年NPP的時(shí)空分布特征進(jìn)行模擬,獲得流域高分辨率的NPP分布數(shù)據(jù),并探討NPP時(shí)空演變特征及其對(duì)地形、土地覆蓋、氣候變化等因子的環(huán)境響應(yīng),模擬結(jié)果對(duì)流域植被保護(hù)、生態(tài)恢復(fù)以及碳循環(huán)特征等均具有重要的借鑒意義。

1 研究區(qū)概況

犟河是堵河的一級(jí)支流,河流總長(zhǎng)50 km,發(fā)源于張家灣區(qū)的花果街道,共有大小支流22條,于黃龍鎮(zhèn)東灣村處與堵河進(jìn)行交匯,隨堵河一起流入丹江口水庫(kù)。犟河流域(圖1)位于湖北省十堰市城區(qū)西部,地理坐標(biāo)為E110°33′50″～110°42′9″,N32°27′27″～32°42′60″,流域面積326 km2,整體地勢(shì)西南高、東北低,流域海拔范圍為139～1 441 m。流域位于亞熱帶季風(fēng)區(qū),降雨以每年6—9月為主,年降雨量800～1 000 m。植被覆蓋度較高,森林覆蓋范圍較廣,主要植被類型有暖溫帶落葉闊葉林、暖溫帶常綠針葉林、針闊混交林等。

2 數(shù)據(jù)與方法

2.1 BEPS-TerrainlabV2.0模型原理

BEPS-TerrainlabV2.0包含冠層輻射傳輸、碳循環(huán)、水循環(huán)、生理調(diào)節(jié)等模塊[10],模型輸入數(shù)據(jù)包括高程、坡度坡向、土壤質(zhì)地、氣象、土地覆蓋類型、葉面積指數(shù)(leaf area index, LAI),空間分辨率為30 m×30 m。模型輸出數(shù)據(jù)包括蒸散發(fā)、總初級(jí)生產(chǎn)力、凈初級(jí)生產(chǎn)力、地表徑流、地下水位、土壤溫度和土壤濕度等生態(tài)水文數(shù)據(jù),以日為時(shí)間步長(zhǎng),空間分辨率與輸入數(shù)據(jù)一致。其在模擬植物光合作用時(shí)將植物葉片劃分為陽(yáng)葉和陰葉[11],分別計(jì)算陽(yáng)葉和陰葉的光合過(guò)程,利用逐日積分方法在時(shí)間上進(jìn)行了擴(kuò)展,獲得植被的總初級(jí)生產(chǎn)力,在扣除植物自身呼吸之后,獲得植物的凈初級(jí)生產(chǎn)力[12]。

1)將冠層分為陽(yáng)葉與陰葉分別計(jì)算植被冠層葉面積指數(shù),公式為:

(1)

,

(2)

式中:LAI為葉面積指數(shù);LAIsun和LAIshade分別為陽(yáng)葉、陰葉的植被冠層葉面積指數(shù);θ為太陽(yáng)高度角; Ω為葉片聚集度指數(shù)。

2)計(jì)算植被冠層總初級(jí)生產(chǎn)力GPPcanopy,公式為:

,

(3)

式中:Asun和Ashade分別為陽(yáng)葉、陰葉的凈同化速率( mol·m-2·s-1)。

3)計(jì)算NPP,NPP通過(guò)GPP扣除植物自養(yǎng)呼吸(Ra)后得到,自養(yǎng)呼吸由生長(zhǎng)呼吸(Rg)和維持性呼吸(Rm)2部分組成,公式為:

,

(4)

,

(5)

,

(6)

,

(7)

式中:GPP為總初級(jí)生產(chǎn)力;i是指植被的不同組分(1,2,3分別代表葉、莖、根)。

2.2 數(shù)據(jù)源及其預(yù)處理

1)高程數(shù)據(jù)。數(shù)據(jù)源自地理空間數(shù)據(jù)云(http: //www.gscloud.cn/)的SRTM1 DEM數(shù)據(jù),空間分辨率為30 m,下載影像行列號(hào)為path110/row032。使用ArcGIS軟件對(duì)DEM數(shù)據(jù)進(jìn)行處理,通過(guò)水文分析等工具可以獲得研究流域邊界,同時(shí)可以獲取研究區(qū)域坡度、坡向數(shù)據(jù)以及裁剪其他數(shù)據(jù)的掩模。

2)遙感數(shù)據(jù)。數(shù)據(jù)源自美國(guó)地質(zhì)調(diào)查局官網(wǎng)(http: //glovis.usgs.gov/)的Landsat5-8影像數(shù)據(jù),選取生長(zhǎng)季(5—9月)影像進(jìn)行下載,空間分辨率為30 m×30 m,下載影像行列號(hào)為path125/row037,2013年前使用Landsat5數(shù)據(jù),2013年后使用Landsat8數(shù)據(jù)。在ENVI中對(duì)下載的遙感影像數(shù)據(jù)進(jìn)行預(yù)處理后使用流域邊界數(shù)據(jù)進(jìn)行批量裁剪,使用神經(jīng)網(wǎng)絡(luò)法進(jìn)行土地覆蓋類型分類,將土地覆蓋類型分為闊葉林、水體、城鎮(zhèn)、針葉林、混交林、裸地、農(nóng)田7種,以5 a為間隔進(jìn)行土地覆蓋分類。模型需要輸入時(shí)間序列LAI,建立常用植被指數(shù)(歸一化植被指數(shù)、比值植被指數(shù)、差值植被指數(shù)、調(diào)整土壤亮度植被指數(shù)、大氣阻抗植被指數(shù))與野外實(shí)測(cè)LAI值的回歸模型,選取了線性方程中R2最高的反演模型,即RVI與LAI的線性擬合方程,公式為:

LAI=1.073 5RVI-5.321 9

。

(8)

在ENVI中使用波段運(yùn)算得到RVI,再通過(guò)此方程計(jì)算LAI 。

3)氣象數(shù)據(jù)。獲取到的十堰市氣象數(shù)據(jù)包括1986—2017逐年間日尺度的最高溫、最低溫、平均氣溫以及降水。使用Magnus-Tentens近似法將缺少的露點(diǎn)溫度補(bǔ)充完整。研究區(qū)域所屬地級(jí)市無(wú)太陽(yáng)輻射數(shù)據(jù),使用反距離加權(quán)法將研究區(qū)域周圍4個(gè)方向的城市氣象站(安康、侯馬、南陽(yáng)、宜昌)的太陽(yáng)輻射數(shù)據(jù)進(jìn)行加權(quán)平均,得到的數(shù)據(jù)來(lái)代替研究區(qū)域的輻射數(shù)據(jù)輸入模型。

4)土壤數(shù)據(jù)。野外采集土壤樣品,通過(guò)馬爾文3000高分辨激光粒度儀測(cè)量土壤粒度,在ArcGIS中通過(guò)克里金插值繪制研究區(qū)域內(nèi)每種粒度的分布圖,在ENVI中通過(guò)波段運(yùn)算將砂土、壤質(zhì)砂土、砂質(zhì)壤土和粉砂壤土4種不同土壤質(zhì)地像元分別賦值為10,20,30和40。

2.3 研究方法

2.3.1 一元線性回歸分析法

本研究采用一元線性回歸分析法逐像元對(duì)犟河流域植被凈初級(jí)生產(chǎn)力在1986—2017年間變化趨勢(shì)進(jìn)行分析,計(jì)算公式[13]為:

,

(9)

式中:n為時(shí)間序列的年份數(shù),應(yīng)用到本研究中n值取32;i為時(shí)間變量,i=1,2,…,n;Ni指第i年的NPP均值;slope為研究因子的變化趨勢(shì),slope>0時(shí)說(shuō)明在研究時(shí)間范圍內(nèi)NPP值呈增長(zhǎng)趨勢(shì),slope<0時(shí)NPP值呈減少趨勢(shì),slope數(shù)值的大小表示NPP值增長(zhǎng)或減少的幅度。

2.3.2 M-K突變檢驗(yàn)

M-K檢測(cè)可判定時(shí)間序列中有無(wú)突變,若有,則判定突變發(fā)生的時(shí)間[14]。通過(guò)定義UF,UB這2個(gè)統(tǒng)計(jì)量來(lái)分析樣本時(shí)間序列的變化趨勢(shì),UF>0說(shuō)明變化趨勢(shì)持續(xù)增長(zhǎng),UF與UB相交于置信區(qū)間 [-1.96,1.96]內(nèi),說(shuō)明該交點(diǎn)對(duì)應(yīng)年份的相關(guān)數(shù)據(jù)在變化趨勢(shì)上產(chǎn)生突變情況。

2.3.3 線性回歸分析

線性回歸的過(guò)程就是把各個(gè)自變量和因變量的值建立回歸方程后,通過(guò)逐步迭代與擬合,找出回歸方程式中的各個(gè)系數(shù),構(gòu)造出一個(gè)函數(shù)式,使用線性回歸參數(shù)判斷自變量和因變量的關(guān)系,由于本研究自變量較多,故而使用偏相關(guān)系數(shù)來(lái)進(jìn)行相關(guān)性分析,能夠在一定程度上減少自變量直接的相關(guān)性對(duì)分析結(jié)果的影響。

3 結(jié)果與分析

3.1 模型精度分析

3.1.1 通量站數(shù)據(jù)驗(yàn)證

由于研究區(qū)域范圍內(nèi)沒(méi)有建立通量站,于是本次研究獲取大九湖盆地的渦度相關(guān)系統(tǒng)NPP通量數(shù)據(jù),獲取的數(shù)據(jù)時(shí)間范圍是2016年1月1日—2016年11月21日。大九湖盆地地理位置處于湖北省神農(nóng)架林區(qū)西部,距離研究區(qū)130余km,通過(guò)臨近地區(qū)通量站的數(shù)據(jù)來(lái)進(jìn)行模型模擬精度的檢驗(yàn),繪制實(shí)測(cè)值和模擬值的變化趨勢(shì)圖見(jiàn)圖2。在圖2中可以看到模擬值和實(shí)測(cè)值變化趨勢(shì)大致相同,且吳宗釩在大九湖區(qū)域采用BEPS-TerrainlabV2.0模型模擬NPP值與大九湖實(shí)測(cè)NPP值擬合R2=0.85[15],可以說(shuō)明該模型在模擬小流域生態(tài)水文過(guò)程上具有適用性,模擬結(jié)果具有分析價(jià)值。

圖2 犟河流域模擬NPP與大九湖實(shí)測(cè)NPP數(shù)據(jù)對(duì)比圖

3.1.2 NPP產(chǎn)品數(shù)據(jù)交叉驗(yàn)證

MODIS的MOD17A3是NPP產(chǎn)品(圖3),時(shí)間分辨率為1 a,空間分辨率為500 m, MODIS數(shù)據(jù)在該研究范圍內(nèi)存在較多空白值,且像元數(shù)量較少(共有1 471個(gè)有值像元),因而MODIS數(shù)據(jù)在小空間尺度研究上存在局限性。將模型模擬的不同覆蓋類型的年總NPP與其他模型模擬的NPP值進(jìn)行對(duì)比,結(jié)果如表1。從表中可以看出,模型模擬得到的值略高,但不同植被覆蓋類型NPP值的大小關(guān)系基本相同。

表1 BEPS-Terrainlab V2.0模型模擬土地覆蓋類型NPP數(shù)據(jù)與其他模擬或?qū)崪y(cè)值對(duì)比

圖3 MOD17A3數(shù)據(jù)空間分布

3.2 時(shí)間變化規(guī)律

3.2.1 時(shí)間分析

模型運(yùn)行輸出后得到的是1986—2017年間日尺度的NPP數(shù)據(jù),依據(jù)不同的時(shí)間單位,可以從年尺度和日尺度上2個(gè)方面上對(duì)NPP進(jìn)行時(shí)間變化的分析。以日為時(shí)間尺度,用來(lái)描述NPP在一年內(nèi)的變化,將具體年份的年內(nèi)變化進(jìn)行匯總得到NPP年內(nèi)變化的多年平均值,見(jiàn)圖4。

圖4 犟河流域NPP年內(nèi)多年平均變化

可以發(fā)現(xiàn)NPP年內(nèi)多年平均隨時(shí)間變化的走向呈單峰式變化,先升后降。冬季(1—60天、330—365天)值最低,春季(60—150天)呈上升的趨勢(shì),夏季(150—240天)時(shí)達(dá)到了峰值,也是NPP數(shù)值最高的階段,秋季(240—330天)NPP呈緩慢下降的趨勢(shì)。從NPP變化的趨勢(shì)可以看出,在夏季(150—240天)內(nèi)有突然下降的波動(dòng),并且在秋季(240—330天)存在短暫波動(dòng)較大的情況。

以年為時(shí)間尺度,統(tǒng)計(jì)1986—2017多年年總NPP,反映多年來(lái)研究區(qū)域NPP變化的情況,見(jiàn)圖5。

圖5 犟河流域NPP多年年總變化

由圖5可知,年總NPP大體變化趨勢(shì)為波動(dòng)上升,最低值出現(xiàn)在1995年,NPP值為911.86 gC·m-2·a-1,最高值出現(xiàn)在2007年,值為1 482.28 gC·m-2·a-1,1986—2017年間NPP的多年平均值為1 232.70 gC·m-2·a-1。

3.2.2 突變檢驗(yàn)

使用DSP數(shù)據(jù)處理系統(tǒng)對(duì)多年年總NPP進(jìn)行M-K檢驗(yàn),結(jié)果如圖6(a)所示,Uf-UB曲線在置信水平區(qū)間內(nèi)相交于2004年,說(shuō)明NPP時(shí)間序列在99%置信水平下突變于2004年,以2004年為分界點(diǎn)繪制多年年際NPP變化趨勢(shì)圖,如圖6(b),可以看到2004年之后NPP的變化幅度明顯大于2004年以前,2004年以前,NPP總體變化趨勢(shì)較為平緩,2004年以后NPP變化趨勢(shì)明顯增大。

(a) M-K檢驗(yàn)結(jié)果 (b) NPP變化趨勢(shì)

3.3 空間變化規(guī)律

在ENVI中將模型輸出的1986—2017年的年總NPP數(shù)據(jù)進(jìn)行波段合成,輸出研究區(qū)域多年平均NPP影像,能夠得到研究區(qū)域多年平均NPP在空間上的分布,如圖7(a)所示,通過(guò)一元線性回歸分析對(duì)1986—2017年NPP數(shù)據(jù)進(jìn)行分析,得到這32 a間NPP在空間上的變化趨勢(shì)值,如圖7(b)所示。

(a) 總NPP空間分布圖 (b) 空間變化趨勢(shì)圖

多年NPP在空間上分布呈現(xiàn)南高北低的分布,南部最高可達(dá)到1 925 gC·m-2·a-1,在主流域附近NPP的值最低,受地勢(shì)及不同土地覆蓋類型的影響,NPP的分布在空間上呈現(xiàn)差異性,圖7(b)為1986—2017年NPP值在空間上變化的趨勢(shì),在大部分區(qū)域NPP值都呈現(xiàn)增加的趨勢(shì),在河流區(qū)域變化不大。

3.4 環(huán)境影響因子分析

3.4.1 土地覆蓋因子

研究區(qū)域土地覆蓋類型分為7類,主要有針葉林、闊葉林、混交林、水體、農(nóng)田、裸地和城市。由于裸地以及水體植被覆蓋度較低且分布面積占比不是很大,于是就其他5類類型數(shù)據(jù)來(lái)分析NPP的分布情況。由圖8(a)可見(jiàn),不同植被覆蓋類型的NPP不同: 闊葉>混交>針葉>農(nóng)田>城市,結(jié)果與薛梅在東北地區(qū)估算的不同植被類型的NPP大小排序結(jié)果相同[19]。其中混交林NPP的值最高,其次針葉林NPP值略高于闊葉林。

(a) 土地覆蓋因子 (b) 土壤因子 (c) 高程因子

3.4.2 土壤因子

將研究區(qū)NPP分布按照不同土壤質(zhì)地進(jìn)行劃分,繪制箱線圖如圖8(b),犟河流域不同土壤類型平均NPP按大小排序?yàn)? 砂土(1 363.8 gC·m-2·a-1)>砂質(zhì)壤土(1 177.1 gC·m-2·a-1)>壤質(zhì)砂土(1 348.6 gC·m-2·a-1)>粉砂壤土(932.9 gC·m-2·a-1)。

3.4.3 地形因子

將研究區(qū)NPP分布按照不同高程進(jìn)行劃分如圖8(c),得到年總NPP隨高程每增加100 m的變化趨勢(shì),NPP隨著海拔的上升而增加,在100～900 m呈上升趨勢(shì),在高程為900 m時(shí)達(dá)到最大值而后變化趨勢(shì)逐漸趨于平穩(wěn),這與植被的垂直分布有關(guān),在高程為900～1 300 m時(shí),土地覆蓋類型為針葉林和混交林,根據(jù)前文中已知混交林和針葉林的NPP值較高于其他土地覆蓋類型。

3.4.4 氣象因子

影響NPP的變化的氣象因素有很多,諸如溫度、輻射以及降水等,采用線性回歸分析分析來(lái)探討NPP與這些氣象因子的相關(guān)性,可以進(jìn)行較長(zhǎng)尺度下的氣象影響因素的分析。在SPSS中添加因變量年內(nèi)NPP以及年際NPP數(shù)據(jù)和自變量氣溫、降水、輻射、風(fēng)速、LAI等自變量數(shù)據(jù),對(duì)這些數(shù)據(jù)進(jìn)行線性回歸分析,由于自變量間可能存在互相影響,探討各自變量對(duì)于因變量的影響力大小時(shí),選擇以偏相關(guān)系數(shù)大小為判斷依據(jù),以標(biāo)準(zhǔn)化回歸系數(shù)作為參考(表2)。

表2 日尺度NPP線性回歸方程參數(shù)

由表2可以看出,影響力最高的是太陽(yáng)輻射,偏相關(guān)系數(shù)為0.647,P值小于0.001,說(shuō)明輻射與NPP顯著性正相關(guān); 其次溫度的偏相關(guān)系數(shù)為0.637,P值小于0.001,具有顯著正相關(guān),說(shuō)明溫度對(duì)NPP影響也較大。降水與風(fēng)速的P值較大,與NPP存在弱負(fù)相關(guān)的關(guān)系,對(duì)NPP的影響較小(表3)。

表3 年尺度NPP線性回歸方程參數(shù)

由表3可知,偏相關(guān)系數(shù)最高的是LAI,P值小于0.001,屬于極顯著正相關(guān),降水偏相關(guān)系數(shù)為0.305,影響力有所提升; 溫度對(duì)年尺度NPP的影響顯著降低,其他因素對(duì)年尺度NPP的影響都較小,偏相關(guān)系數(shù)較低且P值較大。

4 結(jié)論

本研究基于遙感以及氣象等數(shù)據(jù)結(jié)合BEPS-TerrainlabV2.0模型對(duì)湖北省十堰市犟河流域1986—2017年NPP進(jìn)行分析,討論了其時(shí)空分布以及與環(huán)境因子的相關(guān)性分析,結(jié)果表明:

1)與相近地區(qū)通量站數(shù)據(jù)以及MODIS數(shù)據(jù)進(jìn)行分析,驗(yàn)證了該模型對(duì)于較小空間尺度的適用性,同時(shí)與其他數(shù)據(jù)產(chǎn)品進(jìn)行交叉驗(yàn)證,證實(shí)了模型結(jié)果的可用性。

2)犟河流域NPP的年際變化具有波動(dòng)性,空間上總體分布南高北低,不同植被類型的NPP模擬值為闊葉>混交>針葉>農(nóng)田>城市; 不同土壤類型NPP值為砂土>砂質(zhì)壤土>壤質(zhì)砂土>粉砂壤土,同時(shí)NPP值隨著海拔的上升而增加,在海拔為900 m以上時(shí),增加的趨勢(shì)逐漸趨于平緩。

3)日尺度上,氣象因素對(duì)NPP影響較大,太陽(yáng)輻射和溫度與NPP的偏相關(guān)系數(shù)達(dá)到0.6左右,較長(zhǎng)時(shí)間尺度上NPP對(duì)LAI的反響最強(qiáng)。

同時(shí)由于選取區(qū)域較小,無(wú)通量站觀測(cè)數(shù)據(jù),且目前國(guó)內(nèi)存在的NPP數(shù)據(jù)產(chǎn)品大多是500 m分辨率,在對(duì)模型結(jié)果精度驗(yàn)證上有一定的局限性,考慮在后續(xù)的實(shí)驗(yàn)中進(jìn)一步改進(jìn)模型結(jié)果的驗(yàn)證準(zhǔn)確性。