太阳辐射是地球-大气系统最重要的能量来源, 它从根本上决定着地球-大气的热状况^[1]。同时, 它也直接影响着地球上的生物生长和人类活动, 是陆地生态系统生产力和碳收支的重要环境因子^[2]。对湖泊生态系统而言, 太阳辐射也是一个至关重要的影响因子, 它为湖泊藻类的大量生长和水华爆发提供了良好的物理条件^[3], 辐射大小变化将直接影响湖泊生态系统的结构、功能和演变。

太阳辐射作为模型研究的必要参数被广泛利用。一般的作物生长模拟模型、生态系统过程模型、气候模型和水生态动力学模型均需要逐日温度、降水、太阳辐射和日照时数等数据^[4-8], 但目前只有少数站点进行了太阳辐射的逐日监测。研究区域松花湖地区未设太阳辐射站点, 故缺少太阳辐射数据。随着气候、作物和水生态动力学模型的逐步深入, 太阳辐射资料的缺乏成为松花湖地区相关模型研究的限制因素。为弥补松花湖地区日总辐射数据及计算方法的空缺, 对松花湖地区的日总辐射计算模型进行了研究, 探讨了松花湖地区日总辐射与日照时数的相互关系及时空变化特征。

目前, 获取日总辐射的途径有直接观测和模型模拟计算2种, 但由于缺乏观测站, 大部分地区缺少太阳辐射的观测资料, 只能通过模型模拟计算日总辐射。国内外最常用的模拟计算方法是利用常规气象数据来计算日总辐射, 主要包括基于云量和日照时数计算2种方法^[9-12]。已有较多研究采用日照时数计算日总辐射, 如张运林等^[13]对无锡地区太阳辐射进行了计算和分析, 赵军等^[14]利用黄土高原区11个气象站点1992—2000年的日照时数和日总辐射数据, 采用最小二乘法计算得出该区域月总辐射模型参数。从这些研究可以发现, 使用日照时数进行太阳辐射估算是一种方便快捷、具有较高精度且易于推广的方法。但这些研究的计算和数据处理过程均较复杂, 主要是利用日照时数月平均值计算月平均日总辐射, 并未得到逐日日照时数与日总辐射的关系。为此, 运用广义模式搜索算法对长春和延吉地区1961—2000年逐日日照时数和日总辐射数据进行模型参数的率定和优化, 然后利用两地2001—2012年的日照时数和日总辐射数据进行验证，得到两地计算日总辐射的经验参数a_s和b_s。两地a_s和b_s平均值可作为松花湖地区日总辐射模型的经验参数, 为日后松花湖地区需要输入逐日太阳辐射的相关模型研究提供参考。

1 研究区域概况

松花湖即丰满水库, 位于吉林市东南部, 是在松花江小丰满处建水坝而形成的大型人工湖泊。流域面积43 370.8 km², 回水全长180 km, 湖面最宽约5 km, 最大水面面积550.0 km², 平均水深30~40 m, 年入库径流量约137亿m³。松花湖是长春市和吉林市居民生活和生产用水的水源地, 此外还有发电、防洪、灌溉、养殖和航运等功能。由于水土流失严重以及工业和旅游业发展等因素, 湖区生态环境已有一定程度的破坏和污染, 有些区域已呈轻度富营养化^[15]。

