田义超,黄远林,张 强,陶 进,张亚丽,黄 鹄,周国清

1 北部湾大学资源与环境学院·海洋地理信息资源开发利用重点实验室, 钦州 535000 2 北部湾大学广西北部湾海洋生物多样性养护重点实验室, 钦州 535000 3 桂林理工大学广西空间信息与测绘重点实验室, 桂林 541004

全球气候变化与陆地生态系统响应是当前全球变化研究的重要内容,有关地表植被覆盖与环境演变的关系是全球变化中最复杂、最具活力的研究内容[1]。在全球生态系统中,植被是该系统的重要组成部分,植被在保持水土流失、调节大气、维持气候及整个生态系统稳定等方面具有重要的作用[2]。植被净初级生产力是表征全球气候变化和陆地生态系统响应的重要指标[3-4],它作为陆地生态系统中物质循环与能量流动的重要组成部分[5],不仅反映植物群落在自然环境和生态系统中的生产能力,也是判定生态系统的碳源、碳汇和评价陆地调节生态系统生态过程的主要因子[6]。因此,始于1987年的国际地圈—生物圈计划[7]、全球变化与陆地生态系统[8]和全球碳项目[9]都将植被净初级生产力的研究作为其核心研究内容之一。通过对历史时期不同地域的植被净初级生产力时空变化进行动态监测,可以理清区域生态环境变化情况,这对于理解全球气候变化对陆地生态系统植被变化过程的作用机制,掌握研究区域碳源碳汇时空变化特征,促进区域生态建设具有重要的理论与现实意义。

早期在进行区域净初级生产力监测时,主要依靠传统的站点实测方法,该方法观测效率低,难以在大尺度上进行推广,因此采用模型模拟NPP成为一种可行且广为接受的研究方法。自1876年Ebermayer开始对净初级生产力的模型进行研究之后[10],国内外学者和专家提出了众多的计算和观测NPP的研究方法,并在模型模拟方面取得了一系列可喜的成果,推动了净初级生产力模型及其模拟方法的进一步发展。目前,用于估算净初级生产力的模型大体上可分为气候相关统计模型[11],生态系统生态过程模型[12-13]和光能利用模型CASA (Carnegie Ames Stanford Approach)[14]等。CASA模型是光能利用率模型的一种,该模型利用光合有效辐射和光能利用率两个参数来估算区域的净初级生产力,同时引入最适温度、水分利用效率等光合作用胁迫因子,比较适用于区域或者大尺度上的NPP时空动态监测与分析,甚至在全球尺度上的NPP时空动态监测方面得到了广泛的应用。国内很多学者应用CASA模型对NPP时空变化特征及其对气候变化的响应进行了定量研究,部分学者如朴世龙等[15]、朱文泉等[16]、周广胜等[17]、穆少杰等[18]基于CASA模型或改进的CASA模型对中国以及NPP的时空分布特征开展了大量的研究工作。从研究的空间尺度来看,国内的NPP研究集中在全国以及东北、华北等区域,研究尺度多为省域、县域以及乡镇等行政区单元[19],已有的流域尺度上净初级生产力的研究成果也主要集中于北方典型流域,如黄河流域[20]、黑河流域[21]等。然而,对中国南亚热带独流入海河流自20世纪初期以来净初级生产力的演变规律、空间差异特征及其对气候变化和人类活动的响应认识比较薄弱。

南流江流域,位于广西壮族自治区东南部,是广西独流入海第一大河,发源于玉林市北流市大容山南侧,自北向南流经北流、玉东、玉州、福绵、博白,是玉林的母亲河。近年来,由于玉林市及其周边市政府对流域的生态环境保护工作重视不够,导致流域电厂、工业园区以及采矿场侵占流域中的自然保护区的现象时有发生,滥砍乱伐导致区域的水土流失严重,致使河流水质下降,目前针对该流域的生态环境缺乏定量化评估,已有研究关于生态环境的定量化评估多以定性分析为主,而关于南亚热带独流入海河流NPP时空变化规律、气候变化对植被的主控因素以及人类活动通过改造地表土地覆被类型继而影响到流域NPP变化等方面的研究成果较少。另外,由于南流江流域地貌类型多样,山地丘陵地貌分布于上游、河谷平原地貌分布于中游,而下游则以滨海湿地地貌为主,在如此复杂的陆海交错地貌区如何定量识别气候变化和人类活动对植被净初级生产力的影响,在国际上不但缺少相关数据支持,更缺少方法技术上的经验和贡献。鉴于此,本研究以北部湾南流江流域为研究对象,基于光能利用率模型,利用遥感、气象和植被类型等数据估算了南流江流域2000—2015年的净初级生产力,借助于Theil-Sen趋势、Mann-Kendall检验以及Hurst指数等数理统计方法对研究区NPP的时空变化特征、未来趋势及其驱动力因素进行了定量分析。本研究不仅意义重大,同时可为南流江流域新时期水土保持工作,生态环境建设以及流域生态系统可持续性管理等工作的开展提供数据支撑和理论依据。

