田義超,黃遠(yuǎn)林,張 強(qiáng),陶 進(jìn),張亞麗,黃 鵠,周國清

1 北部灣大學(xué)資源與環(huán)境學(xué)院·海洋地理信息資源開發(fā)利用重點(diǎn)實(shí)驗(yàn)室, 欽州 535000 2 北部灣大學(xué)廣西北部灣海洋生物多樣性養(yǎng)護(hù)重點(diǎn)實(shí)驗(yàn)室, 欽州 535000 3 桂林理工大學(xué)廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室, 桂林 541004

全球氣候變化與陸地生態(tài)系統(tǒng)響應(yīng)是當(dāng)前全球變化研究的重要內(nèi)容,有關(guān)地表植被覆蓋與環(huán)境演變的關(guān)系是全球變化中最復(fù)雜、最具活力的研究內(nèi)容[1]。在全球生態(tài)系統(tǒng)中,植被是該系統(tǒng)的重要組成部分,植被在保持水土流失、調(diào)節(jié)大氣、維持氣候及整個(gè)生態(tài)系統(tǒng)穩(wěn)定等方面具有重要的作用[2]。植被凈初級生產(chǎn)力是表征全球氣候變化和陸地生態(tài)系統(tǒng)響應(yīng)的重要指標(biāo)[3-4],它作為陸地生態(tài)系統(tǒng)中物質(zhì)循環(huán)與能量流動的重要組成部分[5],不僅反映植物群落在自然環(huán)境和生態(tài)系統(tǒng)中的生產(chǎn)能力,也是判定生態(tài)系統(tǒng)的碳源、碳匯和評價(jià)陸地調(diào)節(jié)生態(tài)系統(tǒng)生態(tài)過程的主要因子[6]。因此,始于1987年的國際地圈—生物圈計(jì)劃[7]、全球變化與陸地生態(tài)系統(tǒng)[8]和全球碳項(xiàng)目[9]都將植被凈初級生產(chǎn)力的研究作為其核心研究內(nèi)容之一。通過對歷史時(shí)期不同地域的植被凈初級生產(chǎn)力時(shí)空變化進(jìn)行動態(tài)監(jiān)測,可以理清區(qū)域生態(tài)環(huán)境變化情況,這對于理解全球氣候變化對陸地生態(tài)系統(tǒng)植被變化過程的作用機(jī)制,掌握研究區(qū)域碳源碳匯時(shí)空變化特征,促進(jìn)區(qū)域生態(tài)建設(shè)具有重要的理論與現(xiàn)實(shí)意義。

早期在進(jìn)行區(qū)域凈初級生產(chǎn)力監(jiān)測時(shí),主要依靠傳統(tǒng)的站點(diǎn)實(shí)測方法,該方法觀測效率低,難以在大尺度上進(jìn)行推廣,因此采用模型模擬NPP成為一種可行且廣為接受的研究方法。自1876年Ebermayer開始對凈初級生產(chǎn)力的模型進(jìn)行研究之后[10],國內(nèi)外學(xué)者和專家提出了眾多的計(jì)算和觀測NPP的研究方法,并在模型模擬方面取得了一系列可喜的成果,推動了凈初級生產(chǎn)力模型及其模擬方法的進(jìn)一步發(fā)展。目前,用于估算凈初級生產(chǎn)力的模型大體上可分為氣候相關(guān)統(tǒng)計(jì)模型[11],生態(tài)系統(tǒng)生態(tài)過程模型[12-13]和光能利用模型CASA (Carnegie Ames Stanford Approach)[14]等。CASA模型是光能利用率模型的一種,該模型利用光合有效輻射和光能利用率兩個(gè)參數(shù)來估算區(qū)域的凈初級生產(chǎn)力,同時(shí)引入最適溫度、水分利用效率等光合作用脅迫因子,比較適用于區(qū)域或者大尺度上的NPP時(shí)空動態(tài)監(jiān)測與分析,甚至在全球尺度上的NPP時(shí)空動態(tài)監(jiān)測方面得到了廣泛的應(yīng)用。國內(nèi)很多學(xué)者應(yīng)用CASA模型對NPP時(shí)空變化特征及其對氣候變化的響應(yīng)進(jìn)行了定量研究,部分學(xué)者如樸世龍等[15]、朱文泉等[16]、周廣勝等[17]、穆少杰等[18]基于CASA模型或改進(jìn)的CASA模型對中國以及NPP的時(shí)空分布特征開展了大量的研究工作。從研究的空間尺度來看,國內(nèi)的NPP研究集中在全國以及東北、華北等區(qū)域,研究尺度多為省域、縣域以及鄉(xiāng)鎮(zhèn)等行政區(qū)單元[19],已有的流域尺度上凈初級生產(chǎn)力的研究成果也主要集中于北方典型流域,如黃河流域[20]、黑河流域[21]等。然而,對中國南亞熱帶獨(dú)流入海河流自20世紀(jì)初期以來凈初級生產(chǎn)力的演變規(guī)律、空間差異特征及其對氣候變化和人類活動的響應(yīng)認(rèn)識比較薄弱。

