基于多源数据集估算缺资料地区地表净辐射及其时空变化特征

张晓龙，沈 冰，黄领梅，权 全，莫淑红，梁晓燕

(西安理工大学 省部共建西北旱区生态水利国家重点实验室，陕西 西安 710048)

地表净辐射(Rn)是影响陆-气能量交换和再分配过程的重要参数[1-2]，是气候变化的重要驱动力，被广泛应用于气候监测、天气预报和农业气象研究，是当前全球气候变化研究领域的热点之一[3]。同时，Rn也是构建各类生态模式的重要参数之一[4-5]，尤其是在生态系统的蒸散过程中起到非常重要的作用[6]。在全球气候变化的背景下，研究Rn时空变化特征，对于了解过去和当前的气候状态以及对未来的气候变化、蒸散发的估算、植物生长发育过程、生态系统生物量的形成与累积等研究具有重要意义。目前获取地表净辐射的方法主要是将地面实测数据按照某种规则进行空间插值，或通过经验统计方法建立基于地面辐射实测数据和气象资料估算地表辐射收支的经验-统计模型，通过反演得到一定范围内地表净辐射[2, 7-8]。这两种方法虽然在点尺度上估算精度较高，但由于气象站数量有限且分布不均，地表覆被类型复杂，实时气象参数(如气溶胶，云的含量等)难以获取，因此不能满足大范围净辐射估算的需求[12-13]。遥感反演法应是区域由点及面较佳的技术手段。近年来，MODIS数据、Landsat系列数据等被广泛应用于地表辐射收支的研究中[9-11]。

我国气象系统仅有50个地表净辐射观测台站，且均建于1993年以后。由于站点少，资料序列短，开展的相应研究也较少。青藏高原生态环境脆弱，被称为气候的“放大器”[12-13]，但常规气象观测台站分布稀少且极不均匀，为缺资料地区空间范围上净辐射估算带来极大阻碍，远满足不了科研和业务的需求。当前主流遥感反演产品使用地面站点较少，特别是在中国寒区、旱区[18]，数据精度还存在一定误差。利用地面气象站数据对遥感产品进行数据融合得到的再分析常规气象要素的驱动数据集不仅可以满足空间分布的要求，而且在数据精度上也可以大大提高，为进一步研究提供极大的便利。地表反照率决定了太阳辐射能量被地表反射与吸收的比例，是气候系统的重要驱动因子之一[13]。地表发照率的分布与土地利用类型、植被状况、水分状况、积雪、地形等因素密切相关，在时间和空间上呈现高度异质性。因此，本文以黄河源区为例，利用ITPCAS气象要素驱动数据集和GLASS地表反照率数据集，取用联合国粮农组织(FAO)1998年推荐用于计算参考作物蒸散发Penman-Monteith公式估算缺资料地区长时间序列逐日Rn，在再分析数据集估算精度验证的基础上，用GIS空间分析技术与Mann-Kendall趋势分析方法揭示其时间动态特征与空间演变规律，为进一步研究该地区气候变化及陆地生态系统蒸散的影响提供技术支撑和理论依据。

1 研究区概况

研究区为龙羊峡水库以上的黄河源区，位于32.2°～37.1°N，95.9°～103.4°E之间，属于青藏高原东部(见图1)。

图1 研究区概况及气象站点分布Fig.1 General situation of the study area and distribution of meteorological stations

黄河源区总体为西北高，东南低，阿尼玛卿山为最高点，龙羊峡为最低处，海拔在2 508～6 253 m之间，有冰川、盆地、高山和峡谷等地貌。黄河源区属大陆性高原气候，干湿交替、冷热分明，气压低，温差大，辐射强，光热条件较为丰富。受印度洋季风影响，降水集中在6～9月[5]，年平均降水234.6～839.2 mm，空间分布不均[12]。黄河源区3/4的面积分布着高寒植被，其中高寒草甸和高寒草原面积最大。黄河源区为多年冻土和季节性冻土的过渡区域，多年冻土分布下界平均气温在-2.5℃ ～ -3.5℃范围内。由于气候变暖，黄河源区近年来冻土发生显著退化。

2 数据来源与方法

2.1 数据来源及预处理

