2000—2021 年青藏高原地区地表净辐射的时空变化*

崔玉祥 胡斯勒图 李同文 姬大彬 张 颢 施建成

1（中国科学院空天信息创新研究院 北京 100094）

2（中山大学测绘科学与技术学院 珠海 519080）

3（中国科学院国家空间科学中心 北京 100190）

0 引言

地表净辐射是指地球表面的入射和出射辐射之间的差值，包括长波和短波辐射，是衡量地表辐射收支总量的一个重要参数。从地表能量平衡的角度看，地表净辐射可以转化为热通量，包括土壤热通量、感热通量和潜热通量等，地表净辐射是影响整个地气能量交换和再分配的重要参数[1]。同时，作为研究地气间相互作用的关键因素，地表净辐射是地表蒸散、大气热量和水循环的主要驱动力，其时空变化直接影响着水分和能量交换[2-5]，用于气候监测、天气预报和农业气象等领域。

青藏高原地处亚洲大陆中心地带，平均海拔超过4 km，是地球上海拔最高且地形复杂的高原之一，素有地球第三极之称[6]。青藏高原与大气圈、生物圈和水圈等圈层耦合作用，对高原及其周边地区气候格局、亚洲季风形成，乃至北半球大气环流、全球的气候变化和极端天气形成有巨大影响[7-9]。因此，分析研究青藏高原地区地表净辐射的时空变化具有重大意义。

对青藏高原地区地表净辐射的时空变化特征已有很多研究，并取得了一定成果[10-13]。Ciren 等[14]采用青藏高原羊八井地区2009—2010 年期间辐射观测数据，研究发现羊八井地区地表净辐射最高值出现在6 月，最低值出现在12 月。Li 等[15]基于狮泉河地气交换综合观测站和狮泉河气象站辐射观测数据，研究了阿里地区净辐射通量的日循环、日变化、季节变化和年际变化特征。结果表明，阿里地区净辐射平均日变化峰值出现在14:30 LT—15:00 LT，净辐射月均值最大值出现在8 月，且近20 年来净辐射的平均日变化峰值和月均值最大值出现时间延后现象。Zhao 等[16]利用青藏高原羊八井观测站2020—2021 年的太阳短波辐射数据和地气长波辐射数据，发现羊八井地区地表净辐射呈现显著的单峰月际变化特征，月平均值6 月最高，12 月最低。以上研究为青藏高原地区地表净辐射的研究奠定了重要基础，但是大部分地表净辐射变化规律仅从时间或者空间一个维度分析，没有在时空维度综合考虑。另外，大多数研究只在站点位置或者某一局部区域进行研究，涉及青藏高原全域的较少。卫星遥感数据相较于地基站点数据，拥有大面积、长时序等显著优势，特别有利于大范围长时间序列的地表净辐射估算问题的研究[17,18]。

本文利用CERES 获取青藏高原地区2000 年3 月1 日至2022 年2 月28 日的逐日地表辐射各分量数据，计算地表净辐射数据，并结合经验正交函数分析法（Empirical Orthogonal Function，EOF）、Theil-Sen Median 趋势分析和Mann-Kendall 分析法，深入探究青藏高原净辐射通量的时空变化特征，以期为青藏高原地区地气间的能量交换和再分配探究以及长期气候变化研究等提供科学依据。

1 研究区域

青藏高原是世界上海拔最高、面积最大的高原，位于中国西南部。图1 表明研究区** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/东部地势较低，西部地势较高，平均海拔超过4000 m。总体而言，该地区地势高，大气稀薄，太阳辐射能量更容易到达地面，是全球辐射最强的地区之一。本文在基于站点数据验证CERES 卫星数据有效性基础上，结合多种分析方法探究青藏高原地区地表净辐射的时空变化特征。

图1 青藏高原地区高程及地基站点分布Fig. 1 Elevation map of Qinghai-Tibet Plateau

2 数据与方法

2.1 数据来源及预处理

研究采用的数据为2000 年3 月1 日至2022 年2 月28 日CERES 的SYN1 deg-Level 3 数据集*** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/逐日地表辐射各分量数据--地表上行短波辐射、地表下行短波辐射、地表上行长波辐射、地表下行长波辐射，由公式1 计算得到地表净辐射，空间分辨率为1°×1°，数据单位为W·m-2。

使用的地基站点数据来自于中国气象数据网站**** https://cstr.cn/18406.11.Geogra.tpdc.270099** https://ceres.larc.nasa.gov/data/*** http://data.cma.cn/，基于如表1 所示的青藏高原地区所有地基站点（拉萨、西宁、那曲、昌都、玉树），获取2015 年12 月31 日18:00 LT 到2016 年12 月31 日17:00 LT 逐小时地表下行短波辐射，用于评估CERES 辐射参数的数据质量，地表净辐射