南流江流域,位于廣西壯族自治區(qū)東南部,是廣西獨(dú)流入海第一大河,發(fā)源于玉林市北流市大容山南側(cè),自北向南流經(jīng)北流、玉東、玉州、福綿、博白,是玉林的母親河。近年來,由于玉林市及其周邊市政府對流域的生態(tài)環(huán)境保護(hù)工作重視不夠,導(dǎo)致流域電廠、工業(yè)園區(qū)以及采礦場侵占流域中的自然保護(hù)區(qū)的現(xiàn)象時(shí)有發(fā)生,濫砍亂伐導(dǎo)致區(qū)域的水土流失嚴(yán)重,致使河流水質(zhì)下降,目前針對該流域的生態(tài)環(huán)境缺乏定量化評估,已有研究關(guān)于生態(tài)環(huán)境的定量化評估多以定性分析為主,而關(guān)于南亞熱帶獨(dú)流入海河流NPP時(shí)空變化規(guī)律、氣候變化對植被的主控因素以及人類活動通過改造地表土地覆被類型繼而影響到流域NPP變化等方面的研究成果較少。另外,由于南流江流域地貌類型多樣,山地丘陵地貌分布于上游、河谷平原地貌分布于中游,而下游則以濱海濕地地貌為主,在如此復(fù)雜的陸海交錯(cuò)地貌區(qū)如何定量識別氣候變化和人類活動對植被凈初級生產(chǎn)力的影響,在國際上不但缺少相關(guān)數(shù)據(jù)支持,更缺少方法技術(shù)上的經(jīng)驗(yàn)和貢獻(xiàn)。鑒于此,本研究以北部灣南流江流域?yàn)檠芯繉ο?基于光能利用率模型,利用遙感、氣象和植被類型等數(shù)據(jù)估算了南流江流域2000—2015年的凈初級生產(chǎn)力,借助于Theil-Sen趨勢、Mann-Kendall檢驗(yàn)以及Hurst指數(shù)等數(shù)理統(tǒng)計(jì)方法對研究區(qū)NPP的時(shí)空變化特征、未來趨勢及其驅(qū)動力因素進(jìn)行了定量分析。本研究不僅意義重大,同時(shí)可為南流江流域新時(shí)期水土保持工作,生態(tài)環(huán)境建設(shè)以及流域生態(tài)系統(tǒng)可持續(xù)性管理等工作的開展提供數(shù)據(jù)支撐和理論依據(jù)。

1 研究區(qū)概況

圖1 北部灣南流江流域在廣西區(qū)地理位置Fig..1 Location of the Nanliu River basin in the Beibu Gulf in Guangxi Autonomous Region

選擇北部灣南流江流域?yàn)檠芯繀^(qū),該區(qū)位于廣西壯族自治區(qū)東南部,地處北緯21°21′—23°04′,東經(jīng)105°47′—107°41′之間,東西長約139.22 km,南北長約155.06 km,面積9337.26 km2。該地區(qū)氣候?qū)倌蟻啛釒Ъ撅L(fēng)氣候,具有亞熱帶向熱帶過渡性質(zhì)的海洋季風(fēng)氣候特點(diǎn),氣候溫暖,冬短夏長,流域多年平均氣溫在21.5—22.4℃,四季適宜農(nóng)作物生長。均日照1630—1800 h,冬季偶有輕霜,無霜期長達(dá)320 d以上;多年平均降水量為1400—1760 mm。南流江流域是廣西獨(dú)流入海第一大河,發(fā)源于北流市大容山南側(cè),自北向南流,故稱南流江。流域內(nèi)地勢平坦,有玉林盆地、博白盆地和南流江三角洲,土地肥沃,農(nóng)業(yè)發(fā)達(dá)。在下游,南流江分流入海,并且在出海處形成網(wǎng)狀河系,造就了廣西最大的三角洲——南流江三角洲。三角洲地勢低平,是廣西重要的稻谷、甘蔗、花生生產(chǎn)區(qū),其外沿也多島嶼、灘涂,適宜捕撈和海產(chǎn)養(yǎng)殖。

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

本研究所選用的數(shù)據(jù)包括氣象、歸一化指數(shù)NDVI、植被覆被類型以及DEM數(shù)據(jù)。氣象數(shù)據(jù)為2000年—2015年研究區(qū)各個(gè)氣象站點(diǎn)的逐日氣溫和逐日降水量數(shù)據(jù),此數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn)。凈初級生產(chǎn)力的估算主要選擇了氣象因子數(shù)據(jù)(降水量、氣溫和太陽凈輻射量)、植被類型數(shù)據(jù)和歸一化植被指數(shù)等數(shù)據(jù)。研究區(qū)植被覆被類型數(shù)據(jù)來源于全國2000—2010年生態(tài)環(huán)境調(diào)查30 m分辨率的生態(tài)系統(tǒng)分類數(shù)據(jù);NDVI數(shù)據(jù)為研究區(qū)2000—2015年每月的MODIS13Q1遙感影像,來源于美國國家宇航局(http://modis.gsfc.nasa.gov),該數(shù)據(jù)的空間分辨率為250 m,獲取數(shù)據(jù)后利用MRT(Modis Reprojection Tools)工具、ArcGIS 10.2以及ENVI 5.1軟件對下載的遙感數(shù)據(jù)進(jìn)行投影變換、裁剪等一系列操作,所有空間數(shù)據(jù)統(tǒng)一投影為WGS84、Albers Equal Area Conic投影坐標(biāo)系統(tǒng)。1∶5萬DEM數(shù)據(jù)(http://www.gscloud.cn/)來源于中國地理空間數(shù)據(jù)云,主要用于研究區(qū)流域的劃定與提取,南流江流域矢量邊界的確定是通過ArcGIS 10.2中的Hydrology模塊進(jìn)行流域邊界自動提取,依據(jù)《廣西水資源綜合規(guī)劃報(bào)告》南流江流域的上游段以博白水文站以上的區(qū)域?yàn)榧畢^(qū),中游段和下游段的集水區(qū)則以小江水文站為劃分單元,下游段則主要是地勢比較平坦的南流江河口三角洲。

2.2 研究方法

2.2.1凈初級生產(chǎn)力估算CASA模型

凈初級生產(chǎn)力作為陸地地表碳循環(huán)過程的重要組成部分,能反映出陸地生態(tài)系統(tǒng)的質(zhì)量狀況,但是直接或全面的測量出區(qū)域或全球尺度的NPP是很困難的。所以利用模型估算已經(jīng)成為一種重要的研究方法。目前估算NPP的模型分為三類,分別為氣候生產(chǎn)力模型、生態(tài)系統(tǒng)模型和光能利用率模型。CASA(Carnegie-Ames-Stanford Approach)模型是一個(gè)基于過程的遙感模型[15],耦合了生態(tài)系統(tǒng)生產(chǎn)力和土壤碳、氮通量,由網(wǎng)格化的氣候、輻射、土壤和遙感植被指數(shù)數(shù)據(jù)集驅(qū)動,模型在計(jì)算時(shí)主要由植被所吸收的光合有效輻射(APAR)與光能轉(zhuǎn)化率兩個(gè)變量來確定。本研究在對北部灣南流江流域2000—2015年的陸地植被凈初級生產(chǎn)力進(jìn)行估算時(shí),采用的是朱文泉等改進(jìn)的CASA光能利用率模型[18]。

2.2.2Theil-Sen Median趨勢分析

該研究將Theil-Sen趨勢分析方法與Mann-Kendall 檢驗(yàn)方法[22]使用到NPP時(shí)間序列的分析當(dāng)中,可以解釋植被NPP長期的變化的顯著趨勢。

Sen趨勢度(ρ)計(jì)算公式為：

(1)

式中,xj,xi為NPP時(shí)間序列。當(dāng)ρ<0時(shí),表示時(shí)間序列NPP成下降趨勢,當(dāng)ρ>0時(shí),表示時(shí)間序列呈上升趨勢,通過Mann-Kendall方法進(jìn)行顯著性檢驗(yàn)。

Mann-Kendall檢驗(yàn)(MK檢驗(yàn))是在氣象學(xué)和水文學(xué)中,比較常用的時(shí)間序列趨勢檢驗(yàn)方法之一。MK檢驗(yàn)?zāi)軌蛱蕹贁?shù)異常值,樣本不需要遵從一定的分布,適用于非正態(tài)分布的數(shù)據(jù)。