松花湖周边的太阳辐射站点为哈尔滨、长春和延吉, 但与长春和延吉相比, 哈尔滨距离松花湖较远, 所以选取长春和延吉2个站点的逐日太阳辐射与日照时数数据开展研究。长春和延吉2个站点的地理位置分别为43°54′ N, 125°13′ E和42°52′ N, 129°28′ E。松花湖(43°07′~43°50′N, 126°45′~127°38′ E)与2个站点的距离相近, 因此其模型参数可用2个站点模型参数的平均值来估算, 数据来源于中国气象科学数据共享网(http://data.cma.cn/)。

2 研究方法 2.1 基于日照时数的日总辐射模型

日总辐射R_s无法直接测得, 可用下式计算^[16]:

$ {R_{\rm{s}}} = ({a_{\rm{s}}} + {b_{\rm{s}}} \times n/N) \times {R_{\rm{a}}}。$

(1)

式(1) 中, R_s为日总辐射, W·m^-2; R_a为晴天日总辐射, W·m^-2; a_s和b_s均为经验参数; n为日照时数, h;N为日照时长, h;n/N为相对日照时间。

$ \begin{array}{l} {R_{\rm{a}}} = \frac{{24 \times 60}}{{\rm{\pi }}}{G_{{\rm{sc}}}}{d_{\rm{r}}}({w_{\rm{s}}} \times {\rm{sin}}\varphi \times {\rm{sin}}\delta + {\rm{cos}}\varphi \times \\ {\rm{cos}}\delta \times {\rm{sin}}{w_{\rm{s}}}), \end{array} $

(2)

$ {w_{\rm{s}}} = {\rm{arccos}}( - {\rm{tan}}\varphi \times {\rm{tan}}\delta ), $

(3)

$ \delta = 0.409{\rm{sin}}(2{\rm{\pi }}J/365 - 1.39), $

(4)

$ {d_{\rm{r}}} = 1 + 0.033{\rm{cos}}(2{\rm{\pi }}J365)。$

(5)

式(2)~(5) 中, G_sc为太阳常数, 取值1 360.8±0.000 5 W·m^-2; d_r为日地间相对距离的倒数; w_s为太阳时角, rad; φ为地理纬度; δ为太阳磁偏角, rad^[17];J为日序数。

2.2 松花湖日总辐射模型参数的估算

利用长春和延吉2个站点1961—2000年实测日总辐射和日照时数数据, 借助广义模式搜索(GPS)算法^[18], 在MATLAB中进行参数a_s和b_s的率定^[19]。GPS算法是直接搜索算法中的一个特殊子集族, 主要是在特殊方向集上抽取目标函数, 通过比较函数值大小, 找出下降方向进而解决最优化问题^[20]。

GPS算法的搜索步s_k由模式矩阵P_k(确定s_k的方向)和步长控制参数Δk(确定s_k的步长)共同决定。定义模式矩阵为

$ {\mathit{\boldsymbol{P}}_k} = \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{C}}_k} = [\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{M}}_k}-\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{M}}_k}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{L}}_k}\left] = \right[\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{L}}_k}], $

(6)

$ {\mathit{\boldsymbol{C}}_k} = [{\mathit{\boldsymbol{M}}_k}-{\mathit{\boldsymbol{M}}_k}{\mathit{\boldsymbol{L}}_k}\left] = \right[{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}{\mathit{\boldsymbol{L}}_k}]。$

(7)

式(6)~(7) 中, B为基本矩阵, 非奇异, B∈R^n×n, 在迭代过程中一般保持不变; C_k为生成矩阵, C_k∈Z^n×p(p > 2n); M为有限的非奇异矩阵集合, M_k∈M∈Z^n×n; L_k∈Z^n×(p-2n), 至少包含一个0列, 使搜索失败时可以尽快回到初始点, Γ_k∈Γ(Γ为有限的矩阵集合)，为模式的可行性方向, 保证模式有足量的方向进行搜索。由于B和C_k的秩均为n, 故P_k构成Rⁿ的正基。所以, 在任意迭代步k, 定义以x_k为中心的试探点x_kⁱ=x_k+S_kⁱ=x_k+ΔkBc_kⁱ, 其中搜索步长Δk∈R+, 搜索方向为Bc_kⁱ(C_k=[c_kⁱ, c_kⁱ⁺¹, …, c_k^p])。

对GPS算法的探索移动有如下2个假设:s_k∈ΔkP_k≡ΔkBC_k≡Δk[BΓ_kBL_k]; 若min{f(x_k+y), ∀y∈ΔkBΓ_k} < k < f(x_k), 则有f(x_k+s_k) < f(x_k)。

给定初始点x₀∈Rⁿ, 初始步长Δ0 > 0, 则计算步骤如下:计算f(x_k), 其中k=0, 1, 2, …, n; 按一定模式确定搜索步s_k; 计算ρ_k≡f(x_k)-f(x_k+s_k); 判断若ρ_k > 0成立, 则x_k+1≡x_k+s_k, 否则x_k+1≡x_k; 更新生成矩阵C_k和步长Δk。当步长控制参数Δk足够小时, 认为上述算法已经收敛, 算法终止。

GPS算法不仅构思直观, 易于编程实现, 且可以有效收敛到全局最优解, 解决线性约束条件下的非线性复杂水文模型的参数优化问题。该算法既无需计算或近似任何导数信息, 又不用计算或近似任何惩罚因子或拉格朗日乘子, 不强加任何充分下降的概念, 仍然可以实现算法的收敛性, 在非线性规划和非光滑最优化领域中有着广泛应用。

选用Nash-Sutcliffe效率系数(NSE, E_NS)和均方根误差(RMSE, E_RMS)进行模型结果准确性检验^[21-22], 计算公式为

$ {E_{{\rm{NS}}}} = 1 - \sum\limits_{i = 1}^M {{{(R_{{\rm{s, }}i}^{{\rm{sim}}} - R_{{\rm{s, }}i}^{{\rm{obs}}})}^2}/\sum\limits_{i = 1}^M {{{(R_{{\rm{s, }}i}^{{\rm{obs}}} - \overline {R_{\rm{s}}^{{\rm{obs}}}} )}^2}} }, $

