青海三江源区高寒植被地表反照率变化及其辐射温度效应

朱躲萍,叶 辉,王军邦,赵烜岚,左 婵,芦光新,张法伟,李英年

1 九江学院旅游与地理学院,九江 332005 2 中国科学院地理科学与资源研究所 生态系统网络观测与模拟重点实验室 生态系统大数据与模拟中心,北京 100101 3 长江大学园艺园林学院,荆州 434000 4 青海大学农牧学院,西宁 810016 5 中国科学院西北高原生物研究所,西宁 810008

土地利用变化和温室气体排放被认为是全球气候变暖的主要因素,并分别从其生物地球物理过程及生物地球化学过程开展研究[1—3]。表征下垫面反射短波辐射通量与入射短波辐射通量比值的地表反照率,是影响地球能量平衡的关键变量[4—5],是生物地球物理过程中植被变化对气候的反馈机制之一[1, 6—8]。

在植被变化对气候的反馈机制研究中,对地表反照率的影响存在不同的解释[3]。模型模拟研究表明,若全球尺度毁林,其净气候效应是通过改变地表反照率与蒸散而抑制由于毁林产生的温室气体排放导致的气候变暖[9],这也得到了观测数据的支持[10]。过去30年全球植被绿度增加减缓了全球陆地气候变暖增温的12%,其中地表反照率降低贡献了6%[2]。而在青藏高原的研究发现,该地区植被绿度增加,尽管地表反照率降低使地表吸收的辐射增多,但同时也使地表蒸散增强,最终导致了减缓气候变暖的降温效应[11]。然而草地退化与恢复作为土地利用变化的一种渐变方式,所导致的地表反照率的变化及其影响机制,特别是其辐射效应,仍然缺乏研究[2, 12—15]。

作为长江源、黄河源和澜沧江源发源地的青海三江源地区,其分布最广的草地生态系统不仅维持着多种生态功能,也为当地牧民畜牧业发挥重要的生活和生产功能；其敏感而又脆弱的生态系统,是我国重要的生态安全屏障之一[16],在全球气候变化研究领域中占据着重要地位。本文基于2001—2018年中分辨率成像光谱仪(MODIS)地表反照率等数据[17],研究三江源区土地覆盖与利用变化对地表反照率的时空变化的影响,及其地表辐射温度效应,以深入理解区域土地利用与覆被变化的生物地球物理过程与机制,为适应全球气候变化、生态系统管理及科学决策提供支撑。

1 研究区概况

三江源地处青藏高原腹地(图1),地理范围约为31°39′—36°12′N,89°45′—102°23′E,海拔为3335—6564 m,总面积为36.3万km2[16],是长江、黄河和澜沧江的发源地,素有“中华水塔”之美称。三江源具有冷热两季交替、干湿两季分明、年温差小、日温差大、日照时间长、辐射强烈、无明显四季区分的气候特征,为典型的高原大陆性气候[18]。

图1 三江源区地理位置及主要土地覆盖类型

根据2018年30 m空间分辨率土地利用数据一级分类(图1)[19],三江源的植被以草地为主,占全区面积的71.2%,其中高、中、低覆盖草地分别占全区的5.7%、24.9%、40.6%；未利用地占16.7%,森林仅占全区面积的4.3%(表1)。长江源园区、黄河源园区以及澜沧江源园区的草地面积占比分别为85.3%、89.4%和50.8%,未利用地面积分别占40.4%、3.4%和3.6%,森林为0.01%、0.30%和0.90%。其中,长江源园区、黄河源园区以及澜沧江源园区草地高、中、低覆盖草地面积占比分别为1.1%、42.6%、41.6%,0.3%、23.6%、65.5%和2.5%、12.7%、35.6%。

表1 2018年三江源以及三个园区土地利用与土地覆被类型及面积/(×104 km2)

2 研究数据与方法

2.1 研究数据

2.1.1地表反照率

采用2001—2018年每16天0.05°(约5 km)空间分辨率MODIS的地表反照率产品(MCD43C3),是基于分别搭载在Terra和Aqua两个卫星上的传感器数据合成的数据产品[20]。该产品包括太阳时中午的黑空短波反照率、白空反照率、可见光以及近红外波段反照率,精度评估基于3个Ameriflux通量塔及53个FLUXNET塔点观测的反照率数据,均方根误差小于0.013[21]。本研究选择了黑空短波反照率数据进行了分析研究,数据下载自美国地质调查局官网。