Mann-Kendall檢驗(yàn)的公式為：

(2)

(3)

(4)

式中,Q為檢驗(yàn)統(tǒng)計(jì)量;Z為標(biāo)準(zhǔn)化后的檢驗(yàn)統(tǒng)計(jì)量;xj、xi為時(shí)間序列數(shù)據(jù);n為樣本數(shù),當(dāng)n>8時(shí),Q近似為正態(tài)分布,其均值和方差計(jì)算公式為：

E(Q)=0

(5)

(6)

標(biāo)準(zhǔn)化后Z為標(biāo)準(zhǔn)正態(tài)分布,若|Z|>Z1-ɑ/2,表示存在明顯趨勢變化。Z1-ɑ/2為標(biāo)準(zhǔn)正態(tài)分布表在置信度水平ɑ下對應(yīng)的值。本文中置信度水平ɑ為0.05,自由度為15-2=13。根據(jù)研究區(qū)15年NPP年均值數(shù)據(jù),基于Theil-Sen Median trend 和Mann-Kendall 檢驗(yàn)分析原理,借助于MATLAB軟件編程實(shí)現(xiàn)NPP的逐象元柵格計(jì)算。

2.2.3凈初級生產(chǎn)力R/S分析

英國的水文學(xué)家Hurst[23]首先提出了R/S分析方法(Rescaled Range Analysis Method,重新標(biāo)度極差分析法),后來該方法經(jīng)過Mandelbrot和Wallis等學(xué)者的進(jìn)一步修正和完善[24],很快將其發(fā)展成為研究長時(shí)間序列數(shù)據(jù)的經(jīng)典方法,目前該方法已在在水文學(xué)、經(jīng)濟(jì)學(xué)、氣候?qū)W等領(lǐng)域得到了廣泛的應(yīng)用。不同的Hurst指數(shù)H的取值范圍為H(0

(1)若0.5

(2) 若H=0.5,表明NPP時(shí)間序列為互相獨(dú)立的隨機(jī)序列;

(3) 若0

2.2.4偏相關(guān)性分析

本研究采用基于象元分析的空間分析方法[25],在matlab19.0軟件的支持下計(jì)算NPP與氣溫和降水之間的關(guān)系,相關(guān)系數(shù)的分析公式如下：

(7)

偏相關(guān)系數(shù)計(jì)算公式如下：

(8)

式中,Rxy.z為固定自變量z之后因變量x與之變量y的偏相關(guān)系數(shù),偏相關(guān)系數(shù)的顯著性檢驗(yàn)一般采用t檢驗(yàn)法[25]。

3結(jié)果與分析3.1植被凈初級生產(chǎn)力時(shí)間變化特征

2000—2015年北部灣南流江流域植被NPP總體上呈現(xiàn)出上升趨勢,增速為44.03 g C m-2(10 a)-1(圖2),明顯快于廣西自治區(qū)2000—2011年植被凈初級生產(chǎn)力(全區(qū)植被NPP的平均增長速率0.01g C m-2(10 a)-1[26]),未通過顯著性水平為0.01的檢驗(yàn)。流域NPP的變化范圍為911.30—1075.44 g C/m2,平均值為995.25 g C/m2,最大值出現(xiàn)在2011年,其值達(dá)到1075.44 g C/m2,高出平均值80.19 g C/m2,最小值則出現(xiàn)在2006年,低于平均值83.95 g C/m2,從整體上來看,流域植被凈初級生產(chǎn)力呈現(xiàn)“波動―上升”趨勢,其中2011年為NPP的變化閾值點(diǎn),2011年以前植被NPP呈現(xiàn)出上升趨勢,2011年以后則呈現(xiàn)下降趨勢。由不同流域分區(qū)凈初級生產(chǎn)力(圖2.b、c、d)可知：不同分區(qū)植被NPP的多年平均值均呈現(xiàn)出波動增加的趨勢,且NPP波動的幅度及規(guī)律與研究區(qū)整體上的NPP波動規(guī)律保持一致,但是其增長的速度存在顯著差異,具體來說,流域上游和下游地區(qū)植被NPP快于總體,而中游地區(qū)由于植被基質(zhì)較好,年內(nèi)NPP值都保持較高趨勢,增長速率慢于全區(qū)。

圖2 研究區(qū)不同區(qū)域植被NPP的空間變化率Fig.2 Spatial change rate of vegetation NPP in different regions of the study area

3.2 植被凈初級生產(chǎn)力空間變化特征

3.2.1年均值分布特征及其波動性

受南流江流域地形地貌、水汽以及植被覆蓋等自然因子空間格局的綜合影響,北部灣南流江流域植被凈初級生產(chǎn)力多年均值表現(xiàn)出明顯的地域分異規(guī)律。由圖3可知：南流江流域的中游植被NPP最高(1098.99 g C/m2),下游次之(1041.71 g C/m2),而上游最小(1013.22 g C/m2)。受到地形和地貌多種因素的共同影響,位于流域中游地區(qū)的六萬大山地帶、下游浦北縣西部及其西南部五皇山地帶是植被NPP的高值區(qū)域,由于該區(qū)域?qū)偾鹆甑孛差愋?植被資源豐富,天然植被為南亞熱帶常綠季雨林,導(dǎo)致該該區(qū)域NPP較高,其值介于1100—1233 g C/m2之間。位于流域中游博白縣的南部地區(qū)以及流域下游靈山縣的南部地區(qū)NPP次之,其值范圍為900—1100 g C/m2,而玉林市的城郊地區(qū)以及流域下游的出?？诘腘PP則處于500—700 g C/m2之間。由于玉林市轄區(qū)和南流江三角洲合浦縣受到城市化擴(kuò)張等人類活動的影響較大,導(dǎo)致區(qū)域植被退化嚴(yán)重,因此該區(qū)域植被NPP處于較低水平,其值為300—500 g C/m2之間,局部地區(qū)甚至小于300 g C/m2。綜上,流域多年NPP的平均值呈現(xiàn)出中游地區(qū)高,上游和下游地區(qū)低的空間分布格局。