1 研究区概况

图1 北部湾南流江流域在广西区地理位置Fig..1 Location of the Nanliu River basin in the Beibu Gulf in Guangxi Autonomous Region

选择北部湾南流江流域为研究区,该区位于广西壮族自治区东南部,地处北纬21°21′—23°04′,东经105°47′—107°41′之间,东西长约139.22 km,南北长约155.06 km,面积9337.26 km2。该地区气候属南亚热带季风气候,具有亚热带向热带过渡性质的海洋季风气候特点,气候温暖,冬短夏长,流域多年平均气温在21.5—22.4℃,四季适宜农作物生长。均日照1630—1800 h,冬季偶有轻霜,无霜期长达320 d以上;多年平均降水量为1400—1760 mm。南流江流域是广西独流入海第一大河,发源于北流市大容山南侧,自北向南流,故称南流江。流域内地势平坦,有玉林盆地、博白盆地和南流江三角洲,土地肥沃,农业发达。在下游,南流江分流入海,并且在出海处形成网状河系,造就了广西最大的三角洲——南流江三角洲。三角洲地势低平,是广西重要的稻谷、甘蔗、花生生产区,其外沿也多岛屿、滩涂,适宜捕捞和海产养殖。

2 数据来源与研究方法

2.1 数据来源

本研究所选用的数据包括气象、归一化指数NDVI、植被覆被类型以及DEM数据。气象数据为2000年—2015年研究区各个气象站点的逐日气温和逐日降水量数据,此数据来源于中国气象科学数据共享服务网(http://cdc.cma.gov.cn)。净初级生产力的估算主要选择了气象因子数据(降水量、气温和太阳净辐射量)、植被类型数据和归一化植被指数等数据。研究区植被覆被类型数据来源于全国2000—2010年生态环境调查30 m分辨率的生态系统分类数据;NDVI数据为研究区2000—2015年每月的MODIS13Q1遥感影像,来源于美国国家宇航局(http://modis.gsfc.nasa.gov),该数据的空间分辨率为250 m,获取数据后利用MRT(Modis Reprojection Tools)工具、ArcGIS 10.2以及ENVI 5.1软件对下载的遥感数据进行投影变换、裁剪等一系列操作,所有空间数据统一投影为WGS84、Albers Equal Area Conic投影坐标系统。1∶5万DEM数据(http://www.gscloud.cn/)来源于中国地理空间数据云,主要用于研究区流域的划定与提取,南流江流域矢量边界的确定是通过ArcGIS 10.2中的Hydrology模块进行流域边界自动提取,依据《广西水资源综合规划报告》南流江流域的上游段以博白水文站以上的区域为集水区,中游段和下游段的集水区则以小江水文站为划分单元,下游段则主要是地势比较平坦的南流江河口三角洲。

2.2 研究方法

2.2.1净初级生产力估算CASA模型

净初级生产力作为陆地地表碳循环过程的重要组成部分,能反映出陆地生态系统的质量状况,但是直接或全面的测量出区域或全球尺度的NPP是很困难的。所以利用模型估算已经成为一种重要的研究方法。目前估算NPP的模型分为三类,分别为气候生产力模型、生态系统模型和光能利用率模型。CASA(Carnegie-Ames-Stanford Approach)模型是一个基于过程的遥感模型[15],耦合了生态系统生产力和土壤碳、氮通量,由网格化的气候、辐射、土壤和遥感植被指数数据集驱动,模型在计算时主要由植被所吸收的光合有效辐射(APAR)与光能转化率两个变量来确定。本研究在对北部湾南流江流域2000—2015年的陆地植被净初级生产力进行估算时,采用的是朱文泉等改进的CASA光能利用率模型[18]。

2.2.2Theil-Sen Median趋势分析

该研究将Theil-Sen趋势分析方法与Mann-Kendall 检验方法[22]使用到NPP时间序列的分析当中,可以解释植被NPP长期的变化的显著趋势。

Sen趋势度(ρ)计算公式为：

(1)

式中,xj,xi为NPP时间序列。当ρ<0时,表示时间序列NPP成下降趋势,当ρ>0时,表示时间序列呈上升趋势,通过Mann-Kendall方法进行显著性检验。

Mann-Kendall检验(MK检验)是在气象学和水文学中,比较常用的时间序列趋势检验方法之一。MK检验能够剔除少数异常值,样本不需要遵从一定的分布,适用于非正态分布的数据。

Mann-Kendall检验的公式为：

(2)

(3)

(4)

式中,Q为检验统计量;Z为标准化后的检验统计量;xj、xi为时间序列数据;n为样本数,当n>8时,Q近似为正态分布,其均值和方差计算公式为：

E(Q)=0

(5)

(6)

标准化后Z为标准正态分布,若|Z|>Z1-ɑ/2,表示存在明显趋势变化。Z1-ɑ/2为标准正态分布表在置信度水平ɑ下对应的值。本文中置信度水平ɑ为0.05,自由度为15-2=13。根据研究区15年NPP年均值数据,基于Theil-Sen Median trend 和Mann-Kendall 检验分析原理,借助于MATLAB软件编程实现NPP的逐象元栅格计算。

2.2.3净初级生产力R/S分析

英国的水文学家Hurst[23]首先提出了R/S分析方法(Rescaled Range Analysis Method,重新标度极差分析法),后来该方法经过Mandelbrot和Wallis等学者的进一步修正和完善[24],很快将其发展成为研究长时间序列数据的经典方法,目前该方法已在在水文学、经济学、气候学等领域得到了广泛的应用。不同的Hurst指数H的取值范围为H(0

(1)若0.5

(2) 若H=0.5,表明NPP时间序列为互相独立的随机序列;

(3) 若0

2.2.4偏相关性分析

本研究采用基于象元分析的空间分析方法[25],在matlab19.0软件的支持下计算NPP与气温和降水之间的关系,相关系数的分析公式如下：

(7)

偏相关系数计算公式如下：

(8)

式中,Rxy.z为固定自变量z之后因变量x与之变量y的偏相关系数,偏相关系数的显著性检验一般采用t检验法[25]。

3结果与分析3.1植被净初级生产力时间变化特征

2000—2015年北部湾南流江流域植被NPP总体上呈现出上升趋势,增速为44.03 g C m-2(10 a)-1(图2),明显快于广西自治区2000—2011年植被净初级生产力(全区植被NPP的平均增长速率0.01g C m-2(10 a)-1[26]),未通过显著性水平为0.01的检验。流域NPP的变化范围为911.30—1075.44 g C/m2,平均值为995.25 g C/m2,最大值出现在2011年,其值达到1075.44 g C/m2,高出平均值80.19 g C/m2,最小值则出现在2006年,低于平均值83.95 g C/m2,从整体上来看,流域植被净初级生产力呈现“波动―上升”趋势,其中2011年为NPP的变化阈值点,2011年以前植被NPP呈现出上升趋势,2011年以后则呈现下降趋势。由不同流域分区净初级生产力(图2.b、c、d)可知：不同分区植被NPP的多年平均值均呈现出波动增加的趋势,且NPP波动的幅度及规律与研究区整体上的NPP波动规律保持一致,但是其增长的速度存在显著差异,具体来说,流域上游和下游地区植被NPP快于总体,而中游地区由于植被基质较好,年内NPP值都保持较高趋势,增长速率慢于全区。

图2 研究区不同区域植被NPP的空间变化率Fig.2 Spatial change rate of vegetation NPP in different regions of the study area

3.2 植被净初级生产力空间变化特征

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

受南流江流域地形地貌、水汽以及植被覆盖等自然因子空间格局的综合影响,北部湾南流江流域植被净初级生产力多年均值表现出明显的地域分异规律。由图3可知：南流江流域的中游植被NPP最高(1098.99 g C/m2),下游次之(1041.71 g C/m2),而上游最小(1013.22 g C/m2)。受到地形和地貌多种因素的共同影响,位于流域中游地区的六万大山地带、下游浦北县西部及其西南部五皇山地带是植被NPP的高值区域,由于该区域属丘陵地貌类型,植被资源丰富,天然植被为南亚热带常绿季雨林,导致该该区域NPP较高,其值介于1100—1233 g C/m2之间。位于流域中游博白县的南部地区以及流域下游灵山县的南部地区NPP次之,其值范围为900—1100 g C/m2,而玉林市的城郊地区以及流域下游的出海口的NPP则处于500—700 g C/m2之间。由于玉林市辖区和南流江三角洲合浦县受到城市化扩张等人类活动的影响较大,导致区域植被退化严重,因此该区域植被NPP处于较低水平,其值为300—500 g C/m2之间,局部地区甚至小于300 g C/m2。综上,流域多年NPP的平均值呈现出中游地区高,上游和下游地区低的空间分布格局。

2000—2015年植被净初级生产力的空间波动程度如图3所示,1)高波动性区域主要聚集于玉林市辖区及其周边,南流江河口三角洲、合浦-钦州以及灵山的交错地区、浦北县与博白县过渡地区,其中玉林市辖区周边以及南流江河口三角洲主要受到人类活动的干扰较大,而合浦-钦州以及灵山的交错地区,浦北县与博白县过渡地区这两个片区分别为洪潮江水库以及小江水库,该地区水库周边湿地因季节性水位涨落而使周边植被受到季节性淹没的影响,因此净初级生产力波动幅度较大;2)低波动性区域主要分布于六万大山地带以及西南部五皇山地带,该区域植被主要与此区域的花岗岩性有关,属于典型的亚热带季风雨林植被、原始森林等植被遍布于该区域,导致该区域净初级生产力呈现出稳定状态;3)中度波动性的NPP集中分布在合浦县的低海拔地区,其次玉林市的乡郊区周边也有零星分布,该区域植被的中等波动性主要与该区域大片的农田植被有关。总体上来看,南流江流域NPP的多年分布呈现出明显的空间分异规律,主要以低波动性为主,高波动性镶嵌于低波动性周边,两个大型水库周边以及城乡建设用地扩张区域处于高波动状态,而流域的中部十万大山地带则处于低波动状态。因此,未来应该适度加强对NPP高波动区域的保护和管理。

