基于云模型的西北地区辐射时空分布研究

关 静 梁 川 赵 璐 崔宁博 王春懿 姜守政

(1.四川大学水利水电学院， 成都 610065； 2.四川大学水力学与山区河流开发保护国家重点实验室， 成都 610065)

0 引言

地表总辐射量(Rs)，即到达地球表面的太阳辐射量，是地球-大气系统能量平衡的关键因素[1-3]，在作物生长模型、蒸散量估算、灌溉制度制定和气候变化研究等领域都有重要的科学价值和现实意义[4-7]。然而，由于建设和维护成本很高，辐射量观测受到很大限制，其数据并不像常规日照时数、气温等气象数据那样容易获得。我国西北地区地域辽阔，地形复杂，以高原、盆地和山地为主，由于下垫面形式多样，使得各地气候要素分布很不均匀，辐射条件存在较大差异[8-9]。

目前，西北地区共有171个气象站点，太阳辐射观测站点仅35个，尚在运行且连续观测30 a以上的站点仅13个，虽然经过非常严格的质量控制，但有观测条件的Rs数据在时间序列上并不完整，数据缺失现象仍普遍存在。依靠这些稀疏且分布不均的站点的观测资料无法描述该地区Rs时空变化特性。为解决Rs观测资料不足的问题，学者们建立了多种估算Rs的替代方法，如经验模型[5-7]、卫星图像[10]、随机天气模型[11]、机器学习[12]等方法。其中，基于日照时数的Ångström-Prescott(A-P)模型[5-6]是国内外应用最广泛且精度较高的方法。

为提高A-P模型在预测Rs方面的精度，许多学者对其参数取值、率定方法进行了诸多研究。其中，关于A-P模型是否应该分时段率定参数引起了国内外学者的广泛讨论。SOLER[13]发现在欧洲气候环境下使用A-P模型每月参数比固定一个参数估算Rs精度更高；PODESTA等[14]研究阿根廷潘帕斯高原辐射与日照关系后，建议使用月参数或者季度参数来避免时间规律造成的系统误差；HUSSAIN等[15]研究发现A-P模型参数分夏半年(4—9月)和冬半年(9—3月)率定相比多年一次率定可提高Rs预测精度；LI等[16]发现，我国长江流域分月率定要略好于多年一次率定；然而LIU等[17]通过对我国温带季风气候下的20个站点研究发现，分时段率定A-P模型参数相比多年一次率定并不能提高Rs预测精度。云模型是李德毅等[18]在传统模糊集理论和概率统计的基础上提出的一种定性定量不确定性的转换模型，目前已应用于降水量、参考作物蒸散量和干旱等的时空分布描述[19-21]，而关于Rs的时空分布描述还未提及。本文拟在前人研究基础上，基于西北地区16个辐射站点4个分区1995—2015年实测Rs与日照时数(n)率定A-P模型参数a、b，使用A-P模型4种不同参数率定方法(M1：分月率定，M2：分季率定，M3：分半年率定，M4：多年率定)，选取其中精度高且简便的方法计算4个分区1961—2015年共55 a的Rs，利用云模型来描述西北地区Rs时空分布，以期为西北地区各分区A-P模型参数率定方法提供理论基础，提高Rs预测精度，构建完整的Rs时间序列。

1 研究区域和资料

1.1 研究区域

中国西北地区位于73°25′～110°55′E 和31°35′～ 49°15′N 之间，地处青藏高原北部和东北部，包括新疆、青海、甘肃、宁夏、陕西(陕南除外)和内蒙古西部，总土地面积约为中国国土面积的1/3[22]，是中国主要的干旱和半干旱区。该区地处大陆腹地，远离海洋，复杂的地形导致气候要素分布很不均匀，是生态环境极其脆弱的地区[8,22]。

系统聚类分析[23-24]是一种根据多种地学要素对地理实体进行划分类别的方法，其基本思想为距离最近的样本先聚成类，距离远的样本后聚成类，该过程一直进行下去，每个样本最终能聚到适当的类中。降水与气温是气候变化的主导因素，不同地区降水与气温的状况不同，气候特点也就不同[25]。故本文根据西北地区具有30 a以上观测资料的156个气象站点的年平均降水量(P)和年平均气温(T)(其空间插值图如图1a所示)，并采用欧氏距离和最远邻距离的方法对156个站点进行系统聚类分析，同时结合地形地貌，将西北地区划分为4个区域，得到西北地区分区图，如图1b所示。由图1b可以看出：Ⅰ区主要包括新疆北部、甘肃河西走廊中西段、宁夏中北部、内蒙古西部；Ⅱ区主要包括新疆南部；Ⅲ区主要为青海省；Ⅳ区主要包括甘肃东南部、宁夏南部、陕西关中地区与陕北地区。