(8)

$ {E_{{\rm{RMS}}}} = \sqrt {\sum\limits_{i = 1}^M {{{(R_{{\rm{s, }}i}^{{\rm{sim}}} - \overline {R_{\rm{s}}^{{\rm{obs}}}} )}^2}} /M} 。$

(9)

其(8)~(9) 中, R_s^obs和R_s^sim分别为实测和模拟值; $\overline {R_{\rm{s}}^{{\rm{obs}}}} $为实测平均值; M为样本数。NSE反映了模拟值与实测值的拟合程度, NSE值≤1。NSE值=1说明完全拟合; NSE值≤0说明拟合度极差; NSE值越接近1, 说明拟合效果越好, 模型模拟精度更高。RMSE同样能反映模型模拟结果, 其值越小, 说明模拟效果越好。

从长春和延吉日总辐射模型参数a_s和b_s的率定结果中选取最佳参数, 采用2001—2012年日照和辐射数据进行验证, 从而确定长春和延吉日总辐射模型参数a_s和b_s, 松花湖日总辐射模型参数用2个站点模型参数的平均值来估算。

3 结果与分析 3.1 长春和延吉日总辐射模型参数的率定与验证

图 1~2分别为1961—2000年长春和延吉逐日日总辐射的实测和模拟值散点图。由图 1~2可知, 散点较均匀的分布于1:1直线的两侧, 模型率定期间的NSE值均为0.88, RMSE值分别为29.349和28.115 W·m^-2。由此说明, 长春和延吉两地日总辐射模型在率定期模拟所得的日总辐射与实测日总辐射拟合程度较好。

为检验长春和延吉两地日总辐射模型的精度, 采用2001—2012年逐日日总辐射和日照时数实测值进行模型验证。由图 3~4可知, 两地2001—2012年日总辐射实测值散点图均匀地分布在1:1直线两侧, NSE值分别为0.93和0.94, RMSE值分别为22.39和19.53 W·m^-2。与率定期相比, 验证期的NSE值有所提高, 而RMSE值有所降低。NSE值在小于1的情况下越高, 表明模型精度越高, 而RMSE值越低则表明模型的验证结果更好。由此表明, 两地的日总辐射模型均能较好地模拟逐日日总辐射, 模拟结果可行, 模型精度良好。

3.2 模型参数率定结果

表 1为各地日总辐射模型参数率定结果。

其中松花湖模型参数为长春和延吉参数的算术平均值。由表 1得到松花湖地区日总辐射的计算公式为R_s=(0.164 34+0.591 02×n/N)×R_a。

4 松花湖日总辐射与日照时数特征分析

由于缺乏松花湖日总辐射与日照时数数据, 将长春和延吉2个站点的日总辐射和日照时数平均值作为松花湖估算值, 以此来粗略观察松花湖日总辐射与日照时数的多年和年内变化趋势。

从图 5可以看出, 近50 a来松花湖地区日总辐射呈波动减小趋势, 这与全国日总辐射的变化趋势基本保持一致^[23-24]。1961—1968年日总辐射均高于多年平均值, 1969—1981年日总辐射除2个年份外均低于多年平均值, 1992—2010年日总辐射大体上接近多年平均值, 波动较小。松花湖地区日照时数与日总辐射呈正比例关系, 即日照时数增加则日总辐射增加。

图 6为1961—2012年松花湖多年月平均日总辐射与日照时数变化。由图 6可知, 松花湖月平均日总辐射曲线呈上凸的抛物线形式, 夏季最大, 春秋季次之, 冬季最小, 除4、6月外, 其变化趋势与月平均日照时数基本相同。对月平均日总辐射与日照时数进行相关性分析, 发现两者在α＝0.01水平呈极显著相关, 相关系数为0.845。1—5月(除4月外)松花湖地区多年月平均日总辐射与日照时数均呈上升趋势, 其中日总辐射在6月达最高, 月日照时数在5月达最高。分析其可能原因是6月为松花湖丰水期^[25], 降水较多, 所以日照时数较5月少; 但6月太阳对北半球进行直射, 导致其日总辐射增大, 所以6月日总辐射达最高。9月之后, 月平均日总辐射与日照时数均下降, 在12月均达全年最低值。周晓宇等^[26]对1961—2009年东北地区日照时数变化特征的分析发现, 东北地区最高和最低日照时数分别出现在5和12月, 笔者的研究结果与其一致。