2000—2015年植被凈初級生產(chǎn)力的空間波動程度如圖3所示,1)高波動性區(qū)域主要聚集于玉林市轄區(qū)及其周邊,南流江河口三角洲、合浦-欽州以及靈山的交錯(cuò)地區(qū)、浦北縣與博白縣過渡地區(qū),其中玉林市轄區(qū)周邊以及南流江河口三角洲主要受到人類活動的干擾較大,而合浦-欽州以及靈山的交錯(cuò)地區(qū),浦北縣與博白縣過渡地區(qū)這兩個(gè)片區(qū)分別為洪潮江水庫以及小江水庫,該地區(qū)水庫周邊濕地因季節(jié)性水位漲落而使周邊植被受到季節(jié)性淹沒的影響,因此凈初級生產(chǎn)力波動幅度較大;2)低波動性區(qū)域主要分布于六萬大山地帶以及西南部五皇山地帶,該區(qū)域植被主要與此區(qū)域的花崗巖性有關(guān),屬于典型的亞熱帶季風(fēng)雨林植被、原始森林等植被遍布于該區(qū)域,導(dǎo)致該區(qū)域凈初級生產(chǎn)力呈現(xiàn)出穩(wěn)定狀態(tài);3)中度波動性的NPP集中分布在合浦縣的低海拔地區(qū),其次玉林市的鄉(xiāng)郊區(qū)周邊也有零星分布,該區(qū)域植被的中等波動性主要與該區(qū)域大片的農(nóng)田植被有關(guān)?？傮w上來看,南流江流域NPP的多年分布呈現(xiàn)出明顯的空間分異規(guī)律,主要以低波動性為主,高波動性鑲嵌于低波動性周邊,兩個(gè)大型水庫周邊以及城鄉(xiāng)建設(shè)用地?cái)U(kuò)張區(qū)域處于高波動狀態(tài),而流域的中部十萬大山地帶則處于低波動狀態(tài)。因此,未來應(yīng)該適度加強(qiáng)對NPP高波動區(qū)域的保護(hù)和管理。

圖3 研究區(qū)NPP平均值及變異系數(shù)空間分布Fig.3 Spatial distribution of average NPP and its variation coefficient in study area

3.2.2NPP Sen變化趨勢及顯著性檢驗(yàn)

為了揭示北部灣南流江流域凈初級生產(chǎn)力在過去16年的變化趨勢及其變化過程,本研究基于Theil-Sen趨勢以及Mann-Kendall檢驗(yàn)方法,使用ArcGIS 10.2軟件對每年的NPP象元值進(jìn)行提取,之后借助于Matlab 2013軟件對16年的NPP進(jìn)行趨勢分析,將趨勢分析的結(jié)果與MK檢驗(yàn)的結(jié)果進(jìn)行疊加,可得到NPP增加與減少的空間分布情況,依據(jù)該方法可依次將NPP劃分為顯著減少,不顯著減少,不顯著增加以及顯著增加4種類別。由表1和圖4可知：北部灣南流江流域NPP在空間上呈現(xiàn)出增加的趨勢(增加區(qū)域的面積>減少區(qū)域的面積),其中NPP輕微減少的區(qū)域所占的比重最大,其值為46.57%,輕微增加占43.25%,明顯增加的區(qū)域?yàn)?.60%,顯著減少的區(qū)域所占比重最小,其值僅為1.58%。由圖4可知,NPP的sen趨勢值介于-77.10—74.80 g C m-2a-1之間,平均增長速率為9.29 g C m-2a-1。NPP顯著增加的區(qū)域主要集中在玉林市的南部地區(qū)、浦北縣的南部以及合浦縣的西部地區(qū);不顯著增加的區(qū)域則主要分布于上游的河谷地區(qū)、浦北縣與博白縣的交錯(cuò)地帶,另外靈山縣的東部地區(qū)也有大面積分布;不顯著減少的地區(qū)則主要分布于博白縣的六萬大山地區(qū),該地區(qū)雖然植被的基質(zhì)較好,植被凈初級生產(chǎn)力的平均值較高,但是該區(qū)域的NPP在過去16年期間則呈現(xiàn)出微度減少趨勢;顯著減少的區(qū)域則主要集中在玉林市的市轄區(qū)周邊以及合浦縣南流江流域的出?？诟浇?其中,玉林市植被凈初級生產(chǎn)力顯著減少的原因主要與該地區(qū)的城市化水平有關(guān),而南流江出?？贜PP之所以呈現(xiàn)出顯著減少趨勢究其原因則主要與該地區(qū)植被受到大范圍破壞,其植被大面積轉(zhuǎn)變?yōu)槲r塘和魚塘有關(guān)。