表1 青藏高原地区拉萨、西宁、那曲、昌都和玉树5 个地基站点参数Table 1 Parameters of five ground stations in Qinghai-Tibet Plateau region,including Lhasa, Xining, Naqu, Changdu and Yushu

式中，Rsd为地表下行短波辐射；Rld为地表下行长波辐射；Rsu为地表上行短波辐射；Rlu为地表上行长波辐射，单位为W·m-2。

2.2 有效性验证

卫星数据比地基站点观测数据空间覆盖能力强，更利于反映地表净辐射的时空分布和分配过程，但是地基站点的观测具有更高的精度。因此，利用5 个地基站点数据对卫星数据进行精度验证，相关数据采用2015—2016 年共8784 h 地表下行辐射数据，以地基站点值为参考值，分别计算决定系数R2和平均绝对误差E（单位W·m-2），即

其中，N为样本数量，X和X分别为地基站点观测真值及其平均值；Y和Y分别为CERES 卫星数值及其平均值。图2 给出了CERES 卫星数据有效性验证结果，五个站点的R2分别为0.92，0.87，0.88，0.85 和0.87，同时，E分别为44.41，55.51，45.03，56.80 和53.64 W·m-2。上述结果表明，CERES 数据相对稳定且有效性较高，能够满足研究需求。

图2 拉萨（a）、西宁（b）、那曲（c）、昌都（d）和玉树（e）站点CERES地表下行短波辐射和地基站点数据验证Fig. 2 Validation of CERES surface downward shortwave radiation data with ground-based station data at sites in Lhasa （a）, Xining （b）, Naqu （c）, Changdu （d） and Yushu （e）

2.3 研究方法

2.3.1 趋势分析法

Theil-Sen Median 趋势分析[19,20]，是一种稳健性强且不依赖参数的趋势计算方法，该算法摆脱了测量值服从正态分布的缺点，并且在缺失值和异常值影响下依然有效，因而在长时序数据趋势分析中得到了广泛应用。具体计算过程如下：

式中，Median 为中位数运算，Rni和Rnj分别代表第i和j年的地表净辐射，i和j分别为第i和j年，且j＞i；S表征地表净辐射的变化趋势。当S ＞0，表明该研究时段内地表净辐射有加强趋势；当S=0，表明该研究时段内地表净辐射无变化趋势；当S ＜0，则表明该研究时段内地表净辐射有减弱趋势。

2.3.2 Mann-Kendall 检验

Mann-Kendall 显著性检验法[21,22]，是一种不受异常值影响且操作简单的检验长时序变化趋势显著性的方法，已经在水文、气象等领域取得成功应用。本文采用该方法来检验地表净辐射变化的显著性。具体算法如下：

式中，1 ≤i ＜j≤n，n为时间序列长度；fji为Rnj-Rni的差值函数；Svar为S的方差，Z反映地表净辐射变化趋势的显著性。当|Z|≥2.58时，表示地表净辐射变化趋势已通过显著性水平α=0.01的显著性检验；当1.96＜|Z|＜2.58时，则说明地表净辐射变化趋势通过了显著性水平α=0.05的显著性检验；当|Z|≤1.96表示地表净辐射变化趋势未能通过显著性水平α=0.05的显著性检验。

Theil-Sen Median 趋势分析对异常值不敏感，可用于检测数据中是否存在某种趋势，而Mann-Kendall 显著性检验法可以验证被检验的时间序列的趋势性是否通过某置信度下的显著性检验。结合Theil-Sen Median 趋势分析与Mann-Kendall 显著性检测（见表2[22]），即可更有效反应地表净辐射的变化趋势。

表2 地表净辐射变化特征划分标准Table 2 Criteria for classifying changes in surface net radiation

2.3.3 EOF 分析

EOF 分析[23-25]是一种常用的数据分析方法，常用于对海气耦合系统、气象、气候等领域中的大规模环境数据进行降维、分解和重构。其不仅可以降低数据维度，并且能尽可能地保留原始数据的空间和时间变化特征。特征向量所对应的空间样本可以展现其空间分布特点；而主成分则表现时间变化，揭示了随时间变化的不同空间模态的重要性权重。

对青藏高原地区2000 年3 月到2022 年2 月共264 个月逐像素（共299 个像素点）地表净辐射处理成地表净辐射距平值Xm×n，基于EOF 分析求解特征根及相应特征向量，计算各空间模态的方差贡献率并进行North 检验，具体计算公式如下：

式中，n=264，m=299，Xm×n为地表净辐射距平值；Cm×m为协方差矩阵；λ1，λ2,...,λm为协方差矩阵Cm×m对应的特征根；Vm×m为特征向量；Em×m为m×m的对角矩阵，Pm×n为时间系数，A为方差贡献率，R为North 检验容许的误差范围。

3 结果与讨论

3.1 年平均地表净辐射的分布特征