图1 西北地区气象站点分布Fig.1 Meteorological stations distribution in Northwest China

1.2 基础资料

结合资料系列的完整性和站点分布的均匀性，选取西北地区16个代表站点，其1995—2015年日辐射数据序列数据缺失率为0.01%～1.24%，平均缺失率0.32%。站点基本情况如图1b和表1所示。气象数据来自中国气象数据共享服务网(http:∥data.cma.cn/site/index.html)，包括1995—2015年逐日Rs，1961—2015年逐日n。尽管数据发布前已经经过严格的质量控制，但仍存在部分数据缺失的情况，因此在使用前进行了如下处理：①在率定A-P参数时，如果1995—2015年实测Rs和n有1个缺失时，则删除该日所有数据。②如果Rs/Ra(Ra为大气层顶部所接收的太阳辐射量)或n/N(N为最大可能日照时数)大于1，则删除该日全部数据。③利用1961—2015年n计算Rs时，对少量缺测数据n采用多年平均值法补全。

表1 西北地区16个气象站点概况Tab.1 Summary of 16 meteorological stations in Northwest China

2 研究方法

2.1 Ångström-Prescott(A-P)模型

基于n的Rs估算模型是由Ångström首先提出，并由Prescott改进而形成的经验模型，被称为Ångström-Prescott(A-P)模型。其表达式为[5-6]

(1)

其中

(2)

(3)

(4)

ωs=arccos(-tanφtanδ)

(5)

(6)

式中Gsc——太阳常数，取0.082 MJ/(m2·min)

dr——日地相对距离

J——日序数

ωs——太阳高度角，rad

φ——纬度，rad

δ——赤纬，rad

a、b——介于0～1之间的经验系数，二者之和为晴空条件下辐射传输率

A-P模型的参数校正通常使用最小二乘法，建立Rs/Ra与n/N的关系。利用西北地区16个辐射站点1995—2015年实测Rs与n率定参数a、b，使用4种校正方法，将长系列以日为单位的Rs/Ra与n/N多年相同月份(1—12月)进行最小二乘回归，共率定12次，记为M1；将长系列以日为单位的Rs/Ra与n/N分季节(春季3—5月，夏季6—8月，秋季9—11月，冬季12—2月)进行最小二乘回归，共率定4次，记为M2；将长系列以日为单位的Rs/Ra与n/N分半年(夏半年4—9月和冬半年9—3月)进行最小二乘回归，共率定2次，记为M3；将长系列以日为单位的Rs/Ra与n/N进行最小二乘回归，产生一组参数，记为M4。选择适合西北地区不同区域的校正方法，将率定后的参数a、b代入A-P模型，利用1961—2015年实测n估算Rs，构建完整的Rs时间序列。

2.2 统计分析

使用辐射模型常用4个评价指标：平均偏差MBE(MBE)、均方根误差RMSE(RMSE)、归一化均方根误差nRMSE(nRMSE)和t统计量，其计算式为

(7)

(8)

(9)

(10)

式中Si——模拟值Oi——实测值

MBE、RMSE和nRMSE越小，说明模拟效果越好。其中，nRMSE小于等于10%，认为模拟效果较好； nRMSE为(10%,20%]，认为模拟效果良好；nRMSE为(20%,30%]，认为模拟效果合理；nRMSE为(30%,40%]，认为模拟效果一般；nRMSE大于40%，认为模拟效果较差[6]。

2.3 云模型

2.3.1云模型的数字特征

云模型的数字特征用期望Ex、熵En和超熵He3个数值来表征[19-21]。

期望Ex是云滴在论域空间分布的期望。熵En是定性概念随机性的度量，既反映了云滴的离散程度和不均匀程度，又反映了云滴的取值范围，对于论域的定性概念有贡献的云滴，主要落在区间[Ex-3En，Ex+3En]。超熵He是熵的不确定性度量，即熵的熵，由熵的随机性和模糊性共同决定，反映了不确定的凝聚性，超熵He越大，隶属度的随机性越大，云的“厚度”也越大，越不稳定。云模型的3个数字特征值把模糊性和随机性完全集成到一起，构成定性和定量相互之间的映射，为定性与定量相结合提供了有利手段。