圖4 研究區(qū)凈初級生產(chǎn)力空間變化趨勢及顯著性檢驗(yàn)Fig.4 Spatial distribution of NPP sen trend and its Mann-Kendall test in study area

表1 研究區(qū)凈初級生產(chǎn)力變化趨勢

Sen trend：Theil-Sen趨勢;Z：顯著性水平

3.2.3植被凈初級生產(chǎn)力Hurst指數(shù)及未來變化趨勢

由圖5可知,北部灣南流江流域凈初級生產(chǎn)力Hurst的范圍為0—0.99993,平均值為0.70,植被凈初級生產(chǎn)力的反持續(xù)序列僅占到區(qū)域整體的NPP百分比為0.3986%,持續(xù)性序列占到區(qū)域整體的的NPP百分比為99.6014%。從圖6中NPP的正態(tài)分布圖中可知反持續(xù)性、弱反持續(xù)性、弱持續(xù)性和強(qiáng)持續(xù)性的面積百分比分別為0.34%、2.85%、75.73%和21.08%。北部灣南流江流域凈初級生產(chǎn)力Hurst指數(shù)的正態(tài)分布圖的空間形狀表現(xiàn)出單峰右偏趨勢,說明多年NPP的均值持續(xù)性趨勢顯著大于反持續(xù)性,由此可知南流量流域植被凈初級生產(chǎn)力NPP整體的變化趨勢則處于持續(xù)性的狀態(tài)。從凈初級生產(chǎn)力Hurst指數(shù)空間分布狀況(圖5)可以看出,南流江流域Hurst指數(shù)分布的趨勢為：弱持續(xù)性呈現(xiàn)出大面積分布狀態(tài),強(qiáng)持續(xù)性分布于流域的河谷地帶,而強(qiáng)反持續(xù)與弱反持續(xù)性則零星散布于研究區(qū)的上游以及下游地區(qū)。

為了更進(jìn)一步明晰北部灣南流江流域凈初級生產(chǎn)力NPP的變化趨勢及其未來的可持續(xù)性特征,本研究將NPP Sen趨勢的顯著性檢驗(yàn)值(圖4)與Hurst指數(shù)的計(jì)算結(jié)果(圖5)在ArcGIS 10.2中進(jìn)行疊加分析,可以得到北部灣南流江流域植被凈初級生產(chǎn)力增加、減少與未來可持續(xù)性的耦合圖(圖6)。由(圖6和表2)可以看出：在NPP空間尺度上,南流江流域未來植被凈初級生產(chǎn)力持續(xù)性增加的序列所占的百分比(51.64%)高于持續(xù)性減少的百分比(45.18%)。其中強(qiáng)持續(xù)性減少的組合所占的面積比重為5.52%,主要分布于玉林市的城鄉(xiāng)過渡地區(qū)、博白縣的東南部、靈山縣的南部,其次,合浦縣的南流江流域出?？谝灿写竺娣e分布;弱持續(xù)性減少的組合所占的百分比最大,其值為39.66%,主要分布于南流江流域中游的六萬大山地帶,另外,流域下游的五皇山地區(qū)也有大范圍分布;弱持續(xù)性增加的組合所占的百分比與弱持續(xù)性減少的比重相當(dāng),其值為36.07%,主要分布于流域上游河谷地區(qū)向山體的過渡地帶,其次流域下游的河谷地區(qū)也有少量分布;強(qiáng)持續(xù)性增加的組合所占的百分比最小,僅為15.57%,主要分布于流域的上游的河谷地帶,而強(qiáng)反持續(xù)與弱反持續(xù)性則零星散布于研究區(qū)的上游以及下游地區(qū)十萬大山的南麓地帶、浦北縣的東南部,另外,欽州-合浦-靈山的交錯(cuò)地帶也呈現(xiàn)出片狀分布狀態(tài),凈初級生產(chǎn)力Hurst在該類區(qū)域呈現(xiàn)出持續(xù)性增加的區(qū)域,預(yù)示著未來凈初級生產(chǎn)力將有持續(xù)性增加趨勢。

圖5 Hurst指數(shù)和NPP持續(xù)性特征Fig.5 Hurst and sustainability of inter-annual NPP in study area

圖6 Hurst+sen指數(shù)空間分布及指數(shù)正態(tài)分布圖Fig.6 The spatial structure of Hurst+sen index and Normal distribution plot of the Hurst index

3.3 植被凈初級生產(chǎn)力驅(qū)動因素

3.3.1氣象因子的特征及其變化趨勢

由于南流江流域可用的氣象站點(diǎn)數(shù)據(jù)非常稀少,本研究將運(yùn)用泰森多邊形對附近的雨量站進(jìn)行插值,可以得出南流江流域降水和氣溫的變化特征(圖7)及其空間分布狀況(圖8)。南流江流域2000—2015年均降水量介于1361.50—2431.16 mm之間,多年均值為1802.09 mm,其中2008年的年均降水量最大,而最小值則出現(xiàn)在2004年,2000—2015年降水量呈現(xiàn)出顯著增加趨勢(P<0.01),增加的速率為177.35 mm/10 a;南流江流域2000—2015年均氣溫介于21.69—23.17 ℃之間,多年均值為22.53 ℃,其中2015年的年均氣溫最高,而最小值則出現(xiàn)在2011年,2000—2015年氣溫呈現(xiàn)出波動不顯著下降趨勢,下降的速率為0.019 ℃/10 a。

表2 凈初級生產(chǎn)力NPP趨勢以及Hurst指數(shù)統(tǒng)計(jì)特征