图3 研究区NPP平均值及变异系数空间分布Fig.3 Spatial distribution of average NPP and its variation coefficient in study area

3.2.2NPP Sen变化趋势及显著性检验

为了揭示北部湾南流江流域净初级生产力在过去16年的变化趋势及其变化过程,本研究基于Theil-Sen趋势以及Mann-Kendall检验方法,使用ArcGIS 10.2软件对每年的NPP象元值进行提取,之后借助于Matlab 2013软件对16年的NPP进行趋势分析,将趋势分析的结果与MK检验的结果进行叠加,可得到NPP增加与减少的空间分布情况,依据该方法可依次将NPP划分为显著减少,不显著减少,不显著增加以及显著增加4种类别。由表1和图4可知：北部湾南流江流域NPP在空间上呈现出增加的趋势(增加区域的面积>减少区域的面积),其中NPP轻微减少的区域所占的比重最大,其值为46.57%,轻微增加占43.25%,明显增加的区域为8.60%,显著减少的区域所占比重最小,其值仅为1.58%。由图4可知,NPP的sen趋势值介于-77.10—74.80 g C m-2a-1之间,平均增长速率为9.29 g C m-2a-1。NPP显著增加的区域主要集中在玉林市的南部地区、浦北县的南部以及合浦县的西部地区;不显著增加的区域则主要分布于上游的河谷地区、浦北县与博白县的交错地带,另外灵山县的东部地区也有大面积分布;不显著减少的地区则主要分布于博白县的六万大山地区,该地区虽然植被的基质较好,植被净初级生产力的平均值较高,但是该区域的NPP在过去16年期间则呈现出微度减少趋势;显著减少的区域则主要集中在玉林市的市辖区周边以及合浦县南流江流域的出海口附近,其中,玉林市植被净初级生产力显著减少的原因主要与该地区的城市化水平有关,而南流江出海口NPP之所以呈现出显著减少趋势究其原因则主要与该地区植被受到大范围破坏,其植被大面积转变为虾塘和鱼塘有关。