2.3.2正向云发生器

正向云发生器是从定性到定量的映射，其算法的主要过程为输入3个特征值，即Ex、En和He，以及云滴数M；输出M个云滴xi及其确定度ui(i=1,2,…,k)。具体计算步骤如下:

(11)

(12)

(3)计算确定度

(13)

(4)生成具有确定度ui的xi成为数域中的一个云滴。

(5)重复步骤(1)～(4)k次，产生要求的k个云滴为止。

2.3.3逆向云发生器

应用云模型理论分析西北地区Rs在时空分布上的不均匀性，根据云模型的数学特征，分别利用云滴Xi(Rs)在16个站点55 a的定量数值来还原出云的3个特征参数。

逆向云发生器是实现从定量值到定性概念的转换模型，其算法的过程主要为；输入样本点Xi(i=1, 2,…,k)；输出数字特征(Ex，En，He)。Ex、En和He计算式分别为

(14)

(15)

(16)

(17)

(18)

S2——样本方差

3 结果与分析

3.1 A-P模型不同率定方法参数比较

利用西北地区16个辐射站点1995—2015年实测Rs与n，使用A-P模型4种校正方法(M1：分月率定，M2：分季率定，M3：分半年率定，M4：多年率定)率定A-P模型参数a、b。为明确分时段率定A-P模型参数a、b的变化情况，本文以M1为例，所得西北地区各站点A-P模型参数a、b变化情况如表2所示。由表2可以看出：各个站点用M1率定的参数a、b变化显著。其中，乌鲁木齐站分月率定的参数a、b和a+b变幅最大，其a值介于0.10～0.29，变异系数为38%，b值介于0.40～0.66，变异系数为14%，a+b值介于0.66～0.80，变异系数为5%；和田站分月率定的参数a、b和a+b值变幅最小，其a值介于0.24～0.27，变异系数为4%，b值介于0.46～0.51，变异系数为4%，a+b值介于0.71～0.75，变异系数为1%。不同站点参数a、b各月平均值变化也十分显著。其中，乌拉特中旗站参数a的平均值(0.26)最大，伊宁站和延安站参数a平均值(0.13)最小；伊宁站参数b平均值(0.62)最大，阿克苏站参数b平均值(0.48)最小。总的来说，西北地区分月率定的A-P模型参数a值介于0.09～0.39，b值介于0.40～0.66，a+b值介于0.63～0.85，变化幅度a>b>(a+b)。

3.2 A-P模型不同率定方法对Rs的影响

为评价4种校正方法在西北地区的Rs模拟精度，将日尺度上的Rs模拟值和实测值进行线性拟合，结果如表3所示。由表3可以看出：Rs模拟值与实测值拟合结果较好，均呈极显著相关(P<0.01)，Ⅰ区拟合效果最优。其中，Ⅰ区4种校正方法的决定系数(R2)介于0.974～0.978，变异系数为0.17%，Ⅱ区4种校正方法的R2范围为0.960～0.961，变异系数为0.05%，Ⅲ区4种校正方法的R2范围为0.944～0.946，变异系数为0.09%，Ⅳ区4种校正方法的R2范围为0.941～0.942，变异系数为0.05%，由变异系数可以看出4种校正方法在各个分区的R2无明显差异。

表2 M1校正方法下西北地区A-P模型参数Tab.2 A-P model parameters calibrated by M1 in Northwest China

表3 西北地区各分区不同校正方法Rs模拟值与实测值之间的决定系数Tab.3 Determination coefficients between estimated Rs by different calibrated methods and measured Rs in different aeras of Northwest China

注：** 表示极显著相关(P<0.01)。

通过将多年平均月尺度上各方法在西北地区的Rs模拟值与实测值进行比较，计算各方法月尺度Rs模拟值与实测值的平均相对误差，结果如表4所示。由表4可以看出：各方法月尺度Rs模拟值与实测值的平均相对误差由小到大基本表现为M1、M2、M3、M4，相对误差均在较小范围内，4种校正方法拟合精度均较高。

图2 4种校正方法计算Rs日值的统计分析箱线图Fig.2 Statistical analysis box charts of four calibrated methods

校正方法Ⅰ区Ⅱ区Ⅲ区Ⅳ区M10.1150.3360.4220.353M21.0441.2300.7490.743M31.6781.3000.8811.519M42.4061.3960.9541.496