根據(jù)南流江流域2000—2015年的降水量和氣溫的多年均值空間分布特征(圖8)可以看出,氣溫多年平均的空間分布情況則呈現(xiàn)出明顯的地域分異規(guī)律,氣溫較高的地區(qū)主要分布在南流江流域的南部地區(qū)以及流域的河谷地帶,而流域西北部地區(qū)的浦北縣由于地處于六萬大山和五皇山的交接地區(qū),平均海拔較高(六萬大山主峰六萬頂海拔1115 m),導(dǎo)致區(qū)域氣溫偏低。南流江流域由于受到季風(fēng)氣候的影響,降雨量從流域南部到流域中游地區(qū)呈現(xiàn)出減小趨勢,且分布不均勻,東北部地區(qū)以及流域下游南部地區(qū)之間地區(qū)的降雨較多,尤其是下游的合浦縣,降雨量最多,流域的東南部降水量偏少。

圖7 研究區(qū)年降水量和氣溫變化趨勢Fig.7 Variation trend of annual precipitation and temperature in study area

圖8 研究區(qū)多年平均降水量和氣溫空間分布Fig.8 Spatial distribution of annual precipitation and temperature in study area

3.3.2植被NPP與氣象因子之間的關(guān)系

降水和氣溫等氣象因子的變化對陸表植被的生長具有重要的作用,為了定量識別流域凈初級生產(chǎn)力對氣象因子的響應(yīng)特性,本研究逐像元計(jì)算了流域2000—2015年的植被凈初級生產(chǎn)力與年均降水量和年均氣溫之間的偏相關(guān)系數(shù)(圖9),由圖可知,凈初級生產(chǎn)力對氣溫和降水的響應(yīng)特性表現(xiàn)出顯著的空間差異特征,且植被NPP對氣溫的相關(guān)系數(shù)高于對降水的偏相關(guān)系數(shù)。由植被NPP與年均氣溫的偏相關(guān)空間特征(圖9)可知,植被NPP對氣溫的偏相關(guān)系數(shù)介于-0.87—0.95之間,平均值為0.09,由均值可知植被NPP與氣溫之間呈現(xiàn)出正相關(guān)關(guān)系,表明在南流江流域氣溫越高的地區(qū)凈初級生產(chǎn)力越高。具體來說,在北部灣南流江流域,植被NPP與氣溫的偏相關(guān)系數(shù)正負(fù)相關(guān)區(qū)域分別占流域總面積的15.91%和80.09%,其中只有1.95%的正相關(guān)區(qū)域通過了顯著性水平為0.01的檢驗(yàn),主要分布于玉林市的西北部地區(qū),而21.93%的負(fù)相關(guān)區(qū)域通過了顯著性水平為0.01的檢驗(yàn),主要分布于南流江流域的中游地區(qū),其次浦北縣的西南部也有斑塊狀分布。植被NPP對降水的偏相關(guān)系數(shù)介于-0.96—0.86之間,平均值為-0.25,由均值可知植被NPP與氣溫之間呈現(xiàn)出負(fù)相關(guān)關(guān)系,表明在南流江流域降水越高的地區(qū)凈初級生產(chǎn)力越低。具體來說,在北部灣南流江流域,植被NPP與降水的偏相關(guān)系數(shù)正負(fù)相關(guān)區(qū)域分別占流域總面積的62.56%和37.44%,其中18.26%的正相關(guān)區(qū)域通過了顯著性水平為0.01的檢驗(yàn),主要分布于浦北縣與博白縣的交接地帶,而21.93%的負(fù)相關(guān)區(qū)域通過了顯著性水平為0.01的檢驗(yàn),主要分布于浦北縣的五皇山地區(qū)。綜上,植被凈初級生產(chǎn)力NPP與氣溫呈正相關(guān)性,與降水呈負(fù)相關(guān)性的特征。

圖9 研究區(qū)NPP與氣溫和降水的偏相關(guān)系數(shù)Fig.9 The partial correlation coefficient between NPP and temperature and precipitation in study area

3.3.3植被NPP與土地利用之間的關(guān)系

流域凈初級生產(chǎn)力NPP的變化主要受到氣候變化和人類活動的影響,其中氣候變化中的降水和氣溫對植被凈初級生產(chǎn)力的變化具有重要的作用,而人類活動對凈初級生產(chǎn)力的影響則主要通過土地利用方式的改變來體現(xiàn)。本研究為了揭示土地利用變化對植被NPP的損耗與增益效應(yīng),將南流江流域2000—2015年的土地利用圖進(jìn)行空間疊置分析,在北部灣南流江流域土地利用變化圖譜中(圖10),共有36類土地利用圖譜單元,1類未發(fā)生轉(zhuǎn)化,其中南流江流域上游玉林市的北部地區(qū)以及下游的合浦縣是土地利用轉(zhuǎn)移過程最為復(fù)雜的區(qū)域,在上游中耕地轉(zhuǎn)變?yōu)榻ㄔO(shè)用地的區(qū)域所占的比重最大,而下游則以耕地變?yōu)闈竦氐仡悶橹?當(dāng)?shù)卮迕褚詽O業(yè)養(yǎng)殖為主,土地多用于蝦塘和魚塘養(yǎng)殖,分布區(qū)域主要集中在南流江流域的下游地區(qū)。將土地利用發(fā)生變化的圖譜單元按照面積遞降的順序進(jìn)行排序,在GIS中計(jì)算每個(gè)圖譜單元的變化率(表3),其中前15類的土地利用圖譜發(fā)生轉(zhuǎn)移的面積累積百分比達(dá)到總變化面積的95.31%。

圖10 研究區(qū)土地利用變化圖譜及其空間集聚特征Fig.10 Land use changes during the period of 2000—2015 in study area