图4 研究区净初级生产力空间变化趋势及显著性检验Fig.4 Spatial distribution of NPP sen trend and its Mann-Kendall test in study area

表1 研究区净初级生产力变化趋势

Sen trend：Theil-Sen趋势;Z：显著性水平

3.2.3植被净初级生产力Hurst指数及未来变化趋势

由图5可知,北部湾南流江流域净初级生产力Hurst的范围为0—0.99993,平均值为0.70,植被净初级生产力的反持续序列仅占到区域整体的NPP百分比为0.3986%,持续性序列占到区域整体的的NPP百分比为99.6014%。从图6中NPP的正态分布图中可知反持续性、弱反持续性、弱持续性和强持续性的面积百分比分别为0.34%、2.85%、75.73%和21.08%。北部湾南流江流域净初级生产力Hurst指数的正态分布图的空间形状表现出单峰右偏趋势,说明多年NPP的均值持续性趋势显著大于反持续性,由此可知南流量流域植被净初级生产力NPP整体的变化趋势则处于持续性的状态。从净初级生产力Hurst指数空间分布状况(图5)可以看出,南流江流域Hurst指数分布的趋势为：弱持续性呈现出大面积分布状态,强持续性分布于流域的河谷地带,而强反持续与弱反持续性则零星散布于研究区的上游以及下游地区。