图3 给出了青藏高原地区逐年平均地表净辐射空间分布，从空间分布来看，青藏高原年平均地表净辐射的空间分布特征首先表现为纬度水平条带特征，即年平均地表净辐射沿纬度从南到北逐渐降低，最高值出现在喜马拉雅山脉和横断山脉一线。其次，从经度方向来看，年平均地表净辐射从西向东大致呈现先减弱，后增强，再减弱，再增强的高—低—高—低—高模式，这主要与青藏高原的地势有关[26]，青藏高原西部有着平均海拔6000 m 以上的喜马拉雅山脉，唐古拉山脉横卧青藏高原中部，而横断山脉盘亘于青藏高原东部。

从时间变化来看，年平均地表净辐射2000—2016 年大都在100 W·m-2以上，2017—2021 也大都在90 W·m-2以上；青藏高原中部地区2000—2016 年处于90～100 W·m-2之间，2017—2021 处于80～90 W·m-2之间，最低值为塔里木盆地，常年处于50～60 W·m-2之间。其中，值得注意的是，青藏高原大部分地区在2000—2016 年年平均地表净辐射基本保持平稳状态，而在2016—2017 年间发生突变，2017—2021 年年平均地表净辐射也基本维持不变。

3.2 结合趋势分析法与M-K 检验法的地表净辐射趋势分析

图4 给出了2000—2021 年青藏高原地区地表净辐射变化分布，分析可知，2000—2021 年青藏高原地表净辐射整体处于基本不变和不显著减少状态。其中，帕米尔高原、昆仑山、藏北高原、柴达木盆地、藏东南等地区地表净辐射基本不变；冈底斯山脉、唐古拉山脉、昆仑山脉、青海湖以及青藏高原中东部均呈现不显著减少趋势；喜马拉雅山脉附近显著减少；藏北高原附近出现零星不显著增加区域。结合图5 地表辐射各分量变化趋势，可以看到喜马拉雅山脉附近地表下行辐射呈稳定不变或者不显著减少趋势，而地表上行辐射却出现不显著增加，甚至显著增加趋势，因此，此地地表净辐射呈显著减少和不显著减少状态。此外，青藏高原地表上行长波辐射全域呈现基本不显著增加趋势，这在一定程度上表明了青藏高原大部分地区地表温度总体呈增长趋势，响应了全球变暖的气候特征，这与Gu 等[6]探究青藏高原地表辐射通量的气候特征结论高度吻合。

图4 2000—2021 年青藏高原地区地表净辐射变化分布Fig. 4 Distribution of surface net radiation changes in the Qinghai-Tibet Plateau region from 2000 to 2021

图6 给出了2000—2021 年青藏高原地区地表净辐射季节、年均变化趋势以及云覆盖面积分数变化趋势。由图6 分析可知，从年内变化来看（见图6 a），青藏高原地表净辐射有明显季节变化特征。2022 年来青藏高原地表净辐射均呈现出夏季＞春季＞秋季＞冬季的规律。原因可能在于夏季青藏高原地区太阳高度角大，太阳总辐射相应较大，对地表净辐射的增量远大于因大气中水气含量较高以及云量相对较高对地表直接辐射而产生的减少量[27]。这也与Gao 等[1]关于青藏区地表净辐射各个季节总体趋势保持一致。

图6 2000—2021 年青藏高原地区地表净辐射季节趋势（a）、地表净辐射年均变化趋势（b）和云覆盖面积分数变化趋势（c）Fig. 6 Seasonal trends （a）, annual changes （b）, and trends in cloud area fraction （c） for surface net radiation in the Qinghai-Tibet Plateau region from 2000 to 2021

从图6（b） 所示的年际变化来看，2000—2021 年青藏高原地表净辐射在2016 年发生显著变化，突降约6 W·m-2。表现为2000—2016 年和2017—2021年，两个时间段内地表净辐射各季均值变动幅度不大，但是2017 年各季节地表净辐射均明显低于2016 年数值。如图6（c） 所示，查阅资料后发现青藏高原地区2000—2016 年云覆盖度年均值基本维持在47.01%， 而在2017—2020 年，这一数据达到55.87%，提升约18.75%，地表净辐射与云覆盖分数突变一致性较强，进一步利用皮尔逊相关性检验发现：两者相关性系数为-0.90（p＜0.001），具有显著负相关性，猜测可能是云覆盖度增加会阻挡更多太阳辐射，使得地表收到的太阳辐射减少，最终导致地表净辐射下降。

3.3 EOF 方法地表净辐射通量变化分析

表3 给出了前5 个特征向量的特征根、方差贡献率以及特征根误差范围，从表3 可以看出，第一空间模态贡献最大，为75.55%，第二方差贡献急剧减少，仅为7.58%，其后，方差贡献率均较小，前2 个特征向量的累积方差贡献率已经达到83.13%，并且都通过了North 显著性检验，因此前2 个特征向量基本可以反映青藏高原2000—2021 年地表净辐射的分布类型。