影响地区太阳辐射强弱变化的因素有很多, 除各地区地理位置差异外, 还包括大气分子、云量、大气中水汽及悬浮物等。松花湖太阳辐射与日照时数之间的关系十分密切, 日照时数对太阳辐射值的影响较大。松花湖地区附近除吉林外没有其他大型城市, 受城市热岛效应的影响小, 大气中悬浮物相应也较少, 所以大气悬浮物对其影响可能较小。水汽对太阳辐射具有吸收作用, 相关研究表明, 1970—1990年东北地区水汽通量增加, 但在2000年代出现大幅度减小^[27]。由图 6可知, 1970—1990年松花湖日总辐射总体呈减少趋势, 在2000年代则相对稳定, 说明水汽对松花湖地区日总辐射有一定影响, 但不是主要影响因子。云量是另一个影响因素, 云量增加, 太阳总辐射减少。这是因为云层对太阳辐射具有吸收和散射作用^[28], 削弱了到达地面的太阳辐射。但松花湖地区云量数据缺乏, 无法准确了解云量对太阳辐射的影响。谢今范等^[29]研究表明, 云量对东北地区太阳辐射由减少到增加的趋势毫无贡献, 而低云量及气溶胶对其影响较大, 推测低云量及日照时数变化对松花湖日总辐射有较大影响。

5 讨论

目前, 基于日照时数计算太阳辐射的研究很多。ALMOROX等^[30]在没有考虑大气透明度系数的基础上, 认为逐日太阳辐射模拟计算的经验参数a_s和b_s之和为0.75左右, 并通过对西班牙的太阳辐射模拟得到a_s和b_s分别为0.217 0和0.545 3;联合国粮食及农业组织(FAO)建议a_s和b_s分别取值0.25和0.5;张运林等^[13]利用上海、南京和杭州近40 a的实测总辐射与月日照百分率进行回归分析, 得到无锡地区的a_s和b_s值分别为0.135 1和0.570 7;而左大康等^[31]根据我国不同类型地区实测日总辐射、日照百分率月平均值和晴天状态下的月总辐射资料, 计算得到a_s和b_s分别为0.248和0.752。笔者计算得到松花湖地区a_s和b_s分别为0.164 34和0.591 02, 两者之和接近0.75, 与ALMOROX等^[30]的系数之和保持一致。笔者计算的参数值与部分研究存在差异, 分析可能的原因有2个:一是前人主要基于日照百分率的月平均值进行分析和计算, 而笔者是利用逐日太阳辐射与日照时数率定优化而得, 时间尺度不一样, 导致其值可能存在一定差异; 二是研究区域范围不同, 笔者计算得到的是一个相对较小区域的模型参数值。

以往对太阳辐射的研究大多是对全国范围内太阳辐射进行估算。童成立等^[32]通过对我国不同区域9个代表站的模拟分析构建模型, 该模型仅需要输入站点的日照时数和地理信息, 便可以估算出该站点的太阳辐射。赵东等^[33]利用全国范围内的日照时数对直接辐射进行估算, 从而分析我国近50 a来太阳直接辐射资源的基本特征及其变化。以上研究并没有对特定区域进行太阳辐射模型的参数计算和公式反演, 而对于特定环境区域, 模拟参数需要进行调整。该研究通过对长春和延吉逐日太阳辐射模型的模拟和验证, 得到太阳辐射估算模型, 以后只需输入该区域内的日照时数数据便可得到日总辐射值, 一方面可以弥补松花湖区域内日总辐射模型研究的空缺, 另一方面也为日后松花湖区域作物及水生态动力学模型研究提供参考。

由于缺乏松花湖地区实测日总辐射数据, 未对松花湖日总辐射模型公式进行验证。但长春和延吉两地距离松花湖地区较近, 且纬度与松花湖地区相差不大, 因笔者提出的模型公式具有一定的可信度, 对以后松花湖地区日总辐射计算具有一定的参考意义。目前松花湖地区缺少太阳辐射站点, 导致需要这些气象资料的水生态模型研究难度较大, 希望以后能完善松花湖地区的气象站点, 获取更准确的气象资料, 这很大程度上有利于松花湖的管理与利用。

6 结论

(1) 利用长春和延吉2个站点的逐日日总辐射与日照时数数据进行模型参数的率定, 得出2个站点的经验参数a_s和b_s, 松花湖日总辐射模型参数为2个站点参数的算术平均值, 由此得出松花湖地区日总辐射计算公式:R_s=(0.164 34+0.591 02×n/N)×R_a。

(2) 松花湖日总辐射与日照时数呈正比关系, 近50 a来松花湖地区日总辐射总体呈减少趋势, 且1982—2010年间日总辐射变化较小, 日照时数变化趋势与之一致。松花湖多年月平均日总辐射最高和最低值分别出现在6和12月, 多年月平均日照时数则在5月达最高, 在12月达最低。