为了更进一步明晰北部湾南流江流域净初级生产力NPP的变化趋势及其未来的可持续性特征,本研究将NPP Sen趋势的显著性检验值(图4)与Hurst指数的计算结果(图5)在ArcGIS 10.2中进行叠加分析,可以得到北部湾南流江流域植被净初级生产力增加、减少与未来可持续性的耦合图(图6)。由(图6和表2)可以看出：在NPP空间尺度上,南流江流域未来植被净初级生产力持续性增加的序列所占的百分比(51.64%)高于持续性减少的百分比(45.18%)。其中强持续性减少的组合所占的面积比重为5.52%,主要分布于玉林市的城乡过渡地区、博白县的东南部、灵山县的南部,其次,合浦县的南流江流域出海口也有大面积分布;弱持续性减少的组合所占的百分比最大,其值为39.66%,主要分布于南流江流域中游的六万大山地带,另外,流域下游的五皇山地区也有大范围分布;弱持续性增加的组合所占的百分比与弱持续性减少的比重相当,其值为36.07%,主要分布于流域上游河谷地区向山体的过渡地带,其次流域下游的河谷地区也有少量分布;强持续性增加的组合所占的百分比最小,仅为15.57%,主要分布于流域的上游的河谷地带,而强反持续与弱反持续性则零星散布于研究区的上游以及下游地区十万大山的南麓地带、浦北县的东南部,另外,钦州-合浦-灵山的交错地带也呈现出片状分布状态,净初级生产力Hurst在该类区域呈现出持续性增加的区域,预示着未来净初级生产力将有持续性增加趋势。

图5 Hurst指数和NPP持续性特征Fig.5 Hurst and sustainability of inter-annual NPP in study area

图6 Hurst+sen指数空间分布及指数正态分布图Fig.6 The spatial structure of Hurst+sen index and Normal distribution plot of the Hurst index

3.3 植被净初级生产力驱动因素

3.3.1气象因子的特征及其变化趋势

由于南流江流域可用的气象站点数据非常稀少,本研究将运用泰森多边形对附近的雨量站进行插值,可以得出南流江流域降水和气温的变化特征(图7)及其空间分布状况(图8)。南流江流域2000—2015年均降水量介于1361.50—2431.16 mm之间,多年均值为1802.09 mm,其中2008年的年均降水量最大,而最小值则出现在2004年,2000—2015年降水量呈现出显著增加趋势(P<0.01),增加的速率为177.35 mm/10 a;南流江流域2000—2015年均气温介于21.69—23.17 ℃之间,多年均值为22.53 ℃,其中2015年的年均气温最高,而最小值则出现在2011年,2000—2015年气温呈现出波动不显著下降趋势,下降的速率为0.019 ℃/10 a。