從表3和圖10中可以看出,近16年來南流江流域土地利用類型面積變化的最大的圖譜單元為“耕地-建設(shè)用地”,研究期內(nèi)共有4780.62 hm2的耕地轉(zhuǎn)化為建設(shè)用地,占總變化面積的百分比為34.24%,由土地利用面積變化導(dǎo)致的NPP損失值達(dá)到4715.62 t/a,其面積分布區(qū)域主要集中在流域上游玉林市的城鄉(xiāng)過渡地帶。林地轉(zhuǎn)換為建設(shè)用地的面積所占的比重次之,占到變化面積的20.20%,由林地轉(zhuǎn)化為建設(shè)用地所導(dǎo)致的NPP損失值達(dá)到2835.19 t/a,變化圖譜主要分布靈山縣與合浦縣的交接地帶,該地區(qū)主要是受到靈山縣的擴(kuò)張,周邊林地受到人類活動的濫砍亂伐導(dǎo)致林地面積大幅度降低;耕地轉(zhuǎn)換為林地的圖譜單元位居第三位,其面積為1408.59 hm2,占所有變化面積的10.09%,該圖譜單元所導(dǎo)致的NPP損失值達(dá)到2510.03 t/a,空間分布區(qū)域上主要集中在南流江流域中游的博白縣南部的河谷地區(qū),說明退耕還林還草使得人類對于森林干擾減少,使其流域中游地區(qū)植被得到較好恢復(fù);而草地轉(zhuǎn)換為建設(shè)用地的比重最小,其值僅為1.10%,導(dǎo)致NPP減少的損失值僅為184.63 t/a。

3.3.4不同分區(qū)氣候因子及土地利用變化與NPP的關(guān)系

從表4中可以看出,在不同波動區(qū)、不同變化區(qū)以及不同持續(xù)區(qū)多年平均NPP與降水的偏相關(guān)系數(shù)都為負(fù)值,而與氣溫的相關(guān)系數(shù)均為正值,說明在北部灣地區(qū)氣候因子中的氣溫對植被的凈初級生產(chǎn)力起到控制作用。就不同波動區(qū)而言,高波動區(qū)植被對氣溫的相關(guān)系數(shù)高于降水,且高波動性區(qū)域土地利用變化的幅度顯著大于低波動性區(qū)域;從不同變化趨勢區(qū)可以看出,NPP增加的區(qū)域與降水和氣溫的相關(guān)系數(shù)的絕對值都為0.27,但是不同土地利用變化的面積核密度值在NPP增加區(qū)中小于NPP減少區(qū),說明在南流江流域NPP的增加趨勢與土地利用的穩(wěn)定性有著直接的關(guān)系,土地利用變化越劇烈,區(qū)域的NPP增長的趨勢則越慢;就NPP未來變化趨勢而言,持續(xù)性序列與降水和氣溫的相關(guān)系數(shù)都小于反持續(xù)性序列,且土地利用面積變化核密度值在持續(xù)性區(qū)域中小于反持續(xù)性區(qū)域,這與不同變化趨勢區(qū)NPP與土地利用的關(guān)系類似,流域土地利用變化越劇烈則區(qū)域的NPP的可持續(xù)性則呈現(xiàn)出降低趨勢。

表3 研究區(qū)2000—2015年土地利用變化類型變化排序及其導(dǎo)致的NPP減少值

表4 不同分區(qū)氣候因子及土地利用變化與NPP的關(guān)系

4 討論

4.1 植被凈初級生產(chǎn)力模型精度驗(yàn)證

圖11 CASA模型估算NPP與MODIS17數(shù)據(jù)的關(guān)系Fig.11 Correlation of simulated NPP based on CASA model and MODIS17 data

區(qū)域凈初級生產(chǎn)力的精度評價(jià)一直以來都是遙感學(xué)界以及生態(tài)學(xué)中的難點(diǎn)和有爭議的環(huán)節(jié),NPP估算的結(jié)果驗(yàn)證一般來說分為兩種,一種為實(shí)測值驗(yàn)證,這種方法是通過估算值與實(shí)測的NPP數(shù)值進(jìn)行對比驗(yàn)證。另外一種是相對法,即通過將模型的估算結(jié)果與其他模型結(jié)果或者其他遙感產(chǎn)品進(jìn)行對比來評價(jià)。第一種方法是通過野外獲得研究區(qū)的實(shí)測生物量,其數(shù)據(jù)比較可靠,但是目前南流江流域植被生物量的實(shí)地調(diào)查數(shù)據(jù)目前還未開展相關(guān)的研究工作。基于此,本研究采用第二種方法對CASA模型的估算結(jié)果進(jìn)行精度驗(yàn)證,所采用的精度驗(yàn)證的產(chǎn)品為MODIS 17 NPP產(chǎn)品,該產(chǎn)品是基于BIOME-BGC模型計(jì)算出的全球凈初級生產(chǎn)力數(shù)據(jù)[27],國內(nèi)也有大量的研究應(yīng)用該產(chǎn)品驗(yàn)證了CASA模型估算結(jié)果的可靠性,如謝寶妮等[28],孫慶玲等[29]和張繼平等[30]分別使用該產(chǎn)品對黃土高原、武陵山區(qū)以及三江源地區(qū)的凈初級生產(chǎn)力進(jìn)行了精度和可靠性驗(yàn)證。本研究在進(jìn)行模型驗(yàn)證時(shí),由于原始的MODIS17 NPP產(chǎn)品的空間分辨率為1000 m,本研究首先對該數(shù)據(jù)進(jìn)行重采樣到250 m,之后使用GIS的隨機(jī)采樣工具分別在MODIS17數(shù)據(jù)和本文結(jié)果上生成30%的采樣點(diǎn),以此來分析二者之間的相關(guān)關(guān)系。從圖11可以看出,本文估算結(jié)果與MODIS17數(shù)據(jù)存在著顯著的相關(guān)性(P<0.01),相關(guān)系數(shù)為0.7943,但是本研究通過CASA模型估算的結(jié)果略低于MODIS17的平均值,這可能與MODIS17采用的是16天合成的植被指數(shù)中可能存在陰天或多云等情況,使得光合有效輻射吸收系數(shù)比實(shí)際值偏大,導(dǎo)致光合有效輻射的結(jié)果偏大,進(jìn)而增大NPP的估算值。