再分析气象数据集使用由中国科学院青藏高原研究所开发的ITPCAS(中国区域高时空分辨率地面气象要素)驱动数据集(http://westdc.westgis.ac.cn/)[13]。该数据集以GEWEX-SRB 辐射资料、Princeton 再分析资料、TRMM 降水数据以及GLDAS数据为背景场，并融合中国气象观测数据制作而成。气象要素包括降水量、温度、比湿、风速、大气压力、向下短波和向下长波辐射等，时间分辨率为3 h，空间分辨率为0.1°。该数据集在黄河源区数据采样点1 421个，可满足黄河源区气象要素空间分布的要求(见图1)。本研究中所涉及的日尺度参数是通过3 h时段参数计算得到。为了分析Rn的时空变化特性，利用ArcGIS软件CUBIC重采样技术将0.1°数据集生成空间分辨率为1 km的数据集。气象站数据源自国家气象科学数据共享服务平台(http://data.cma.cn/)。黄河源区再分析数据集数据采样点和气象台站点分布见图1。

地表反照率使用全球陆表参量GLASS数据产品，取值范围0～1，被用于计算短波净辐射，来源与北京师范大学(http://glass-product.bnu.edu.cn)。该产品2000-2012年反照率数据是基于MODIS数据利用AB(Angular Bin)算法开发的，其空间分辨率为1 km，时间分辨率为8天。Liu等利用地面实测数据和MODIS地表反照率数据MOD43对GLASS反照率数据产品进行对比验证，结果显示两者具有很好的一致性，均方根误差小于0.05[14]，所以该数据集可以直接使用。因数据源的时间限制，本研究计算得到2000-2012年白天地表反照率每一日的多年平均值，然后线性插值得到年内每天的反照率数据集。

DEM数据源自地理空间数据云(http://www.gscloud.cn/)中90 m 分辨率SRTM 产品。为了确保数据空间尺度的一致性，将90 m DEM数据利用ArcGIS软件重采样生成1 km分辨率的数据。

2.2 地表净辐射估算方法

本研究采用FAO推荐的标准方法[15]，由净短波辐射Rns和净长波辐射Rnl之差得到Rn(MJ·m-2·d-1)：

Rn=Rns-Rnl

(1)

Rns=(1-α)Rs

(2)

(3)

Rso=(0.75+2z×10-5)Ra

(4)

(5)

cos(φ)cos(δ)sin(ωs)

(6)

式中：α为地表反照率；Rs和Rso为向下短波辐射和净空向下短波辐射(MJ·m-2·d-1)；Ra为地外辐射(MJ·m-2·d-1)，n为日照时数，N为可照时数，z为海拔高度(m)；Tmax,k和Tmin,k分别为日最大、最小绝对温度(K)；σ为Stefan-Boltzman常数(4.903×10-9MJ·m-2·d-1)；ea为实际水汽压(kPa)；as=0.25，as+bs=0.75；Gsc为日辐射常数(0.082 MJ·m-2·min-1)，dr为日地距离，ωs为日落时角，φ为太阳时角，δ为太阳赤纬角。

需要指出的是，基于中国区域高时空分辨率地面气象要素驱动数据集估算Rn_ITP时，Rs直接采用驱动数据集数据。而基于气象台数据估算Rn_m时，Rs采用式(5)计算得到。

2.3 地表净辐射时空变化分析方法

采用Mann-Kendall趋势分析方法对多站点和区域平均长时间序列的Rn进行时间趋势变化检验和突变分析。基于秩的Mann-Kendall趋势分析方法是一种非参数统计检验法，与传统参数方法相比，样本可不遵从特定分布，亦不受个别异常值干扰，可较客观地反映样本序列变化趋势，该方法被广泛应用于序列突变分析和趋势检验。

利用线性回归分析可得到每个栅格的Rn在一定时间序列的变化趋势和强度[6]。本文通过IDL程序得到每个栅格多年变化趋势b，进而得到黄河源区Rn年际变化率空间分布图。计算公式为：

(7)

式中：b为变化趋势，如果b<0则表示减少趋势，反之b>0则表示增加趋势；n为时间长度；i为年序号；xi为第i年的栅格值。

3 结果与分析

3.1 再分析数据集精度检验

黄河源区处于高寒区域，气象站点稀少，且分布极不均匀。所以使用气象站点进行空间插值的结果并不能较为真实地反应研究区实际情况。

本研究使用ITPCAS驱动数据集估算黄河源区Rn的时空变化。将再分析数据集代入式(1)～(6)得到的1979—2015年逐日Rn_ITP与气象站数据代入式(1)～(6)得到的对应的逐日Rn_m。

为了定量评价Rn的估算结果精度，将7个典型气象站的气象数据估算的Rn_m与对应的1 km气象要素数据集计算的Rn_ITP进行对比分析，Rn_ITP与Rn_m的散点图见图2。虚线为散点的拟合直线，对角线为1∶1等值线。

从图2可以看出整体上Rn_ITP与Rn_m的相关系数R为0.98，其显著性在P<0.01范围内显著，表明散点分布有一定的稳定性，二者有显著的相关性。大部分的点分布在1:1等值线上方，表明Rn_ITP较Rn_m略微偏大。

图2 Rn_ITP与Rn_m的散点图Fig.2 Scatter plots of Rn_ITPvaluesand Rn_mvalues

对逐日Rn的估算精度进行检验的评价指标为：均方根误差(RMSE)、相关系数(R)和误差，在各站点上计算结果见表1，其中区站号为国家气象站统一编号，N为样本数。从表1中可知，Rn_ITP与Rn_m的相关系数在0.98以上，且均通过显著性检验(P<0.01)，各站点RMSE在1.29～1.77之间，值得注意的是，Rn_ITP均比Rn_m偏大12%～19%，主要原因可能是在气象数据估算Rs时外辐射到达地表的比例参数选用默认值(as= 0.25和bs= 0.50)导致ITPCAS驱动数据集中Rs比气象站基于日照时数计算的Rs偏高，也可能是气象采样点和气象站并不完全重合，存在着坡度、海拔和坡向等差异。综上，通过该数据集估算的Rn可以较好的反映区域变化规律，该再分析数据集在黄河源区有较好的适用性。

表1 再分析数据集估算Rn的结果验证

3.2 黄河源区地表净辐射的变化特征

3.2.1 地表净辐射时间变化特征

精度检验后将ITPCAS驱动数据集和GLASS地表反照率数据集代入公式，得到1979―2015年黄河源区逐日Rn。对逐日数据进行处理，得到各月、季、年Rn数据。季节采用国内常用的气象季节划分法，即上年12月至当年2月为冬季、3月至5月为春季、6月至8月为秋季。各月Rn为月内每日平均得到，各季节Rn为相应月份Rn的平均值，年平均Rn为逐月数据平均得到。

图3为1979-2015年黄河源区年、春季、夏季、秋季、冬季平均Rn变化曲线。1979-2015年黄河源区年、春季、夏季、秋季、冬季多年平均Rn分别为9.63、12.11、13.70、7.68和5.02 MJ·m-2·d-1；其最大值分别为10.00(1997年)、12.80(1980年)、14.59(1994年)、8.07(1997年)和5.38(1983年)MJ·m-2·d-1；最小值分别为9.12(2009年)、11.42(1989年)、12.63(2009年)、7.24(1988年)和4.76(2009年)MJ·m-2·d-1。

利用Mann-Kendall趋势分析方法得到黄河源区年、四季Rn的变化趋势和统计值，见表2。从表2可知，黄河源区年、四季Rn均呈下降趋势，每10年变化幅度为-0.067 MJ·m-2·d-1；春、冬季变化幅度最大，每10年变化幅度为-0.075 MJ·m-2·d-1；秋季变化幅度最小，每10年变化幅度为-0.039 MJ·m-2·d-1；夏、秋季没有通过显著性检验。通过Mann-Kendall非参数检验分析研究区年、春季、夏季、秋季、冬季Rn突变点，得到结果：年Rn和冬季Rn序列存在突变点，分别为2004年和2003年左右；而其他季节Rn突变点不明显。这可能与当地气候变化和人类活动有关。2000年以后黄河源区呈现“变暖变湿”的趋势，云量上升，降雨/雪增加，温度也增加，植被NDVI也呈上升趋势。另外三江源自然保护区于2000年成立，实施了一系列的生态保护与建设工程，青海省也实施了人工增雨、治理草原鼠害、退耕/牧还草和生态移民等 “综合性”工程(2004年)。 这些都可能是年Rn和冬季Rn序列在2004年和2003年左右存在突变的原因。

图3 1979-2015年研究区年、春季、夏季、秋季、冬季平均Rn变化曲线Fig.3 Curves of annual and seasonal mean Rn during 1979-2015 in study area

时间序列起始年份终止年份 样本数量Test Z显著性每10年变化幅度/(MJ·m-2·d-1)年Rn1979201537-2.68**-0.067春季Rn1979201537-2.31*-0.075夏季Rn1979201537-1.14-0.074秋季Rn1979201537-1.03-0.039冬季Rn1979201537-3.60***-0.075

注：***通过显著性水平为0.001的检验；**通过显著性水平为0.01的检验；*通过显著性水平为0.05的检验。

图4为黄河源区年内各月平均Rn的变化趋势，可知，区域平均Rn在年内呈单峰型，1～4月上升迅速，5月达到最大值(14.31 MJ·m-2·d-1)，9月份以后Rn开始快速下降，12月份最小(3.34 MJ·m-2·d-1)。通过最大、最小值误差线可知，1979-2015年Rn变化幅度夏季大于冬季，其中7月份Rn变化幅度最大，变化范围为11.79～15.86 MJ·m-2·d-1。

3.2.2 地表净辐射空间变化特征

对1979-2015年黄河源区逐栅格计算得到区域上年、春、夏、秋、冬季Rn和其变化幅度的空间分布及分类面积统计，结果见图5。

从图中可知，年和各季节Rn的空间分布相似，即不同季节对Rn的空间分布影响不大，最大值分布在扎陵、鄂陵湖及龙羊峡等大面积水域附近，最小值分布在中部的阿尼玛卿山附近和共和盆地附近。年和各季节Rn的变化趋势的空间分布也相似，研究区Rn以下降趋势为主，主要分布在研究区西部山区、龙羊峡南部及共和盆地附近。研究区Rn上升趋势主要分布在中部的阿尼玛卿山附近和研究区东部的零星区域，年、春、夏、秋、冬季Rn上升趋势分别占流域的8.4%、15.8%、4.6%、44.1%和0.0%。

3.2.3 海拔对地表净辐射的影响

Rn的空间分布与海拔分布呈现一定的统计规律。将DEM数据重采样到1 km栅格，然后将该DEM分成5类对多年平均的年、春、夏、秋、冬季Rn和其变化比例的进行统计分析，见表3。从表3中可知，Rn在4 000～4 500 m的区域内最大，而在4 500 m以上最小，这是因为4 500 m以上区域存在积雪和冰川，地表反照率高，从而Rn较小。

Rn的变化幅度中除了秋季在3 500～4 000 m区域内呈0.12%的上升趋势外，其余皆呈现下降趋势，年、春季和夏季Rn下降幅度最明显的区域为4 500 m以上区域，秋季和冬季下降幅度最明显的区域为3 000 m以下区域。

图4 研究区Rn年内变化 Fig.4 Changes of monthly mean Rn in study area (The bars indicate the maximum and minimum monthly mean values)

DEM分类面积比例Rn/(MJ·m-2·d-1)年春季夏季秋季冬季每10年Rn变化比例/%年春季夏季秋季冬季3 000 m以下2.6%9.3611.8213.707.344.59-0.66-0.21-0.66-0.52-2.013 000～3 500 m16.0%9.3911.9513.297.474.85-0.78-0.67-0.79-0.31-1.723 500～4 000 m26.1%9.7312.3713.307.885.36-0.57-0.68-0.610.12-1.274 000～4 500 m37.4%9.8212.3213.927.965.07-0.79-0.80-0.84-0.13-1.644 500 m以上17.9%0.921.131.400.700.46-0.96-0.83-1.07-0.31-1.84

4 讨 论

Rn的时空变化改变了区域地表天气和气候系统的多样性、地表蒸散的速率及不同区域生态系统的结构和功能。刘新安等[1]对几种推算地表净辐射的方法进行了评价；陈征等[3]、叶晶等[16]利用遥感资料推算了个别地区的地表净辐射，高扬子等[2]计算了全国699个站点的地表净辐射，并分析了各分量的空间分布特征，对中国地表净辐射做了较为系统的分析。

有研究表明，1980年以来全球Rn具有普遍的下降趋势[17]， 国内近50年来站点平均Rn在年、季节均呈现较明显的下降趋势，每10年均降幅为0.74 W/m2(0.064 MJ·m-2·d-1)[2]，本研究得到与此结论一致。

影响特定区域Rn时空变化特征的因素主要有天文辐射、云、气溶胶、温度、土壤湿度和土壤植被等[2]。

有研究表明，云量不是造成中国区域短波辐射下降的主要因素，大气中持续增多的气溶胶是主导因素[18]。大气气溶胶主要来自于化石燃料的燃烧、工业生产等人为活动及火山爆发等自然现象。

自上世纪80年代以来，随着经济快速发展，化石燃料的大量燃烧造成大气气溶胶的快速增加，这可能是区域乃至全球Rn持续下降的主要原因。

有研究表明，青藏高原积雪变化具有确定的长期增加的趋势[19]，从而导致的青藏高原冬季净辐射的显著降低，另外近年来由于环境污染和气候变化等原因导致冬季的雾霾天气偏低，这也可能造成冬季Rn显著下降的原因之一。

图5 1979-2015年研究区年和四季平均Rn以及对应的每10年Rn变化幅度的空间分布Fig.5 Distribution of mean Rn values and its change rates every 10 years during 1979-2015 in study area

5 结 论

本文以黄河源区为例，利用ITPCAS气象要素驱动数据集和GLASS地表反照率数据集，基于FAO1998年推荐的Penman-Monteith公式估算缺资料地区长时间序列逐日Rn，在再分析数据集估算精度验证的基础上，用 GIS 空间分析技术与Mann-Kendall趋势分析方法揭示其时间动态特征与空间演变规律。本文得到以下主要结论。

1) 各站点Rn_ITP与Rn_m的相关系数均在0.98以上，且均通过显著性检验(P<0.01)，各站点RMSE在1.29～1.77之间，表明ITPCAS气象要素驱动数据集在黄河源区有较好的适用性，通过该数据集估算的Rn可以较好地反映区域变化规律，该数据集可用于中国西北高寒山区缺资料地区。

2) 在1979-2015年期间年、春、夏、秋、冬季区域平均Rn的每10年下降幅度分别为-0.067、-0.075、-0.074、-0.039和-0.075 MJ·m-2·d-1，其中夏、秋季Rn趋势没有通过显著性检验；年Rn的变化趋势在2004年左右存在突变点，冬季Rn的变化趋势在2003年左右存在突变点，其他季节突变点不明显；区域平均Rn在年内呈单峰型，1～4月上升迅速，5月达到最大值(14.31 MJ·m-2·d-1)，9月份以后Rn开始快速下降，12月份最小(3.34 MJ·m-2·d-1)。

3) 年和各季节Rn值及变化趋势的空间分布相似；年和各季节Rn值最大值分布在扎陵、鄂陵湖及龙羊峡等大面积水域附近，最小值分布在中部的阿尼玛卿山附近和共和盆地附近；研究区91.6%的面积的年平均Rn呈下降趋势，在季节尺度中，冬季Rn下降最明显，几乎100%的区域处于下降趋势之中；秋季Rn下降最不明显(约55.9%)。

4) 地形对Rn的空间分布产生一定的影响，Rn在4 000～4 500 m的区域内最大，而在4 500 m以上最小；Rn的变化幅度中除了秋季在3 500～4 000 m区域内呈0.12%的上升趋势外，其余皆呈现下降趋势；年、春季和夏季Rn下降幅度最明显的区域为4 500 m以上区域，秋季和冬季下降幅度最明显的区域为3 000 m以下区域。

本文对长时间序列的逐日逐栅格Rn进行计算，并进行时空趋势变化分析，为无资料地区今后开展地表辐射研究、气候变化预测、生态系统蒸散发过程和机理研究具有重要意义。由于缺少实测的Rn观测数据，在计算精度上存在一定的不确定性。另外本文未定量分析造成黄河源区Rn下降的主要驱动因子，Rn的影响过程与机理需在今后研究工作中进一步探讨。