表2 净初级生产力NPP趋势以及Hurst指数统计特征

根据南流江流域2000—2015年的降水量和气温的多年均值空间分布特征(图8)可以看出,气温多年平均的空间分布情况则呈现出明显的地域分异规律,气温较高的地区主要分布在南流江流域的南部地区以及流域的河谷地带,而流域西北部地区的浦北县由于地处于六万大山和五皇山的交接地区,平均海拔较高(六万大山主峰六万顶海拔1115 m),导致区域气温偏低。南流江流域由于受到季风气候的影响,降雨量从流域南部到流域中游地区呈现出减小趋势,且分布不均匀,东北部地区以及流域下游南部地区之间地区的降雨较多,尤其是下游的合浦县,降雨量最多,流域的东南部降水量偏少。

图7 研究区年降水量和气温变化趋势Fig.7 Variation trend of annual precipitation and temperature in study area

图8 研究区多年平均降水量和气温空间分布Fig.8 Spatial distribution of annual precipitation and temperature in study area

3.3.2植被NPP与气象因子之间的关系

降水和气温等气象因子的变化对陆表植被的生长具有重要的作用,为了定量识别流域净初级生产力对气象因子的响应特性,本研究逐像元计算了流域2000—2015年的植被净初级生产力与年均降水量和年均气温之间的偏相关系数(图9),由图可知,净初级生产力对气温和降水的响应特性表现出显著的空间差异特征,且植被NPP对气温的相关系数高于对降水的偏相关系数。由植被NPP与年均气温的偏相关空间特征(图9)可知,植被NPP对气温的偏相关系数介于-0.87—0.95之间,平均值为0.09,由均值可知植被NPP与气温之间呈现出正相关关系,表明在南流江流域气温越高的地区净初级生产力越高。具体来说,在北部湾南流江流域,植被NPP与气温的偏相关系数正负相关区域分别占流域总面积的15.91%和80.09%,其中只有1.95%的正相关区域通过了显著性水平为0.01的检验,主要分布于玉林市的西北部地区,而21.93%的负相关区域通过了显著性水平为0.01的检验,主要分布于南流江流域的中游地区,其次浦北县的西南部也有斑块状分布。植被NPP对降水的偏相关系数介于-0.96—0.86之间,平均值为-0.25,由均值可知植被NPP与气温之间呈现出负相关关系,表明在南流江流域降水越高的地区净初级生产力越低。具体来说,在北部湾南流江流域,植被NPP与降水的偏相关系数正负相关区域分别占流域总面积的62.56%和37.44%,其中18.26%的正相关区域通过了显著性水平为0.01的检验,主要分布于浦北县与博白县的交接地带,而21.93%的负相关区域通过了显著性水平为0.01的检验,主要分布于浦北县的五皇山地区。综上,植被净初级生产力NPP与气温呈正相关性,与降水呈负相关性的特征。

图9 研究区NPP与气温和降水的偏相关系数Fig.9 The partial correlation coefficient between NPP and temperature and precipitation in study area

3.3.3植被NPP与土地利用之间的关系

流域净初级生产力NPP的变化主要受到气候变化和人类活动的影响,其中气候变化中的降水和气温对植被净初级生产力的变化具有重要的作用,而人类活动对净初级生产力的影响则主要通过土地利用方式的改变来体现。本研究为了揭示土地利用变化对植被NPP的损耗与增益效应,将南流江流域2000—2015年的土地利用图进行空间叠置分析,在北部湾南流江流域土地利用变化图谱中(图10),共有36类土地利用图谱单元,1类未发生转化,其中南流江流域上游玉林市的北部地区以及下游的合浦县是土地利用转移过程最为复杂的区域,在上游中耕地转变为建设用地的区域所占的比重最大,而下游则以耕地变为湿地地类为主,当地村民以渔业养殖为主,土地多用于虾塘和鱼塘养殖,分布区域主要集中在南流江流域的下游地区。将土地利用发生变化的图谱单元按照面积递降的顺序进行排序,在GIS中计算每个图谱单元的变化率(表3),其中前15类的土地利用图谱发生转移的面积累积百分比达到总变化面积的95.31%。

图10 研究区土地利用变化图谱及其空间集聚特征Fig.10 Land use changes during the period of 2000—2015 in study area

