1 引言
IPCC第六次评估报告指出， 2011 -2020年全球平均温度相比工业化革命前升高了1.09 ℃， 且未来全球每个地区都面临更显著的气候变化， 极端事件的发生频率也将增加（IPCC， 2021）。青藏高原具有特殊和复杂的地理特征（姚檀栋等， 2017）， 被称为“世界屋脊”“第三极”， 该地区生态环境脆弱， 极易受到气候变化的影响， 是气候变化的敏感和脆弱区（秦大河和丁永建， 2009； 陈虹举等， 2021； 樊威伟等， 2021）。目前， 青藏高原气候变化研究已成为国内外关注的热点（杨耀先等， 2022）。此外， 青藏高原地势高， 风速较同纬度其他地区大， 但空气密度较低， 长期以来其风资源被认为不利于开发利用。为更好地评估青藏高原地区风速的长期变化特征， 及其为双碳背景下的风资源发展规划等提供重要参考， 有必要深入分析地表风速的变化及其规律， 获得定性的结论。
以往基于观测事实的研究表明， 最近几十年全球多个地区的地表风速均呈现下降趋势， 如北美洲、 南美洲、 欧洲和亚洲等地（Pryor et al， 2009； Jiang et al， 2010； Wan et al， 2010； Guo et al， 2011； McVicar et al， 2012； Hueging et al， 2013； Kim and Paik， 2015； Li et al， 2020a； Wu and Shi， 2021）。同样， 针对中国不同区域地表风速变化的研究结果也显示， 从20世纪中期至21世纪初期， 整个中国及其大部分地区地表风速明显下降（Jiang et al， 2010； Fu et al， 2011； Guo et al， 2011； Lin et al， 2013； Wu and Shi， 2021）。此外， 一些研究也指出， 近年来， 中国地区地表风速有所增加（Lin et al， 2013； 张志斌等， 2014； Li et al， 2020a； Yang et al， 2021）， 青藏高原增加的趋势也较为明显（Wu and Shi， 2021）。然而， 由于青藏高原地面气象观测站点主要集中在东部地区， 西部地区站点稀少， 数据缺乏， 难以为青藏高原地区风速变化给出全面和准确的结论。此外， 目前高原地区地表风速的研究较集中于1960年代至2010年之间， 很少分析至2020年。因此， 有必要基于高分辨率格点化观测数据以及再分析等多种资料对青藏高原1961 -2020年的风速变化进行深入分析。另外， 除了平均风速外， 关于不同等级风速概率变化的认识对于地表通量估算、 区域风险评估以及风电气候学中的许多应用都至关重要（He et al， 2010， 2012）， 如强风容易对建筑物和交通造成损害， 而低速风则有利于空气中污染物的增加等。然而， 目前很少有研究探讨1961 -2020年期间青藏高原不同等级风速的概率变化。
除了观测事实研究之外， 利用数值模式对风速变化的研究也广泛开展， 其中给出了许多检测和评估的结果， 为中国区域气候变化的未来预估奠定了基础。基于全球海气耦合模式对中国地区风速变化的模式评估结果表明， 耦合模式比较计划第五阶段（CMIP5）和第六阶段（CMIP6）的所有模式都倾向于低估中国地表风速的年际变化， 无法再现观测中的下降趋势（Jiang et al， 2017； 江滢等， 2010， 2018； Wu et al， 2021）。研究还指出， 尽管CMIP6比CMIP5具有更高的空间分辨率和更完善的物理过程， 但其捕捉局部和区域强迫的能力仍然不足， 特别是在中国等地形复杂的地区， 而高分辨率区域气候模式（RCM）在中国区域小尺度强迫和地形方面有很大改进（Gao et al， 2013； Zou and Zhou， 2013； Qian and Leung， 2017； Gutowski et al， 2020； Li et al， 2020b）， 对于青藏高原气候变化研究尤为重要（Ji et al， 2013， 2015； 吉振明， 2018）。然而， 高分辨率区域气候模式动力降尺度对青藏高原地表风速的模拟能力如何， 还鲜少有深入的分析和探讨。
基于以上原因， 本研究使用1961 -2020年青藏高原的观测和再分析的5种风速资料， 分析了近60年来青藏高原月、 季节和年平均地表风速的变化及趋势， 以及不同等级风速概率的变化及趋势。随后， 使用区域气候模式动力降尺度的5个模拟结果及其集合， 检验评估了区域模式对1986 -2005年地表风速月、 季节和年平均、 年内变化和长期变化趋势的模拟能力， 为双碳背景下的青藏高原风资源发展规划以及未来预估、 灾害风险评估等提供依据。
2 数据来源和方法介绍
2.1 观测数据
本研究中使用的逐日10 m风速观测数据为CN05.1（吴佳和高学杰， 2013）， 数据时段为1961年1月1日至2020年12月31日， 水平分辨率为0.25×0.25（纬度×经度）。数据集插值使用的站点数据来源于中国2416个地面观测台站资料， 并首先经过基本的质量控制（Wu et al， 2017）。随后使用“距平逼近方法”插值到网格上（New et al， 2002）。具体而言， 就是首先使用薄板样条插值方法将30年的气候态插值到网格上， 该方法引入协变量子模型可以反映插值气候要素与地形的关系。随后， 使用角距权重方法， 将逐日距平场插值到相同的网格上。最后将以上产生的气候场和距平场相叠加得到最终的数据集。目前CN05.1被广泛应用于中国地表风速的检验评估中（Yu et al， 2019； Chen et al， 2020； Wu et al， 2021； Wu and Shi， 2021）， 表现出与其他风速数据类似的空间分布和时间变化， 同时也在再现中国地区地表风速历史演变方面显示出较好的适用性（Wu and Shi， 2021）。
此外， 选择NCEP-1（Kalnay et al， 1996）、 JRA-55（Kobayashi et al， 2015）和ERA5（Hersbach et al， 2020）再分析资料与CN05.1进行比较， 希望获得基于多源数据的青藏高原风速变化的定性结论， 之所以选择这些数据是因为它们同样覆盖了1961 -2020年的长时间序列。另外注意到， 以往的研究指出， JRA-55在青藏高原地区的适用性比NCEP资料更强（姚慧茹和李栋梁， 2016）， 但其对比时段仅为1971 -2012年， 也有研究指出格点数据插值基于的青藏高原站点较少， 导致插值结果具有一定的不确定性， 因此将1961 -2020年CN05.1格点风速与再分析风速进行对比可望得到其在高原地区适用性的更可靠结论。
风速分级标准使用Zha et al（2016）的方法， 将日平均风速分成6个不同的等级， 分别为无风（0.0~0.2 m·s-1）、 软风（0.3~1.5 m·s-1）、 轻风（1.6~3.3 m·s-1）、 微风（3.4~5.4 m·s-1）、 和风（5.5~7.9 m·s-1）以及大风（≥8 m·s-1）。不同等级风速概率的计算方法为： 某一等级风速在1961 -2020年的日数占总日数的比率（Wu and Shi， 2021）。值得注意的是， 模拟和再分析数据的风速是连续的， 而观测风速是保留一位小数的， 本文其他计算中均采用连续的模拟和再分析风速结果， 只有计算风速等级时， 为了和观测结果一致， 将模拟和再分析数据处理为只保留1位小数， 并进行风速分级。
2.2 模式数据及评估方法
本研究用于动力降尺度的区域气候模式为意大利理论物理中心开发的RegCM4（Giorgi et al， 2012， http： //gforge.ictp.it/gf/project/regcm/）。驱动场选取CMIP5的5个全球模式， 分别为CSIRO-Mk3-6-0、 EC-EARTH、 HadGEM2-ES、 MPI-ESM-MR和NorESM1-M（以下分别称为CdR， EdR， HdR、 MdR和NdR）， 进行了5组水平分辨率为25 km的动力降尺度模拟。全球模式的选择基于CORDEX-CORE框架（International Coordinated Regional Climate Downscaling Experiment Coordinated Output for Regional Evaluations； Gutowski et al， 2016）。模式陆面过程方案选用CLM3.5， 对流参数化方案为Emanuel， 大气辐射传输过程使用NCAR CCM3， 行星边界层方案为Holtslag et al （1990）， 大尺度降水方案采用SUBEX， 关于参数化配置的详细介绍可以参考（Gao et al， 2017）。用于模式评估的时段为1986-2005年。模拟区域为中国及周边地区， 分析区域为青藏高原（图1 ）， 覆盖范围为73°E -104°E， 26°N -39°N， 格点数为125×53。CN05.1使用的青藏高原台站资料为133个， 西藏西北部地区仅有零星站点分布［图1 ， 文中涉及的地图是基于中华人民共和国自然资源部地图技术审查中心标准地图服务系统下载的审图号为GS（2021）7297的中国地图制作， 底图无修改］。
5个模拟的集合平均（ENS）采用等权重方法（Han et al， 2019； Wu et al， 2021）。高原的地表风速在各种资料中均统一使用10 m风速结果， 为便于比较， 使用双线性方法将NCEP-1、 JRA-55和ERA5再分析资料以及动力降尺度模拟结果统一插值到与CN05.1相同的网格上。另外为更好地估计地表风速的变化趋势， 趋势检测采用了Mann-Kendall方法（P<0.05）（Mann， 1945； Kendall， 1948）。
评估指标包括模式与CN05.1的相关系数、 偏差和均方根误差（RMSE）。在评估模式对青藏高原的月、 季和年平均风速的模拟效果时， 还采用了Gleckler et al （2008） 定义的与相对误差类似的相对均方根误差（RMSE′）， 具体定义如下：
首先计算每个模式（或多模式集合）的RMSE。RMSE median是特定时间段所有模式的中值。RMSE′的值越低表示性能越好， 反之亦然。具体地说， 一个模拟在某个时间尺度上的RMSE′小于1表示RMSE小于所有模拟结果在同一时段的RMSE中值。
3 结果分析
3.1 观测到的青藏高原地表风速的变化
图2 1961 -2020年青藏高原区域平均地表风速的变化（相对于1961 -2020年）Fig.2 Time series of the regional mean surface wind speed over Qinghai-Xizang Plateau form 1961 to 2020 （relative to 1961 -2020） |
表1 1961 -2002年和2002 -2020年青藏高原区域平均地表风速变化趋势Table 1 Trends of the regional mean surface wind speed over Qinghai-Xizang Plateau for the periods of 1961 -2002 and 2002 -2020 |
| 数据来源 | 1961 -2002年风速变化/（m·s-1） | 2002 -2020年风速变化/（m·s-1） | ||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| 年平均 | 春季 | 夏季 | 秋季 | 冬季 | 年平均 | 春季 | 夏季 | 秋季 | 冬季 | |
| CN05.1 | -0.23* | -0.18* | -0.22* | -0.18* | -0.22* | 0.12* | 0.03 | 0.12* | 0.19* | 0.15* |
| 站点 | -0.17 * | -0.15 * | -0.17 * | -0.14 * | -0.17 * | 0.11 * | 0.03 | 0.12 * | 0.15 * | 0.15 * |
| NCEP-1 | -0.05* | -0.06 | 0.00 | -0.08 * | -0.08 | 0.16* | -0.01 | 0.06 | 0.30* | 0.16 |
| JRA-55 | -0.08* | -0.06 | -0.07* | -0.07* | -0.12* | -0.06* | -0.19* | -0.02 | 0.07 | -0.07 |
| ERA5 | 0.01 | 0.03 | 0.01 | 0.01 | -0.02 | -0.00 | -0.10 | 0.02 | 0.08 | -0.07 |
*表示变化趋势通过P<0.05的显著性检验， 粗体表示与CN05.1变化符号一致， 斜体表示与CN05.1变化数值差异最小 |
2002 -2020年CN05.1、 站点数据和NCEP年平均风速呈现增加趋势［图2 （a）］且变化趋势均通过P<0.05显著性检验， 分别增加0.12 m·s-1·（10a）-1、 0.11 m·s-1·（10a）-1和0.16 m·s-1·（10a）-1， 而JRA-55则表现为减少0.06 m·s-1每十年（通过P<0.05的显著性检验）， ERA5的年平均风速则无明显变化［0.00 m·s-1·（10a）-1］。从四个季节来看， 秋季平均风速的增加趋势较明显［图2 （d）］， 5种资料均在秋季表现为增加， 但数值上存在差异， NCEP-1增加幅度最大， 为0.30 m·s-1·（10a）-1， CN05.1和站点数据接近， 分别增加0.19 m·s-1·（10a）-1和0.15 m·s-1·（10a）-1， 3者均通过P<0.05的显著性检验， 而JRA-55和ERA5分别增加0.07 m·s-1·（10a）-1和0.08 m·s-1·（10a）-1。春季平均风速除了CN05.1和站点数据表现为弱地增加外， 其他3种再分析资料均表现为减少趋势， 尤其是JRA-55减少趋势达 -0.19 m·s-1·（10a）-1且通过显著性检验（表1 ）。总体来说， 2002 -2020年， CN05.1、 站点和NCEP的变化趋势值更为接近（表1 粗斜体）。
进一步比较多种风速数据在高原地区的年际变率（表2 ）可见， CN05.1与站点数据的相关性达到0.95以上， 与其他风速年际变化的相关系数， 也可以看到， 1961 -2020年， NCEP-1和JRA-55的年平均风速逐年变化与CN05.1最为接近， 相关系数分别为0.48和0.47， 均通过P<0.05的显著性检验。JRA-55的春季、 夏季和冬季平均风速的年际变化与CN05.1最为相似， 相关系数均超过0.5以上（通过P<0.05的显著性检验）， 而NCEP-1的秋季平均风速年际变化与CN05.1的相关系数值略高于JRA-55。注意到， 尽管JRA-55的年际变率与CN05.1在三个季节更为接近， 但其年平均风速年际变率与CN05.1的相关系数并不比NCEP-1高， 其原因是NCEP-1的春、 夏和冬季平均风速变化幅度偏大［图2 （b）， （c）， （e）］， 与CN05.1的相关系数值较JRA-55要小， 标准差也较大， 而年平均则变化幅度正负相抵消使得相关系数较高。总体来看， CN05.1地表风速的年际变率在高原地区的表现和站点数据吻合， 与JRA-55和NCEP-1一致性也较好， 且与JRA-55更为接近。
表2 1961 -2020年CN05.1与NCEP-1、 JRA-55和ERA5风速年际变化的相关系数Table 2 Correlation coefficients of surface wind speed anomaly between the CN05.1 and the NCEP-1、 JRA-55 and ERA5 over Qinghai-Xizang Plateau during 1961 -2020 |
| 数据来源 | 年平均 | 春季 | 夏季 | 秋季 | 冬季 |
|---|---|---|---|---|---|
| 站点观测 | 0.98* | 0.98* | 0.95* | 0.96* | 0.98* |
| NCEP-1 | 0.48* | 0.53* | 0.43* | 0.40* | 0.38* |
| JRA-55 | 0.47* | 0.59* | 0.51* | 0.38* | 0.54* |
| ERA | 0.21 | 0.35* | 0.23 | 0.10 | 0.46* |
*表示变化趋势通过P<0.05的显著性检验 |
随后， 对多种风速资料在青藏高原地区的年内变化也进行了对比（图3 ）。从图3 中可以看到， 1961 -2020年NCEP-1的逐月风速值明显较其他资料要大， 尤其是与ERA5的差异最显著， 除7、 8月外， 其他各月均超过1 m·s-1， 且差异最大的季节为冬季。CN05.1和站点数据的年内变化最为吻合。与JRA-55的风速年内变化也较为一致， 青藏高原地区风速峰值均出现在3 -4月， 但两者低值出现时间存在差异， CN05.1的低值在11 -12月， 而JRA-55则出现在8 -9月， 两个数据差异最大的月份出现在7月和8月， 但其偏差值不超过0.5 m·s-1， 其他月份的差异则在0.1~0.3 m·s-1。可以看到， CN05.1风速的年内变化与以往基于观测数据的结果基本一致， 但数值上普遍高1 m·s-1左右， 这与以往的研究结果一致（姚慧茹和李栋梁， 2016）， 原因可能与青藏高原上的观测台站主要位于东部及南部， 北部至昆仑山北麓等地势高的地区则基本没有观测站点（吴佳和高学杰， 2013）， 这也决定了直接使用较少的观测台站来代表青藏高原平均风速会存在一定程度的低估。
另外， 除了平均风速外， 图4 还给出了不同等级风速概率的变化。从图4 中可以看到， 青藏高原地区不同等级风速变化存在明显差异， 并且总体上再分析资料的年际变率较CN05.1要小， 这与中国区域平均的结果一致（Wu and Shi， 2021）。无风（0.0~0.2 m·s-1）的概率变化较小， 在1961 -2002 年和2002 -2020年以微弱的减少趋势为主（表3 ）。3种以上的数据显示软风（0.3~1.5 m·s-1）和轻风（1.6~3.3 m·s-1）在1961 -2002 年增加， 在2002 -2020年减少， 并且具有显著性变化（通过P<0.05的显著性检验）。微风（3.4~5.4 m·s-1）、 和风（5.5~7.9 m·s-1）的概率则表现为在1961 -2002年减少， 在2002 -2020年增加， 并且两种以上的数据显示出显著性变化。大风（≥8 m·s-1）在1961 -2002年和2002 -2020年均以减少趋势为主。总体来说， CN05.1与NCEP-1、 JRA-55数据的变化趋势比较一致， 相比之下， NCEP-1不同等级风速概率变化正负符号与CN05.1比较接近， 但变化数值上JRA-55与CN05.1更加接近（表3 ）。此外， 注意到， CN05.1与台站数据的不同等级风速概率变化存在明显偏差， 原因主要是插值引起的， 站点在地形高的地方分布稀少， 总体风速值偏低， 而CN05.1经地形订正后风速变大， 从而导致两者不同等级风速概率的显著差异。另外， JRA-55在60年代中期有个明显的突变值， 可能与其使用的卫星资料有关， 而ERA5与其他三种数据的变化差异最大， 在使用时需要予以注意。
图4 1961 -2020年青藏高原区域年平均不同等级风速概率的变化（单位： %， 相对于1961 -2020年）Fig.4 Time series of annual mean surface wind speed probabilities in different grades over Qinghai-Xizang Plateau form 1961 to 2020 （unit： %， relative to 1961 -2020） |
表3 青藏高原不同等级风速概率在1961 -2002年和2002 -2020年的变化趋势Table 3 Trends of the annual mean surface wind speed probabilities in different grades over Qinghai-Xizang Plateau for the periods of 1961 -2002 and 2002 -2020 |
| 数据来源 | 无风/［%·（10a）-1］ | 软风/［%·（10a）-1］ | 轻风/［%·（10a）-1］ | 微风/［%·（10a）-1］ | 和风/［%·（10a）-1］ | 大风/［%·（10a）-1］ | |||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1961 -2002年 | 2002 -2020年 | 1961 -2002年 | 2002 -2020年 | 1961 -2002年 | 2002 -2020年 | 1961 -2002年 | 2002 -2020年 | 1961 -2002年 | 2002 -2020年 | 1961 -2002年 | 2002 - 2020年 | ||
| CN05.1 | -0.02* | -0.02 | 0.60* | -1.21* | 3.77* | -1.54* | 0.12 | 1.50* | -2.85* | 1.74* | -1.03* | -0.18* | |
| 站点观测 | 0.36* | -2.64* | 4.11* | -3.98* | -0.44* | 6.74* | -2.48* | 0.05 | - 1.07 * | -0.24* | -0.26* | - 0.10 * | |
| NCEP-1 | 0.02* | -0.06* | 0.59 * | - 1.33 * | 0.36 | - 1.07 | -0.21* | 0.26 | -0.43* | 1.29* | - 0.29 | 0.92 | |
| JRA-55 | 0.00 | - 0.01 | 0.91* | -0.74 | 0.78 * | 1.00* | -0.77* | 1.22 * | -0.72* | -0.21 | -0.08 | -0.52* | |
| ERA5 | 0.00 | 0.00 | -0.36 | 0.03 | 0.03 | -0.25 | 0.16 | 0.54 | 0.01 | -0.20 | 0.01 | -0.08* | |
*表示变化趋势通过P<0.05的显著性检验， 粗体表示与CN05.1变化符号一致， 斜体表示与CN05.1变化数值差异最小 |
比较多种风速数据的年际变率（表4 ）可以发现， 1961 -2020年NCEP-1、 JRA-55与CN05.1不同等级风速概率的逐年变化最为接近， 尤其是软风、 轻风、 和风及大风（相关系数均通过P<0.05的显著性检验）， 并且NCEP-1与CN05.1的软风和轻风的年际变率更一致， JRA-55与CN05.1的和风和大风的年际变率更一致， 但两者均无法再现CN05.1中无风和微风的年际变率特征。ERA5与CN05.1的轻风和大风年际变率较为接近， 相关系数通过P<0.05的显著性检验， 但其他等级风速概率的年际变化则与CN05.1相关不显著。
表4 1961 -2020年NCEP-1、 JRA-55和ERA5与CN05.1的不同等级风速概率年际变化的相关系数Table 4 Correlation coefficients of annual mean surface wind speed probabilities anomaly in different grades between the CN05.1 and NCEP-1、 JRA-55 and ERA5 over Qinghai-Xizang Plateau during 1961 -2020 |
| 数据来源 | 无风 | 软风 | 轻风 | 微风 | 和风 | 大风 |
|---|---|---|---|---|---|---|
| NCEP-1 | -0.05 | 0.47* | 0.41* | -0.04 | 0.46* | 0.45* |
| JRA-55 | 0.08 | 0.27* | 0.39* | -0.02 | 0.51* | 0.63* |
| ERA | 0.05 | 0.13 | 0.37* | 0.25 | 0.25 | 0.29* |
*表示变化趋势通过P<0.05的显著性检验 |
除此以外， 对于多种格点资料的空间分布差异也进行了对比（图5 ）。从图5 中可以看到， CN05.1中高原地区年平均地表风速呈现西高东低的分布型， 高值区位于高原西北部， 年平均风速值达5 m·s-1以上， 低值区位于高原东部及南部， 数值在1~3 m·s-1［图5 （a）］。NCEP-1的年平均风速［图5 （b）］则基本无法再现出CN05.1中的分布型， JRA-55和ERA5的表现则明显比NCEP-1要好， 与CN05.1中的分布型基本吻合， 但数值上普遍偏小［图5 （c）~（d）］。3种再分析数据与CN05.1年平均风速的相关系数分别为0.08、 0.84和0.55， 其中JRA-55、 ERA5和CN05.1的相关系数均通过P<0.05的显著性检验。总体上， JRA-55与CN05.1的空间分布最为接近。
综上， CN05.1对青藏高原年平均风速和季节风速变化速率、 年际变率、 年内变化、 不同等级风速概率变化以及空间分布等方面都表现出与再分析数据相一致的特征， 并且与站点数据最为吻合， 表明该格点化观测风速数据在高原地区的适用性较好。此外， 由于地形订正插值的考虑， CN05.1风速比站点普遍偏大， 导致两者不同等级风速概率之间的显著差异。而再分析资料与观测之间的偏差可能与其采用不同的数据源和同化方法有关。
3.2 动力降尺度模拟能力评估
随后基于CN05.1格点化观测数据， 对5组水平分辨率为25 km的动力降尺度模拟（CdR， EdR， HdR、 MdR和NdR）的地表风速进行综合评估。首先图6 给出了CN05.1和区域模式集合的1986 -2005年青藏高原平均风速分布及差值。可以看到， CN05.1中的青藏高原年平均地表风速［图6 （a）］大值区位于高原西北部， 数值在6 m·s-1以上， 而高原东部及南部风速较小， 基本在2~4 m·s-1。区域模式集合［图6 （b）］基本上再现了观测中年平均风速的分布型， 尤其是柴达木盆地与周边地区风速的差异得到了很好的体现， 但仍存在对高原西北部风速高值区的低估（0.5~2 m·s-1）及其对其他区域风速高估的误差（南部高估达2 m·s-1）。区域模式集合平均与观测的年平均风速空间相关系数值为0.44， 通过P<0.05的显著性检验。观测中高原地区春、 冬季平均风速较大， 夏、 秋季较小， 集合结果基本再现了这种季节差异， 但仍然存在与年平均类似的模拟偏差。高原西北部夏季平均风速存在明显的低估， 而整个高原的冬季平均风速存在显著的高估。春、 夏、 秋、 冬季四个季节， 区域模式集合平均与观测的空间相关系数值分别为0.41、 0.25、 0.45和0.46， 均通过P<0.05的显著性检验。动力降尺度模拟结果对高原地表风速模拟偏高的误差与以往全球模式结果较为一致， 可能与全球模式对青藏高原复杂下垫面以及物理过程的描述存在不足（陈说等， 2021； Wu and Shi， 2021）， 导致风速的模拟产生偏差， 这种偏差会引入区域模式， 从而对区域模式结果造成一定影响（江滢等， 2018； Wu and Gao， 2020）。此外， 高原西北部风速模拟偏低的误差可能与该地区观测站点稀少， 仅为133个（图1 ）， 插值所用的台站多源于周边风速较低的站点， 导致格点观测数据本身的风速偏小（吴佳和高学杰， 2013； Wu et al， 2017）。
图6 CN05.1（左列）和区域模式集合ENS（中列）的1986 -2005年青藏高原年平均和季节平均风速分布及其差值（右列）（单位： m·s-1）Fig.6 Spatial distribution of the annual and seasonal mean surface wind speed from CN05.1 （left column）， ensemble mean of downscaling simulations （middle column） and their differences （right column） over Qinghai-Xizang Plateau from 1986 to 2005.Unit： m·s-1 |
随后， 绘制了不同区域模式模拟结果及其集合平均的1986 -2005年青藏高原区域平均风速泰勒图（图7 ）。泰勒图分析能够给出5个模拟的青藏高原地区季节和年平均地表风速与CN05.1的空间相关系数、 标准化标准差（横纵轴）和标准化均方根误差（实线）， 结果离CN05.1（图中REF）越近代表模拟结果越接近观测（Taylor， 2001）。从图7 中可以看出， 5个降尺度模拟结果均能够模拟出CN05.1中青藏高原地表风速的基本分布型， 相关系数值普遍在0.3至0.5之间（均通过P<0.05的显著性检验）， 但较气温和降水变量明显较差（Wu et al， 2020）， 相比中国区域平均的结果也要差（Wu et al， 2021）。其中夏季的相关系数最小， 不到0.3， 标准化标准差在0.75左右，而冬季尽管相关系数在0.4~0.6之间， 但标准化标准差较大， 在1~1.5之间， 表明模式对夏季平均风速的空间分布、 强度以及冬季风的强度等模拟存在明显偏差。在5个模拟中， HdR对青藏高原地表风速的空间分布模拟相对较好， 相关系数较其他模拟结果要高， NdR对空间分布模拟较差。而MdR的标准化标准差在春季、夏季、 冬季及年平均均比其他模拟要小， 表明其对高原地表风速强度的模拟较好。值得注意的是， 不同模拟结果之间通常彼此接近， 表明5个降尺度模拟对再现青藏高原季节及年平均风速分布方面表现一致。
为进一步评估各个模拟结果对月、 季节和年平均风速等的模拟效果， 图8 给出了模式在青藏高原地区月、 季节和年平均地表风速的相对均方根误差。从图8 中可以看到， 单个动力降尺度结果中， HdR在1 -3月的表现最好， RMSE’的值小于0.95， 其次是MdR， RMSE’的值基本在0.95~1.05， 而CdR在1 -3月、 12月以及冬春季表现不佳。与大多数单个模拟结果相比， 集合平均的模拟性能更加稳定， RMSE’的值在0.95~1.05， 结果也一定程度上证实了Tebaldi and Knutti（2007）的发现， 即多峰平均值或中值可以部分抵消误差， 从而产生更好的性能。
进一步评估动力降尺度模拟结果对青藏高原地表风速年内变化的模拟效果， 结果如图9 所示。从图9 中可以看到， CN05.1中青藏高原地表风速最大值出现在春季， 这与以往基于台站的研究结果一致（姚慧茹和李栋梁， 2016）， 但5个模拟结果以及集合平均都显示其峰值出现在冬季， 主要原因是模式在冬半年（1 -3月及11 -12月）的高估较大， 超过1 m·s-1， 在夏半年（4 -9月）所有模拟值均接近于CN05.1的结果， 偏差较小。此外， 不同模拟结果之间的范围（最大最小值区间）即模拟的不确定性由图中柱体给出， 可以看到模拟结果之间的差异在冬半年较大， 下半年较小， 意味着冬半年不确定性较大， 而下半年的不确定性相对较小。
图9 区域模式模拟及其集合平均对1986 -2005年青藏高原平均风速年内变化点代表CN05.1， 柱体代表不同模拟结果的范围， 柱体中的线条代表集合平均 Fig.9 Boxplots of monthly mean surface wind speed during 1961 -2005 over Qinghai-Xizang Plateau from five downscaling simulations and their ensemble mean.Dots indicate the monthly mean wind speeds from the CN05.1， and the boxes indicate the simulation spread among five simulations， while the lines in the boxes indicate the ensemble mean |
最后还评估了各个模拟及其集合平均结果对青藏高原地表风速变化趋势的模拟效果。从表5 可以看到， 1986 -2005年CN05.1的年平均及季节平均地表风速均显示为显著下降的趋势（通过P<0.05的显著性检验）， 变化趋势值在-0.35~-0.22 m·s-1·（10a）-1， 春季下降最显著， 冬季次之， 秋季最小， 这与以往基于站点的结果基本一致（You et al， 2010； 姚慧茹和李栋梁， 2016）。5个动力降尺度模拟结果与CN05.1变化趋势值之间存在明显偏差， 变化趋势总体偏低。其中， 模式对秋季风速变化趋势的模拟最好， CdR、 HdR、 MdR以及ENS的模拟效果均较好， 与CN05.1的结果最为接近， 模式对冬季风速变化趋势也有较好的模拟能力。相比之下， MdR对青藏高原地表风速变化趋势比其他模式更接近观测， 其次是CdR和HdR， 其中MdR的变化符号与CN05.1均完全一致， 并且秋季和冬季风速的变化趋势值与CN05.1也最为接近。另外， 注意到， 尽管动力降尺度对青藏高原地表风速变化趋势的模拟还存在一定偏差， 但总体比以往的全球模式模拟结果要好， 尤其在秋、 冬季（江滢等， 2018； Wu et al， 2020）。
表5 CN05.1及区域模式模拟的青藏高原1986 -2005年区域平均风速变化趋势Table 5 Trends of the regional mean surface wind speed over Qinghai-Xizang Plateau from CN05.1， five simulations and their ensemble mean during 1986 -2005 |
| 模式 | 平均风速变化趋势/［m·s-1·（10a） -1］ | ||||
|---|---|---|---|---|---|
| 年平均 | 春季 | 夏季 | 秋季 | 冬季 | |
| CN05.1 | -0.27* | -0.35* | -0.26* | -0.22* | -0.29* |
| CdR | -0.05 | -0.08 | - 0.03 | -0.20 | 0.03 |
| EdR | -0.04 | 0.10 | - 0.03 | 0.02 | -0.10 |
| HdR | - 0.11 * | - 0.10 | 0.01 | -0.40* | -0.07 |
| MdR | -0.04 | -0.02 | -0.01 | - 0.22 | - 0.20 |
| NdR | -0.04 | 0.10 | - 0.03 | 0.02 | -0.10 |
| ENS | -0.06 | 0.01 | -0.02 | -0.18* | -0.09 |
*表示变化趋势通过P<0.05的显著性检验， 粗体表示与CN05.1变化符号一致， 斜体表示与CN05.1变化数值差异最小 |
4 讨论和结论
4.1 讨论
青藏高原地形复杂， 观测台站主要位于东部及南部， 北部至昆仑山北麓等地势高的地区则基本没有观测站点分布（吴佳和高学杰， 2013）， 这也决定了直接使用较少的观测台站来代表青藏高原平均风速会存在一定程度的低估。而本研究结果显示CN05.1与在青藏高原地区适用性较好的NCEP-1、 JRA-55再分析数据较一致， 并且与以往基于台站的研究结果也吻合（You et al， 2010； 姚慧茹和李栋梁， 2016）， 对台站数据偏低的误差有一定的改善。表明格点化观测资料在高原地区适用性较高。
尽管由于分辨率的提高， 区域模式相比以往的全球模式结果有一定的改进（江滢等， 2018； Wu et al， 2020）， 但其对高原地区地表风速的模拟能力还有较大的提高空间。首先， 全球模式对青藏高原复杂下垫面以及物理过程的描述存在的问题（陈说等， 2021； Wu and Shi， 2021）也会引入区域模式， 从而对区域模式结果造成一定影响（江滢等， 2018； Wu and Gao， 2020）。其次， 区域模式本身物理过程的限制也是一个方面（Gao et al， 2017； Wu et al， 2021）。另外， 高原西北部风速模拟偏低的误差可能与该地区位势太高， 观测站点稀少， 插值所用的台站多源于周边风速较低的站点， 导致格点观测数据本身也具有一定的不确定性（吴佳和高学杰， 2013； 韩海涛等， 2021）。未来可以通过提高区域气候模式分辨率、 完善内部物理过程等改进模式对高原风速的模拟效果（Zou and Zhou， 2013； Wu et al， 2021； Wu and Shi， 2021）， 也可以针对高原地区地表风速开展模式数据订正工作， 这对青藏高原未来气候变化预估、 影响评估与风险应对等研究具有重要的意义（陈说等， 2021； Wu et al， 2021）。
另外需要说明的是， 区域模式对青藏高原地表风速变化趋势的评估存在一定的问题。风速变化趋势同时受外强迫和气候系统内部变率影响。本研究采用了全球海气耦合模式驱动区域气候模式开展动力降尺度模拟， 该模拟无法再现观测气候系统内部变率的位相变化特征。（Chen et al， 2013； IPCC， 2021； 左志燕等， 2021）， 因此仅外强迫引起的那部分变化是可比较的。其模拟的变化趋势更接近观测也不足以说明该模式模拟能力更好， 可能仅仅说明其全球模式驱动场的低频变率位相偶然地与观测比较一致。从评估结果看， 不同降尺度模拟的风速均表现为下降趋势， 可能说明外强迫的影响占较大比重。另外， 模拟的趋势值比观测普遍偏小， 可能与全球模式地表边界条件没有考虑粗糙度的变化有关， 研究表明粗糙度增大对地表风速减少的贡献占了25%~60%（Vautard et al， 2010； Chen et al， 2012）。评估气候系统内部变率和外强迫对风速变化趋势的定量贡献需要借助数值模拟、 检测归因等方法实现（Vautard et al， 2010； 左志燕等， 2021）。
4.2 结论
基于格点化观测资料CN05.1的1961 -2020年青藏高原地表风速， 以及台站数据、 NCEP-1、 JRA-55和ERA5再分析风速资料， 分析了近60年来青藏高原月、 季节和年平均地表风速以及不同等级风速概率的变化及趋势。随后， 使用区域气候模式动力降尺度的5个模拟结果及其集合， 检验评估了区域模式对1986 -2005年地表风速空间分布、 逐月、 季节平均、 年内变化和长期变化趋势的模拟能力， 主要结论如下：
（1） 多种资料均显示1961 -2002年青藏高原地区的季节和年平均地表风速为减小趋势， 但数值上存在差异， CN05.1和NCEP-1的风速在2002 -2020年为上升趋势。对于不同等级风速概率的变化， 无风变化较小， 在1961 -2002 年和2002 -2020年均以微弱的减少趋势为主， 软风和轻风在1961 -2002 年增加， 在2002 -2020年减少。微风（3.4~5.4 m·s-1）、 和风（5.5~7.9 m·s-1）则基本表现出与平均风速类似的变化特征。大风（≥8 m·s-1）在1961 -2002 年和2002 -2020年均以减少趋势为主。
（2） CN05.1对青藏高原年平均风速和四个季节风速变化速率、 年际变率、 年内变化、 空间分布、 年际变率等方面均与站点数据最为吻合， 也表现出与NCEP-1和JRA-55再分析数据相一致的特征， 并且与JRA-55更加接近， 表明该格点化观测风速数据在高原地区的适用性较好。此外， 由于地形订正插值的考虑， CN05.1风速比站点普遍偏大， 导致两者不同等级风速概率之间的显著差异。CN05.1在1961 -2002年年平均及四季风速变化的趋势值与JRA-55资料最为接近， 在2002 -2020年与NCEP再分析更为接近； 在风速年际变率方面， CN05.1的表现与JRA-55和NCEP-1一致， 且与JRA-55更为接近； 在年内变化上CN05.1和JRA-55更加吻合； 在不同等级风速变化趋势上， CN05.1与NCEP-1正负变化符号更为一致， 但变化数值上与JRA-55更加接近； 此外， CN05.1与NCEP-1的软风和轻风的年际变率更一致， 与JRA-55的和风和大风的年际变率更加接近。
（3） 区域气候模式动力降尺度可以再现青藏高原平均风速的基本分布， 并且可以很好地模拟出青藏高原复杂地形区的风速差异， 但仍存在对高原西北部风速模拟偏小以及对其他区域风速模拟偏大的问题； 模式对夏季平均风速的空间分布以及冬季风强度等模拟存在明显偏差； 模式无法很好地模拟出高原地表风速的年内变化， 但对在夏半年（4 -9月）的风速变化模拟较好； 所有模拟结果均低估了青藏高原地表风速变化趋势； MdR和HdR对各个评估指标的模拟均要好， 集合平均的模拟性能则较单个模拟更为稳定。