4种校正方法计算Rs日值的统计分析箱线图如图2所示。图2显示Ⅰ区4种校正方法的MBE中位数均为0，离散程度由小到大为M1、M2、M3、M4，且RMSE、nRMSE分布相似，除M4外，其余3种方法的t检验值(排除异常值)均小于1.65(P<0.05)(注：若t检验值小于1.65(P<0.05)，说明该方法率定下Rs模拟值和实测值无明显差异)。Ⅱ区M1、M2和M3的MBE中位数均为0，M4的MBE值均小于0，4种校正方法RMSE、nRMSE分布相似，M1、M3的t检验值均小于1.65(P<0.05)，M2、M4的t检验值部分大于1.65。Ⅲ区M1、M2的MBE中位数均为0，且分布集中，M3的MBE均小于0，M4的MBE中位数大于0，分布较离散，4种校正方法RMSE、nRMSE分布相似，M1、M2、M3的t检验值均小于1.65(P<0.05)，M4的t检验值部分大于1.65，且离散程度较高。Ⅳ区M1、M2的MBE中位数均为0，且分布集中，M3、M4的MBE均小于0，4种校正方法RMSE、nRMSE分布相似，M1、M2的t检验值均小于1.65(P<0.05)，M3、M4的t检验值部分大于1.65，且离散程度较高。

综上所述，日尺度和月尺度上4种方法Rs模拟值与实测值拟合结果均较好。各分区4种校正方法计算Rs的RMSE、nRMSE分布相似，MBE略有差异，t检验结果表明Ⅰ区和Ⅲ区使用M1、M2、M3计算的Rs值与实测值无明显差异，即可认为M1、M2、M3在Ⅰ区和Ⅲ区的预测精度等效；Ⅱ区使用M1、M3计算的Rs值与实测值无明显差异，即可认为M1、M3在Ⅱ区的预测精度等效；Ⅳ区使用M1、M2计算的Rs值与实测值无明显差异，即可认为M1、M2在Ⅳ区的预测精度等效。选取其中精度高且简便的方法，即Ⅰ区、Ⅱ区和Ⅲ区各站点用M3率定A-P模型参数a、b，Ⅳ区各站点用M2率定A-P模型参数a、b，各站点具体使用的A-P模型参数见表5、6。

表5 Ⅰ区、Ⅱ区、Ⅲ区各辐射站点具体使用的A-P模型参数Tab.5 A-P model parameters of radiation sites fromⅠ area, Ⅱ area and Ⅲ area

表6 Ⅳ区各辐射站点具体使用的A-P模型参数Tab.6 A-P model parameters of radiation sites from Ⅳ area

3.3 Rs时间分布特性的云模型分析

Ⅰ区、Ⅱ区、Ⅲ区各站点用M3率定A-P模型参数a、b，Ⅳ区各站点用M2率定A-P模型参数a、b，将率定后的参数a、b代入A-P模型，利用1961—2015年实测n估算Rs，构建完整的西北地区Rs时间序列。平均西北地区4个分区各站点春、夏、秋、冬及多年平均日Rs，得到各区域不同时间尺度的Rs值。以西北地区4个分区的1961—2015年共55 a的Rs作为研究对象，根据逆向云发生器的算法计算各自云模型的数字特征(表7)，然后根据正向云发生器的算法计算云滴并分别绘制西北地区不同时间尺度的隶属云图(图3)。云模型中期望Ex是Rs的平均值，反映了该区1961—2015年春、夏、秋、冬及多年的平均日Rs值。熵En体现了Rs相对于平均值的确定度，也反映Rs相对于平均值的离散度，其值越大，Rs在55 a分布越不均匀。超熵He反映了熵的离散程度，其值反映不均匀性的稳定性。