从表3和图10中可以看出,近16年来南流江流域土地利用类型面积变化的最大的图谱单元为“耕地-建设用地”,研究期内共有4780.62 hm2的耕地转化为建设用地,占总变化面积的百分比为34.24%,由土地利用面积变化导致的NPP损失值达到4715.62 t/a,其面积分布区域主要集中在流域上游玉林市的城乡过渡地带。林地转换为建设用地的面积所占的比重次之,占到变化面积的20.20%,由林地转化为建设用地所导致的NPP损失值达到2835.19 t/a,变化图谱主要分布灵山县与合浦县的交接地带,该地区主要是受到灵山县的扩张,周边林地受到人类活动的滥砍乱伐导致林地面积大幅度降低;耕地转换为林地的图谱单元位居第三位,其面积为1408.59 hm2,占所有变化面积的10.09%,该图谱单元所导致的NPP损失值达到2510.03 t/a,空间分布区域上主要集中在南流江流域中游的博白县南部的河谷地区,说明退耕还林还草使得人类对于森林干扰减少,使其流域中游地区植被得到较好恢复;而草地转换为建设用地的比重最小,其值仅为1.10%,导致NPP减少的损失值仅为184.63 t/a。

3.3.4不同分区气候因子及土地利用变化与NPP的关系

从表4中可以看出,在不同波动区、不同变化区以及不同持续区多年平均NPP与降水的偏相关系数都为负值,而与气温的相关系数均为正值,说明在北部湾地区气候因子中的气温对植被的净初级生产力起到控制作用。就不同波动区而言,高波动区植被对气温的相关系数高于降水,且高波动性区域土地利用变化的幅度显著大于低波动性区域;从不同变化趋势区可以看出,NPP增加的区域与降水和气温的相关系数的绝对值都为0.27,但是不同土地利用变化的面积核密度值在NPP增加区中小于NPP减少区,说明在南流江流域NPP的增加趋势与土地利用的稳定性有着直接的关系,土地利用变化越剧烈,区域的NPP增长的趋势则越慢;就NPP未来变化趋势而言,持续性序列与降水和气温的相关系数都小于反持续性序列,且土地利用面积变化核密度值在持续性区域中小于反持续性区域,这与不同变化趋势区NPP与土地利用的关系类似,流域土地利用变化越剧烈则区域的NPP的可持续性则呈现出降低趋势。

表3 研究区2000—2015年土地利用变化类型变化排序及其导致的NPP减少值

表4 不同分区气候因子及土地利用变化与NPP的关系

4 讨论

4.1 植被净初级生产力模型精度验证

图11 CASA模型估算NPP与MODIS17数据的关系Fig.11 Correlation of simulated NPP based on CASA model and MODIS17 data

区域净初级生产力的精度评价一直以来都是遥感学界以及生态学中的难点和有争议的环节,NPP估算的结果验证一般来说分为两种,一种为实测值验证,这种方法是通过估算值与实测的NPP数值进行对比验证。另外一种是相对法,即通过将模型的估算结果与其他模型结果或者其他遥感产品进行对比来评价。第一种方法是通过野外获得研究区的实测生物量,其数据比较可靠,但是目前南流江流域植被生物量的实地调查数据目前还未开展相关的研究工作。基于此,本研究采用第二种方法对CASA模型的估算结果进行精度验证,所采用的精度验证的产品为MODIS 17 NPP产品,该产品是基于BIOME-BGC模型计算出的全球净初级生产力数据[27],国内也有大量的研究应用该产品验证了CASA模型估算结果的可靠性,如谢宝妮等[28],孙庆玲等[29]和张继平等[30]分别使用该产品对黄土高原、武陵山区以及三江源地区的净初级生产力进行了精度和可靠性验证。本研究在进行模型验证时,由于原始的MODIS17 NPP产品的空间分辨率为1000 m,本研究首先对该数据进行重采样到250 m,之后使用GIS的随机采样工具分别在MODIS17数据和本文结果上生成30%的采样点,以此来分析二者之间的相关关系。从图11可以看出,本文估算结果与MODIS17数据存在着显著的相关性(P<0.01),相关系数为0.7943,但是本研究通过CASA模型估算的结果略低于MODIS17的平均值,这可能与MODIS17采用的是16天合成的植被指数中可能存在阴天或多云等情况,使得光合有效辐射吸收系数比实际值偏大,导致光合有效辐射的结果偏大,进而增大NPP的估算值。

