香格里拉高山松地上生物量及遥感因子的时空变化特征研究

王书贤,张加龙,鲍 瑞,韩东阳,廖 易

(1.西南林业大学 林学院,云南 昆明 650224;2.云南省遥感中心,云南 昆明 650034;3.国家林业和草原局 西南调查规划院,云南 昆明 650031)

森林地上生物量(abovegroundbiomass,AGB)是估算森林生态系统碳储量、碳通量的重要参数[1-2],森林AGB空间分布变化能够更好地量化碳储量的空间分异特征,对于研究森林生态系统碳循环具有重要意义。现阶段主要运用ArcGIS软件中的空间插值法来对森林AGB的空间异质性进行分析研究,而用遥感影像数据反演森林AGB,然后对其进行空间异质性分析的研究较少[3-6]。基于遥感与地统计方法,能够有效揭示属性变量在空间上的分布特征及其变异性[7-13],当前运用空间地统计学方法的林业研究主要集中于森林土壤空间变异、物种空间分布、森林干扰、林分因子和种子资源这五个方面[7,14],而森林AGB空间的统计分析研究还较少。在森林AGB的估算中,AGB的空间关系未引起研究者的重视。当前估算森林AGB的方法还是以统计分析为主,这种方法不但不能准确、及时地得到森林AGB的空间分布信息,而且也不能有效反映森林AGB的空间自相关性和异质性[13,15]。

香格里拉地处青藏高原地区东南部,其气候变化是全球气候变化的敏感响应区[16]。高山松(Pinusdensata)是该地区的优势树种之一,探究其AGB的时空变化特征对研究高原地区生态系统的碳循环,提高对该类型地区森林生态价值的认识具有重要意义[15]。本研究以香格里拉市高山松区域7个时期(1987-2017年)的Landsat遥感影像为基础,对不同阶段高山松AGB采用地理信息系统(GIS)和地统计学2种方法进行时空变异分析,探究该地区高山松AGB的空间变化特征,为日后监测AGB空间变化和模拟预测提供理论依据。

1 研究区概况

研究区位于云南省迪庆藏族自治州香格里拉市(99°20′-100°19′E,26°52′-28°52′N),地处云南省西北部,迪庆州东北部。境内海拔最高5 545 m,最低1 503 m,相对高差4 042 m。香格里拉市位于高海拔低纬度地带,属山地寒温带季风气候。境内地貌复杂多样,地势北高南低,立体气候明显,植被南北分布差异明显。研究区植被类型主要为寒温性针叶林,其树种有高山松、丽江云杉(Picealikiangensis)等,其中高山松森林覆盖面积占优势树种面积的22.7%[15,17]。

2 材料与方法

2.1 数据来源

2.1.1 样地数据 样地数据来自国家森林资源连续清查项目,时间为1987-2017年,每隔5 a进行一次样地实测,共7期。每块调查样地面积为0.08 hm2(28.28 m×28.28 m),共有136块高山松样地,各年份样地数量各异,部分样地为复测样地。

2.1.2 遥感数据 采用美国地质勘探局(United states geological survey,USGS)和美国国家航空航天局(national aeronautics and space administration,NASA)已经过校准的Landsat地表反射率数据(SR)作为时间序列遥感数据,运用可提高数据精度的最优斜率算法(the best slope temporal segmentation filtering,BSTS)对数据进行了时间序列滤波重构[17]。