4.2 植被凈初級生產(chǎn)力貢獻(xiàn)率識別及其未來研究方向

由于本研究中的南流江流域氣象站點(diǎn)偏少,尤其是太陽輻射站點(diǎn)較少,本文雖然采用了插值方法對氣象數(shù)據(jù)進(jìn)行了處理,但是由于該流域廣泛存在高山地貌、丘陵地貌以及河谷地貌。因此還需要發(fā)展更為精確的插值方法以提高模型估算的精度。另外,本研究雖然采用偏相關(guān)分析方法量化了氣候因子對植被NPP的響應(yīng)關(guān)系,利用土地利用變化轉(zhuǎn)移矩陣分析了流域土地利用變化導(dǎo)致的NPP增加和減少值,但是未對氣候變化和人類活動對流域的貢獻(xiàn)率和控制區(qū)域進(jìn)行定量識別,因此未來可以從兩方面進(jìn)行入手,一方面可以采用相關(guān)的數(shù)據(jù)模型,如殘差分析模型[31]建立回歸分析方程,對每年的NPP殘差序列進(jìn)行Sen趨勢度分析,以得到的趨勢值為標(biāo)準(zhǔn),如果趨勢值為正值,說明NPP象元值受到人類活動的影響,反之,則NPP象元值則受到氣候變化的影響,利用正負(fù)值所占的比例可以計(jì)算氣候變化和人類活動對植被NPP的貢獻(xiàn)率。另一方面,可以通過統(tǒng)計(jì)研究區(qū)牲口的數(shù)量、草地、坡耕地的面積在研究區(qū)植被恢復(fù)方面的貢獻(xiàn)比例。因此,關(guān)于氣候變化和人類活動對凈初級生產(chǎn)力的貢獻(xiàn)率將是本研究以后重點(diǎn)的研究方向。

5 結(jié)論

本文以北部灣南流江流域?yàn)檠芯繉ο?基于光能利用率模型(CASA),利用遙感、氣象和植被等數(shù)據(jù)估算了研究區(qū)2000—2015年流域的凈初級生產(chǎn)力,借助于Theil-Sen趨勢、Mann-Kendall檢驗(yàn)以及Hurst指數(shù)等數(shù)理統(tǒng)計(jì)方法對研究區(qū)NPP的時(shí)空變化特征、未來趨勢及其驅(qū)動因素進(jìn)行了定量化分析。研究結(jié)果表明：

(1)時(shí)間尺度上,2000—2015年北部灣南流江流域植被凈初級生產(chǎn)力NPP總體上呈現(xiàn)出波動上升趨勢,增速為44.03 g C m-2(10 a)-1,明顯快于廣西自治區(qū)2000—2011年植被凈初級生產(chǎn)力,未通過顯著性水平為0.01的檢驗(yàn)。流域上游和下游地區(qū)植被NPP快于全區(qū),而中游地區(qū)慢于全區(qū),其中,流域上游、中游和下游的凈初級生產(chǎn)力年平均增加速度分別為：54.17g C m-2(10 a)-1、31.74g C m-2(10 a)-1和46.1g C m-2(10 a)-1。多年平均凈初級生產(chǎn)力介于911.30—1075.44 g C/m2之間,平均值為995.25 g C/m2,最大值出現(xiàn)在2011年,其值為1075.44 g C/m2,高出平均值80.19 g C/m2,最小值則出現(xiàn)在2006年,低于平均值83.95g C/m2。(2)空間尺度上,流域植被凈初級生產(chǎn)力的分布呈現(xiàn)出明顯的地域分異規(guī)律,中游植被凈初級生產(chǎn)力最高(1098.99 g C/m2),下游次之(1041.71 g C/m2),而上游最小(1013.22 g C/m2)。NPP的Sen趨勢度介于-77.10—74.80 g C m-2a-1之間,流域凈初級生產(chǎn)力在空間上呈現(xiàn)出增加的趨勢,其中NPP顯著減少的區(qū)域所占1.58%,輕微減少占46.57%,輕微增加占43.25%,而明顯增加占8.60%。

(3)空間波動性上,流域凈初級生產(chǎn)力的變異系數(shù)較大,其值介于0.01—0.70。流域NPP的多年平均值呈現(xiàn)出明顯的空間分異規(guī)律,主要以低波動性為主,其中,洪潮江水庫、小江水庫周邊以及玉林市的城鄉(xiāng)建設(shè)用地?cái)U(kuò)張區(qū)域處于高波動狀態(tài),而流域的中部六萬大山以及五皇山地帶則處于低波動狀態(tài)。

(4)未來變化趨勢上,流域凈初級生產(chǎn)力Hurst的范圍為0—0.99,平均值為0.70,植被凈初級生產(chǎn)力的反持續(xù)序列僅占到區(qū)域整體的NPP百分比為0.3986%,持續(xù)性序列占到區(qū)域整體的的NPP百分比為99.6014%,Hurst指數(shù)正態(tài)分布圖呈現(xiàn)單峰右偏分布,即NPP均值的持續(xù)性序列顯著大于反持續(xù)性序列,預(yù)示著流域NPP未來處于持續(xù)增加的趨勢。

(5)驅(qū)動機(jī)制上,流域NPP與多年平均氣溫呈正相關(guān)關(guān)系,與年均降水量呈負(fù)相關(guān)關(guān)系,表明溫度是影響該流域植被NPP的主要?dú)夂蛞蜃?。由耕地轉(zhuǎn)化為建設(shè)用地,導(dǎo)致的NPP損失值最大,其值達(dá)到4715.62 t/a,林地轉(zhuǎn)換為建設(shè)用地所導(dǎo)致的NPP損失值次之,其值為2835.19 t/a,而草地轉(zhuǎn)換為建設(shè)用地導(dǎo)致NPP損失值最小,其值僅為184.63 t/a。