4.2 植被净初级生产力贡献率识别及其未来研究方向

由于本研究中的南流江流域气象站点偏少,尤其是太阳辐射站点较少,本文虽然采用了插值方法对气象数据进行了处理,但是由于该流域广泛存在高山地貌、丘陵地貌以及河谷地貌。因此还需要发展更为精确的插值方法以提高模型估算的精度。另外,本研究虽然采用偏相关分析方法量化了气候因子对植被NPP的响应关系,利用土地利用变化转移矩阵分析了流域土地利用变化导致的NPP增加和减少值,但是未对气候变化和人类活动对流域的贡献率和控制区域进行定量识别,因此未来可以从两方面进行入手,一方面可以采用相关的数据模型,如残差分析模型[31]建立回归分析方程,对每年的NPP残差序列进行Sen趋势度分析,以得到的趋势值为标准,如果趋势值为正值,说明NPP象元值受到人类活动的影响,反之,则NPP象元值则受到气候变化的影响,利用正负值所占的比例可以计算气候变化和人类活动对植被NPP的贡献率。另一方面,可以通过统计研究区牲口的数量、草地、坡耕地的面积在研究区植被恢复方面的贡献比例。因此,关于气候变化和人类活动对净初级生产力的贡献率将是本研究以后重点的研究方向。

5 结论

本文以北部湾南流江流域为研究对象,基于光能利用率模型(CASA),利用遥感、气象和植被等数据估算了研究区2000—2015年流域的净初级生产力,借助于Theil-Sen趋势、Mann-Kendall检验以及Hurst指数等数理统计方法对研究区NPP的时空变化特征、未来趋势及其驱动因素进行了定量化分析。研究结果表明：

(1)时间尺度上,2000—2015年北部湾南流江流域植被净初级生产力NPP总体上呈现出波动上升趋势,增速为44.03 g C m-2(10 a)-1,明显快于广西自治区2000—2011年植被净初级生产力,未通过显著性水平为0.01的检验。流域上游和下游地区植被NPP快于全区,而中游地区慢于全区,其中,流域上游、中游和下游的净初级生产力年平均增加速度分别为：54.17g C m-2(10 a)-1、31.74g C m-2(10 a)-1和46.1g C m-2(10 a)-1。多年平均净初级生产力介于911.30—1075.44 g C/m2之间,平均值为995.25 g C/m2,最大值出现在2011年,其值为1075.44 g C/m2,高出平均值80.19 g C/m2,最小值则出现在2006年,低于平均值83.95g C/m2。(2)空间尺度上,流域植被净初级生产力的分布呈现出明显的地域分异规律,中游植被净初级生产力最高(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之间,流域净初级生产力在空间上呈现出增加的趋势,其中NPP显著减少的区域所占1.58%,轻微减少占46.57%,轻微增加占43.25%,而明显增加占8.60%。

(3)空间波动性上,流域净初级生产力的变异系数较大,其值介于0.01—0.70。流域NPP的多年平均值呈现出明显的空间分异规律,主要以低波动性为主,其中,洪潮江水库、小江水库周边以及玉林市的城乡建设用地扩张区域处于高波动状态,而流域的中部六万大山以及五皇山地带则处于低波动状态。

(4)未来变化趋势上,流域净初级生产力Hurst的范围为0—0.99,平均值为0.70,植被净初级生产力的反持续序列仅占到区域整体的NPP百分比为0.3986%,持续性序列占到区域整体的的NPP百分比为99.6014%,Hurst指数正态分布图呈现单峰右偏分布,即NPP均值的持续性序列显著大于反持续性序列,预示着流域NPP未来处于持续增加的趋势。

(5)驱动机制上,流域NPP与多年平均气温呈正相关关系,与年均降水量呈负相关关系,表明温度是影响该流域植被NPP的主要气候因子。由耕地转化为建设用地,导致的NPP损失值最大,其值达到4715.62 t/a,林地转换为建设用地所导致的NPP损失值次之,其值为2835.19 t/a,而草地转换为建设用地导致NPP损失值最小,其值仅为184.63 t/a。