2.1.2归一化植被指数

归一化植被指数(Normalized Difference Vegetation Index,NDVI),是利用2001—2018年250 m空间分辨率的MODIS波段反射率数据产品(MOD09Q1)计算而得。在计算得到NDVI后,通过Timesat软件进行了SG滤波[22—23],以消除云等对NDVI数据的影响。

2.1.3气象数据

以中国气象科学数据共享网(http://date.cma.cn/)提供的研究区及周边气象台站观测的降水量、最高和最低气温日值数据,采用气象数据空间插值软件ANUSPLIN,经空间插值得到三江源8天时间分辨率250 m空间分辨率的栅格数据,用于分析。ANUSPLIN软件是基于普通薄膜和局部薄膜样条函数多变量数据进行内插的工具,能同时进行多表面空间插值[24],插值时以90 m空间分辨率SRTM数字高程模型(DEM),经空间重采样为250 m的DEM数据作为协变量。

2.1.4土地利用和覆盖(LUCC)数据

本研究所用2018年三江源地区LUCC数据,来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。该数据空间分辨率为30 m,是根据土地资源的利用属性和自然属性,以及地貌位置等信息通过人工目视解译的方式进行分类；其二级类型包括水田、旱田、有林地、疏林地、灌丛、高覆盖草地、中覆盖草地、低覆盖草地、水域、建设用地、湿地、荒漠及未利用地等[19]。

2.1.5地表反照率地面观测数据

以通量观测塔测定的上行和下行短波辐射,计算得到的反照率作为地表反照率地面观测数据,对地表反照率遥感数据进行了精度评价。通量塔测定数据来自中国通量观测网络(ChinaFLUX)海北高寒草地长期定位观测站(37.6°N,101.3°E)。通量塔所提供的数据为2015—2017年每30分钟的数值,在精度评价时,以每天11：00—15：00时范围内数值计算均值,作为每日地表反照率,然后计算月平均值,与对应时段MCD43C3月平均数据进行对比分析。

2.2 研究方法

2.2.1地表反照率数据质量控制

对于16天和0.05°时空分辨率MCD43C3地表反照率数据产品,经投影转换、范围裁剪、数据空间重采样以及缺值插补,处理为8天250 m时空分辨率数据,以与植被指数及气象数据时空匹配。针对本研究所用的黑空短波辐射反照率数据,根据其质量标识数据(QA)进行了如下质量控制：选择QA小于等于3的值作为有效值,对QA大于等于4的像元值,作为缺值进行插补。当QA等于4,地表反照率值小于三倍有效值的标准差时,保留原值,否则以QA小于3的像元有效值的平均值进行插补[25]。

2.2.2太阳短波辐射

到达地面的太阳总辐射的计算,一般需要考虑地形、云量、大气透明度、太阳光线入射角以及太阳高度角等因素[26],这里根据文献采用下式计算：

(1)

式中,Ra为地面晴空太阳总辐射(MJ m-2d-1),N为可照时数,n为实照时数,n/N为日照百分率,as和bs的值由地区确定[27]。

2.2.3归一化差植被指数(NDVI)

归一化差植被指数(NDVI),MODIS波段反射率数据计算得到[25],计算公式如下：

(2)

式中,NR为近红外波段,RR为可见光波段,即近红外波段和可见光波段之差除以近红外与可见光波段和。

2.2.4随机森林回归算法

随机森林算法最早是在2001年由Breiman提出,主要用于解决分类和回归的问题,是通过集成学习思想将多棵决策树集成的一种机器学习算法[28]。随机森林算法能提高预测能力,具有抗拟合能力,模型的泛化能力强,参数优化简单[29—30]。本研究基于地表反照率逐像元的年际变化影响因子分析,选择NDVI、最高气温、最低气温以及降水4个变量,构建了地表反照率年际变化的随机森林回归模型,对变量重要性进行评估,来判断像元尺度对地表反照率年际变化影响因子。本文通过Python中RandomForestRegressor库实现了随机森林回归模型,变量重要性X按式(3)计算：

X=∑(E1-E2)/Ntree

(3)

式中,E1为每一棵树对应的袋外数据误差,E2为对袋外数据的变量加入噪声干扰,重新排序和再次计算的袋外数据误差,Ntree随机森林测试样本。

2.2.5一元线性回归分析

本研究采用最小二乘法拟合2001—2018年地表反照率的变化趋势,并进行像元尺度和分土地覆盖类型的年际变化分析,计算给出统计显著性水平P值。

2.2.6地表反照率变化的辐射温度效应计算

在陆面能量收支对气候变化的贡献计算中,不考虑各种反馈过程,仅考虑地表反照率对辐射平衡的影响,辐射平衡状态时[31],计算公式如式(4)：

(4)

式中,RS是太阳短波辐射(MJ m-2d-1),a是地表反照率,ε是比辐射率(灰体系数,取1.0),δ是斯蒂芬-玻耳兹曼常数(为5.67×10-8W m-2K-4)。根据式(4),草地覆盖度变化引起的辐射温度效应,采用下式计算[27]：

(5)

式中是草地覆盖度变化引起的辐射温度效应,Q为研究时段内平均年太阳短波辐射,a1、a2分别为前后两个研究时段不同覆盖度草地的平均地表反照率。本文中分别以2001—2004年期间和2005—2018年期间为草地覆盖度变化前后时段,分别组合计算高、中、低覆盖度草地变化的辐射温度效应。

2.2.7植被覆盖度

根据李苗苗等像元二分模型计算植被覆盖度[32]：

(6)

式中,NDVIS生长季平均归一化差植被指数NDVIveg和NDVIsoil分别为植被完全覆盖时和无植被覆盖时NDVIS值,以土地利用单元内频率为95%和5%的NDVIS值分别确定。

3 结果与分析

3.1 地表反照率空间分布

如图2所示,卫星遥感的地表反照率数据产品MCD43C3,在以海北站通量塔测定的下行和上行短波辐射计算的地表反射率进行精度评价时,二者之间存在显著的季节变化一致性,复相关系数R2为0.79,均方根误差为0.022,显著性水平P小于0.001,表明该数据产品在三江源区域具有较高的数据精度和可靠性。

图2 基于2015—2017年海北通量塔观测数据对MCD43C3产品的精度评价

如图3所示,2001—2018年三江源区平均反照率空间分布具有较大的空间变化,总体上呈自东南向西北逐渐增大的变化格局。区域平均值为0.163±0.027,主要集中在0.12—0.18之间,占总面积的76.4%；部分地表反照率大于0.18的地区,占总面积的21%,主要分布在三江源东北部和中西部；最高值达到0.56,分布于研究区西部地区戈壁、沙地以及裸土地等未利用地；研究区东南部地表反照率在0.12以下。长江源、黄河源和澜沧江源三个园区的地表反照率差异较大,以长江源园区最高(0.177±0.036),澜沧江源园区(0.156±0.002)略高于黄河源园区(0.153±0.037)。

在图3中,三江源区各土地覆盖类型的地表反照率,未利用地的最大,为(0.176±0.006),林地的最小(0.136±0.02)；不同覆盖草地间,低覆盖草地的地表反照率高于高覆盖草地,具体地高、中、低覆盖草地地表反照率分别为(0.154±0.003)、(0.155±0.003)、(0.163±0.004)。值得注意的是长江源园区的未利用地的地表反照率,为(0.182±0.008),相比黄河源园区(0.162±0.006)和澜沧江源园区(0.156±0.004)更高,因为未利用地主要集中在长江源园区,而黄河源园区和澜沧江源园区主要分布着低覆盖草地(表1),在澜沧江源园区的森林(0.9%)的地表反照率较长江源园区(0.01%)和黄河源园区(0.3%)更低,可能是因为林地的冠层会影响地表反照率,冠层可以吸收可见光波段,其面积的增加会使其对反照率的吸收增加,导致地表反照率下降。

图3 2001—2018三江源生长季年均地表反照率空间分布和三江源及三个园区各LUCC类型生长季平均地表反照率柱状图

3.2 地表反照率年际变化趋势

如图4所示,2001—2018年三江源区生长季地表反照率变化整体呈每年不显著下降趋势。全区平均下降速率为(0.152±0.763)%/10a(P=0.47),趋势速率的空间变化范围为-0.010—0.007之间。呈上升和下降趋势的地区分别占全区总面积的41.42%和58.58%,其中,显著下降区域占总面积的8.4%,主要分布在三江源区的东北部。呈显著增加的面积仅占全区1.9%,零散分布在三江源东部。长江源园区、黄河源园区以及澜沧江源园区平均地表反照率的年际变化为统计不显著趋势下降(P>0.05),每10年分别下降(0.078±0.900)%、(0.215±0.740)%、(0.152±0.450)%的速率下降。

如表2中,三江源不同土地覆盖与利用类型生长季地表反照率均以不同趋势下降。其中高、中、低覆盖草地分别以-0.0003%/10a、-0.0013%、0.0001%/10a速率下降,而林地和未利用地的下降速率分别为每10年0.0003%和0.0014%,但下降趋势均统计不显著(P>0.05)。

表2 三江源植被类型年际变化趋势、相关系数及其显著性

3.3 地表反照率年际变化素

应用随机森林以2001—2018年三江源区域气候因子(最高气温、最低气温、降水)和植被因子(NDVI)为自变量,以地表反照率为应变量,自变量的权重如图5所示,NDVI、降水、最高气温、最低气温为主导因子的区域,分别占全区面积的40.7%、16.4%、29.0%和13.8%,表明植被是影响地表反照率年际变化的主要因素,且各自变量的权重值分别为(0.31±0.149)、(0.22±0.106)、(0.21±0.098)、(0.26±0.131),如图6所示,表明植被因子对地表反照率的年际变化具有最高的贡献。NDVI的权重值,在三个园区间相比时,在黄河源园区贡献最高(0.319±0.148)、其次是澜沧江源园区(0.316±0.151),在长江源园区贡献最低(0.283±0.131)。

图5 2001—2018年三江源生长季地表反照率年际变化的贡献因子权重：NDVI、降水、最低气温和最高气温

图6 2001—2018年三江源区及长江源园区、黄河源园区、澜沧江源园区(NDVI、降水、最低气温和最高气温)各因子对地表反照率年际变化的贡献频数分布

3.4 辐射温度效应

以2001—2004年期间为基准,采用2001—2018年生长季平均太阳短波辐射,以消除太阳辐射变化的影响,定量分析2005—2018年期间生长季地表反照率变化的辐射温度效应,结果如图7,期间植被覆盖度的变化全区平均值为(0.01±0.08),即增加了1%[33],增加区域占全区面积的57.1%；地表反照率变化整体上呈辐射增温效应,增温(0.11±0.42)℃,辐射温度变化值主要集中在-0.3—0.6℃之间,占全区面积的84.8%；空间上表现为中西部以及东北部为辐射增温效应,东南和西北部为辐射降温效应,分别占全区面积54.86%和45.14%。

空间叠加分析表明,植被覆盖度增加区域既存在辐射增温效应,也存在辐射降温效应,分别占全区32.6%和24.5%。辐射增温效应主要分布在东北部和黄河源园区,植被覆盖度增加,地表反照率降低,地表吸收太阳净辐射增多,进而对气候产生了变暖的正反馈[34—35]；植被覆盖度增加伴随着的辐射降温效应主要分布在东南部,主要土地覆被类型为森林,说明这些地方森林绿度增加导致了地表反照率上升,从而导致净辐射吸收减少表现出辐射降温效应,已有研究也发现森林变绿会导致地表反照率增加[25]。

三江源42.9%的区域植被覆盖度下降,也是存在辐射增温和辐射降温效应,分别占全区面积的21.7%和21.2%,主要分布于三江源中西部和西部地区。植被覆盖度下降导致的地表反照率下降表现出的辐射增温效应,可能是生长季土壤冰层融化、土壤水分增多,而土壤水分增多会引起地表反照率下降[36—37],从而净辐射吸收增多表现为辐射增温效应,西北部植被覆盖度下降导致地表反照率上升,而地面吸收净辐射减少表现为辐射降温效应。在黄河源园区植被覆盖度的增加导致了辐射增温效应,值为(0.24±0.37)℃,而长江源园区和澜沧江源园区植被变化引起的辐射温度效应不显著。

图7 2005—2018年相比于2001—2004年基准期三江源生长季植被覆盖度变化及地表反照率变化的辐射温度效应

三江源生长季2015—2018年与2001—2004年高覆盖草地渐变到低覆盖草地NDVI值之差以及低覆盖草地过渡到高覆盖草地NDVI值之差与相应时段地表反照率变化之差引起的辐射温度效应关系均呈正相关,如图8所示,上升趋势分别为0.060(R=0.34,P<0.001)和0.058(R=0.30,P<0.001),其中低覆盖草地渐变到高覆盖草地引起的草地覆盖度增加(草地覆盖度变化大于0)导致的辐射增温效应(大于0℃),占选取随机点总数的57.3%,高覆盖草地过渡到低覆盖草地引起的草地覆盖度降低(草地覆盖度变化小于0)导致辐射降温效应(小于0℃),占选取随机点总数的56.3%,从而得出生长季草地覆盖度增加主要导致的地表反照率减少的辐射温增效应,而植被覆盖度变低主要导致的地表反照率增加的辐射降温效应,这与地表反照率上升和下降引起的地面吸收的净辐射下降和上升导致的辐射温度效应相一致。

图8 2005—2018年与2001—2004年三江源生长季高低覆盖草地对应的NDVI值之差与辐射温度效应散点回归图

4 讨论

土地覆被变化导致地表辐射特性参数改变进而影响地表辐射收支,其中作为关键物理参数的地表反照率的较小改变,将影响地气系统的能量收支,进而影响到区域气候[38]。三江源区域地表反照率在空间上呈现从东南到西北上升的分布特征,空间上差异性明显,这与赵之重等人得出的三江源地表反照率空间分布规律具有较好的一致性[39]；其中,地表反照率较低的东南部,分布着森林[40],而森林较其他覆盖类型具有较低的地表反照率[41—44]。三江源区生长季地表反照率变化整体呈现不显著下降趋势,且相对较为稳定,这与2001—2010年期间整体不显著增加趋势的研究结论基本一致[39]；对整个青藏高原的研究表明,约占77.5%的区域地表反照率变化趋势统计显著,仅有20.1%(2.4%)的区域呈显著降低(增加)趋势[39]。三江源区域地表反照率年际变化的主导因子是植被因子,其中NDVI对地表反照率权重值最大占全区比接近全区面积的一半,且NDVI对地表反照率影响最大的区域与地表反照率高于0.15的区域(东北部和中西部)相重合,这些区域主要是沙地等低植被覆盖度的未利用地,即与三江源区域地表反照率与NDVI负相关的结论一致[39],而降水、最高气温和最低气温对地表反照率的影响区域分布相对较零散。

除植被和气候因子外,短波反照率的变化还受云覆盖和积雪覆盖及其变化的影响[44]。本研究对反照率数据质量(QA)做了控制处理,把填充值在25%以上的(QA=4、5)像元用临近像元数值平均值代替,进而达到去云处理作用[25]；对于地表积雪对地表反照率的影响[45],本研究选择了6—8月作为生长季,以减弱积雪对地表反照率的影响。已往研究指出三江源部分区域存在多年冻土[46],而对此还需要进一步开展观测和机理模拟研究和分析。

本研究进一步分析了三江源区域植被覆盖度和地表反照率变化的地表辐射温度效应,表明植被恢复可能会导致辐射增温效应,从而对区域气候变暖产生正反馈。但值得注意的是,本研究是基于卫星遥感观测数据的统计分析,在未考虑短波辐射本身年际变化前提下,仅仅计算了由地表反照率变化导致的辐射温度效应,所得到的结论值得参考和进一步思考。当然,今后研究需进一步考虑由地表净辐射增加后,地表辐射在潜热和显热间的分配,这就需要基于陆表过程模式的模拟和多源数据整合的方法从机制过程方面开展研究,进而更准确地确定其辐射温度效应。

5 结论

三江源生长季地表反照率在空间上呈现从东南到西北上升的分布特征,空间上差异性明显,全区平均地表反照率为0.163±0.027,其中76.4%区域的值集中在0.12—0.18。2001—2018年生长季三江源地表反照率年际变化趋势范围为-0.010—0.007,显著下降区域占总面积的8.4%,显著上升区域占总面积1.9%；总体上三江源地表反照率以每10年(0.152±0.763)%速率不显著下降,而植被是这种变化的主导因素。三江源生长季植被覆盖度增加整体上导致辐射增温效应,升温(0.11±0.42)℃。由植被覆盖度的增大导致的辐射增温占总面积的32.6%,导致的辐射降温占总面积的24.5%,同时覆盖度下降导致的辐射增温效应占总面积的21.7%,导致的辐射降温效应占总面积的21.2%。且草地覆盖度增加,引起地表反照率下降,导致辐射增温效应；而草地覆盖度降低,引起地表反照率上升,导致辐射降温效应。今后需基于陆表过程模式模拟和多源数据整合方法从机制和过程方面开展研究,并考虑积雪覆盖、永久冻土等对地表辐射收支的影响；本研究基于卫星遥感数据的统计分析,有助于对草地退化与恢复过程中的生物地球物理效应的理解和进一步思考。