由表7中的4个分区云模型的数字特征可以看出，春季Rs值由小到大为Ⅱ区、Ⅳ区、Ⅰ区、Ⅲ区，夏季Rs值由小到大为Ⅳ区、Ⅲ区、Ⅱ区、Ⅰ区，秋季Rs值由小到大为Ⅳ区、Ⅰ区、Ⅱ区、Ⅲ区，冬季Rs值由小到大为Ⅰ区、Ⅱ区、Ⅳ区、Ⅲ区。春、夏、秋、冬均表现为Ⅰ区、Ⅱ区、Ⅲ区Rs时间分布的不均匀性较小且较稳定，Ⅳ区Rs时间分布的不均匀性较大且不稳定。年尺度上，Rs强度由小到大为Ⅳ区、Ⅱ区、Ⅰ区、Ⅲ区，Ⅰ区Rs时间分布的不均匀性较小但不稳定，Ⅱ区、Ⅲ区Rs时间分布的不均匀性较小且较稳定，Ⅳ区Rs时间分布的不均匀性较大且不稳定。整个西北地区春、夏、秋、冬及年尺度上Rs平均值分别为18.799、21.515、13.376、9.540、15.808 MJ/(m2·d)。其中，春、夏、秋和年尺度均表现为Rs在时间分布上的不均匀性较小且较稳定，冬季Rs在时间分布上的不均匀性较小但不稳定。这些特征也可由数字特征值经逆向云发生器算法得到的隶属云图(图3)中直观看出。

表7 西北地区及各个分区不同时间尺度Rs时间分布云模型的数字特征Tab.7 Cloud model digital characteristics of time distribution of each area at different time scales of Rs in Northwest China MJ/(m2·d)

图3 西北地区不同尺度Rs时间分布隶属云图Fig.3 Time distribution clouds of Rs at different scales in Northwest China

3.4 Rs空间分布特性的云模型分析

根据16个代表气象站点在55 a春、夏、秋、冬及年尺度的多年平均Rs，分析Rs在空间上的变化特性，西北地区Rs的空间分布如图4所示。

图4 西北地区不同尺度Rs空间分布图Fig.4 Space distributions of Rs at different scales in Northwest China

由图4可以看出，西北地区Rs空间分布不均匀且具有季节性差异。夏季Rs值最大，其次是春季，秋季和冬季Rs值较小。4个季节均表现为Ⅲ区(青海省)Rs值较大，这与青海省特殊的地形地貌和气候区域条件有关。青海省位于青藏高原东北部，是“世界屋脊”青藏高原的一部分，海拔1 650～6 860 m，光能资源极其丰富[26]，故其常年日照辐射强烈。同时通过计算1961—2016年春、夏、秋、冬及年尺度Rs在空间分布上的云模型数字特征可知，春、夏、秋、冬及年尺度空间分布上的熵En分别为1.628、1.456、1.634、1.809、1.352 MJ/(m2·d)，远大于时间尺度上的0.597、0.457、0.333、0.300、0.236 MJ/(m2·d)，这说明各站点的多年Rs相对于其多年平均值较为离散，同时春、夏、秋、冬及年尺度空间分布上的超熵He分别为0.423、0.688、0.300、0.475、0.326 MJ/(m2·d)，远大于时间分布上的0.035、0.015、0.012、0.107、0.016 MJ/(m2·d)。因此，与Rs在时间上的分布特性相比，其在空间上的分布特性更不均匀、且不均匀性更不稳定。其中，夏季Rs空间分布的不均匀性最不稳定。这也可以通过图5和图3的比较直观得出。

图5 西北地区不同尺度Rs空间分布隶属云图Fig.5 Space distribution clouds of Rs at different scales in Northwest China

4 结论

(1)各站点分月率定的参数a、b变化显著；日尺度上4种方法Rs模拟值与实测值线性拟合结果均较好；多年平均月尺度上各方法计算Rs的相对误差由小到大基本表现为M1、M2、M3、M4，相对误差均在较小范围内。

(2)各分区4种校正方法计算Rs的RMSE、nRMSE分布相似，MBE略有差异，t检验结果表明Ⅰ区使用M1、M2、M3计算的Rs值与实测值无明显差异，Ⅱ区使用M1、M3计算的Rs值与实测值无明显差异，Ⅲ区使用M1、M2、M3计算的Rs值与实测值无明显差异，Ⅳ区使用M1、M2计算的Rs值与实测值无明显差异。

(3)Ⅰ区Rs时间分布的不均匀性较小但不稳定，Ⅱ区、Ⅲ区Rs时间分布的不均匀性较小且较稳定，Ⅳ区Rs时间分布的不均匀性较大且不稳定。整个西北地区春、夏、秋和年尺度均表现为Rs在时间分布上的不均匀性较小且较稳定，冬季Rs在时间分布上的不均匀性较小但不稳定。

(4)西北地区Rs空间分布不均匀，4个季节均表现为Ⅲ区(青海省)Rs较强；与Rs在时间上的分布特性相比，其在空间上的分布特性更不均匀、不稳定。