表3 青藏高原地区地表净辐射时空场EOF 分解前5 个模态特征值与方差贡献率Table 3 Characteristic values and variances contribution rate of the first five EOF modes for the spatial-temporal pattern of surface net radiation in the Qinghai-Tibet Plateau region

第1 模态特征向量的方差贡献率为75.55%，远高于其他模态的方差贡献率，可以反映青藏高原地表净辐射分布最重要的典型场，代表区域分布的平均特征。从图7（a）所示的第一特征向量空间分布特征来分析，第1 模态在青藏高原地区的空间系数全部为正值，这表明2000—2022 年青藏高原地表净辐射变化趋势具有整体一致性，即青藏高原地区地表净辐射均呈现出普遍较强或者普遍较弱的特征。同时，第1 模态系数大致呈以喀喇昆仑山脉为低值中心向四周扩大的趋势，高值区位于喀喇昆仑山脉、喜马拉雅山脉东南部、巴颜克拉山脉、横断山脉附近，表明这些区域地表净辐射变化大。在图7（b）所示的第1 模态时间系数中，110 个月为正值，154 个月为负值，第一时间系数均在零值上下周期浮动，且呈现准正弦振动，年周期变化明显，这主要是受到太阳高度角的周期变化影响。峰值主要出现在每年夏季的7 月和8 月，谷值出现在每年冬季的2 月和3 月，恰好与青藏高原夏冬季节气候特征吻合。

图7 2000—2021 年青藏高原地区地表净辐射第1 模态（a）特征向量分布及对应月时间分布（b） ，第2 模态（c）特征向量分布及对应月时间分布（d）Fig. 7 Distribution of the first mode （a） and corresponding monthly time distribution （b） of surface net radiation, distribution of the second mode （c） and corresponding monthly time distribution （d）in the Qinghai-Tibet Plateau region from 2000 to 2021

第2 模态特征向量的方差贡献率为7.58%，也属于青藏高原地区较为典型的地表净辐射表现形式。从图7（c）可以看出，这种分布格局大致以喀喇昆仑山脉和喜马拉雅山为高值中心向四周减小，青藏高原中部大范围为正值，而在兴都库什山脉、祁连山脉和喜马拉雅山脉南麓出现负值中心，大体呈现西南-东北反向分布模式。除喜马拉雅山脉外，特征向量值从西南向东北依次减小，反映出青藏高原地区地表净辐射变化也是从西南向东北依次递减。由图7（d）分析表明，第2 模态对应的时间系数也有一定的年周期变化，这可能是地表综合作用的结果。其中，时间系数为负值有123 个月，地表净辐射位相与图7（c） 相反，空间分布呈现喀喇昆仑山脉和喜马拉雅山低值、而周围高值的情形，典型时间为2004 年8 月。时间系数为正的月份有141 个月，呈中间高四周低分布，典型时间为2002 年12 月。

4 结论

青藏高原高原气候复杂且大气圈、生物圈和水圈等圈层耦合作用，导致其地表净辐射的时空分布和变化复杂多样。基于2000 年3 月至2022 年2 月共22 年的CERES 卫星数据，采用Theil-Sen Median 趋势分析、Mann-Kendall 检验以及EOF 分析多种方法，分析22 年间青藏高原地表辐射能量平衡的时空变化规律，得出结论如下。

（1）青藏高原地表净辐射通量总体呈现南高北低的分布特征，最高值在喜马拉雅山脉和横断山脉一线，最低值在柴达木盆地。其中，青藏高原大部分地区在2000—2016 年地表净辐射通量明显更强，而2017—2021 年地表净辐射通量明显降低。

（2）Theil-Sen Median 趋势分析和Mann-Kendall检验表明，青藏高原地表净辐射通量整体处于基本不变和不显著减少状态，但是喜马拉雅山脉附近地区的地表净辐射通量显著减少。其他年份地表净辐射通量大体处于平稳状态，2016—2017 年，出现近6 W·m-2的突变，这可能是对云覆盖分数增加的响应。

（3）EOF 分解结果表明，青藏高原地表净辐射通量场主要有两种模态。第1 模态方差贡献率为75.55%，可以代表地表净辐射通量变化的总体空间特征，地表净辐射通量变化趋势具有高度一致性，青藏高原南部地表净辐射通量波动程度远高于其北部，对应的时间系数显示，年周期变化明显，第一峰值出现在2018 年8 月，第一谷值出现在2001 年3 月；第2 模态大体呈现西南-东北反向分布模式，除喜马拉雅山脉外，特征向量值从西南向东北依次减小，其时间系数也有一定的年周期变化，这可能是地表综合作用的结果。