2.1.3 其他数据 DEM(数字高程模型)数据来自地理空间数据云(http://www.gscloud.cn/),分辨率为30 m×30 m。

2.2 基于随机森林回归的生物量反演结果

2.2.1 遥感特征变量的选择 遥感特征变量的选取参考了香格里拉高山松AGB估测研究中常用的遥感因子,纹理因子是基于ENVI 5.3软件对Landsat的6个波段分别计算了7种纹理因子,总共计算了3×3、5×5、7×7、9×9、11×11、13×13、15×15、17×17、19×19共 9 种不同的窗口大小。由于纹理特征种类过多,单独对 378种纹理特征进行了筛选。将全部的 378个特征因子数据输入随机森林模型中,通过迭代筛选,选出15个遥感与纹理特征(表1)。

表1 15个遥感因子特征变量Table 1 Characteristic variables of 15 remote sensing factors

2.2.2 多期高山松AGB估测数据 使用筛选出的15个遥感与纹理特征变量,基于随机森林估测模型对1987、1992、1997、2002、2007、2012、2017年香格里拉高山松地上生物量进行估测反演,得到各年份的高山松地上生物量反演数据结果。

2.3 地统计分析方法

在ArcGIS 10.7中对7期的高山松AGB随机取点,样本量取值为200个随机点,进行标准化异常值处理后,各年份样本量见表2。应用GS+9.0地理空间分析软件对高山松AGB的随机点计算Moran'sI系数和半变异函数,然后对其进行空间自相关性和空间异质性分析。

表2 随机点在各年份中的分配Table 2 Distribution of randomplots in each year

2.3.1 空间自相关性分析方法 通过Moran'sI系数计算高山松AGB和遥感因子的空间自相关度,并用Moran'sI系数图进行分析,计算公式[5,9,11,18-20]如下。

(1)

Moran'sI系数的取值范围在( 0,1],说明高山松AGB和遥感因子的空间分布呈正相关是聚集化状态,有整体性分布; Moran'sI系数的取值范围在[-1,0),说明高山松AGB和遥感因子的空间分布呈负相关,是破碎化状态,没有整体性分布。

2.3.2 空间异质性分析方法 运用半变异方差函数对高山松AGB和遥感因子进行空间异质性分析,计算公式[5,9,11,18-20]如下。

(2)

式中:N(h)表示某一空间方位上相距距离为h的固定样地对数;Z(xi)和Z(xi+h)分别表示固定样地xi和与xi相隔距离h处的生物量值。

采用4个参数进行空间异质性分析,分别为块金值C0、基台值(C0+C)、变程A0、结构比C0/(C0+C)。块金值C0为区域变量的方差;基台值(C0+C)反映区域化变量在研究范围内变异的强度[21];变程A0表示变量空间自相关的有效范围,距离越近,空间特征的相似程度越高,超出该值范围则变量没有自相关关系;结构比C0/(C0+C)判断随机因子对空间总变异的影响程度[9],也可表示变量的空间相关性程度,若比例<25%,说明变量具有强烈的空间自相关;若比例在25%～75%,变量具有中等空间自相关;比例>75%时,变量空间自相关很弱[22]。

3 结果与分析

3.1 高山松AGB的时间变化

根据AGB反演结果(表3),1987-1997年的高山松分布面积减小,由17.156万hm2减少到17.059万hm2,期间的AGB也在快速减少,由951.911万t减少到810.696万t;1997年后香格里拉市高山松的分布面积开始增加,截至2017年,高山松分布面积达18.484万hm2,但该区域的高山松AGB并未随之相应增加,而是呈减少趋势,2017年的高山松AGB总量为741.084万 t。

表3 1987-2017年高山松AGB反演结果Table 3 Inversion results of P.densata AGB during 1987-2017

3.2 高山松AGB的空间自相关性

从图1可以看出,1987年高山松AGB在55 000 m范围内为聚集状态,呈连片化、整体性空间分布,在南北方向的聚集状态为最好,空间格局尺度达85 000 m左右;1992年高山松AGB在60 000 m内为聚集状态,在南北方向的聚集状态为最好,空间格局尺度为75 000 m,但总体的AGB空间分布较为破碎;1997年高山松AGB空间自相关性较差,整体呈破碎化空间分布,不具有整体性;2002年和2007年的高山松AGB在20 000～50 000 m都呈现破碎化空间分布,在南北方向和东西方向也在该空间距离内呈现破碎化的空间分布;2012年在20 000 m内和50 000～70 000 m具有聚集状态的高山松AGB,其南北方向、东北-西南方向和西北-东南方向也在该空间距离内具有同样空间分布状态的高山松AGB;2017年的高山松AGB在不同空间距离上具有不同的空间分布,全局和各方向的空间分布趋势大体一致。

图1 1987-2017年高山松AGB的Moran's I系数Fig.1 Moran's I coefficient of P.densata AGB from 1987 to 2017

3.3 高山松AGB的空间异质性

从表4可以看出:经过GS+9.0的筛选,1987-2012年的AGB的最优半变异函数模型均为指数模型,其结构比C0/(C0+C)都接近0。除1992年和1997年的变程A0>10 km外,其余年份的变程A0都在10 000 m内。经过变程A0筛选出与各年份AGB具有空间分异特征相似的遥感因子,分别是NBR、B3/Albedo、FR11B2VA和FR17B7DS,其中,NBR和B3/Albedo为出现频率最高的遥感因子;由于2007年AGB与遥感因子的变程A0相差较大,因此该年没有与其具有相似空间分异特征的遥感因子;2017年高山松AGB的半变异函数模型的最优模型为球形模型,其结构比C0/(C0+C)是近30 a来最接近0的,说明2017年香格里拉高山松AGB的空间自相关程度非常高。根据变程A0的相近程度,具有与2017年AGB相似空间分异特征的遥感因子为B3/Albedo。

表4 1987-2017年半变异参数Table 4 Semi-variation parameter table in 1987-2017

4 结论与讨论

1987-1992年,南北方向的高山松AGB在60 000 m内的空间自相关性最好,呈聚集状态,高山松为连片化分布;1997-2007年的高山松AGB全局空间自相关性差,呈破碎化空间状态;而2017年的高山松AGB在不同空间距离上呈现出不同的空间分布状态。

森林AGB的空间异质性产生的主要原因可能是立地条件的空间异质性[9-10]。1987年和2002-2017年的结构比C0/(C0+C)都接近0,变程A0都在10 000 m内,说明该时期的高山松AGB和遥感因子变量的异质性变化是由空间自相关引起的,即其异质性变化是由结构性因素(立地条件,海拔、坡度等)及水热条件等自然因素引起的。该时期的AGB空间变异受人为干扰因素(砍伐、污染等)的影响较小,也说明1987年和2002-2017年香格里拉高山松地区的森林保护良好。

1992年和1997年的C0/(C0+C)虽然接近于0,但变程A0>10 000 m,说明这2个时期的高山松AGB和遥感因子变量的异质性变化可能是由人为干扰等外界因素引起的。结合空间自相关可判断该时期的AGB不具有空间整体性。

追溯香格里拉市各年高山松变化情况可知:20世纪80-90年代,当地为了发展经济盲目新建了许多森林木材工厂企业,从而加剧了森林资源的消耗,因此可解释引起该变化的因素主要是人为干扰[16]。1997年后国家出台全面禁止天然林采伐的政策,同时实施了退耕还林等一系列森林生态保护工程和森林保护措施,使得香格里拉高山松开始恢复,从2002-2017年的研究结果也反映了该恢复情况[15,22,23]。研究结果表明,基于遥感与地统计学相结合的方法能够高效地获得AGB及其相关的遥感因子的时空变异特征,即可在一定变程A0范围内判断高山松AGB的分布,以及可大体判断受何种类型因素的影响,该方法可为日后高山松AGB的时空模拟预测和反演制图奠定一定的理论依据。

相较于空间插值法受样本点的影响[22],本研究采用经过BSTS滤波算法重构时间序列数据来估测反演高山松AGB提供了更可靠的数据源[3,4,11,25-26];受地统计方法样本量的取值大小,以及获取样地点实测数据局限性的影响,样地点的样本量采用随机取点的方式满足研究条件。样本量的取值多少在此类研究中暂无具体的科学依据,后续对该部分内容值得探究;本研究运用了30 m分辨率的遥感影像数据,若用分辨率更高的遥感影像数据,可能会获得更加精确的结果。