变差函数理论在湿季降水量推算全年降水量中的应用

赵楠楠 杨振京 杨庆华

摘要 基于石家庄市16个气象站点1981—2010年逐日降水资料，根据变差函数理论分析并获取了时间变差函数参数，利用该参数作为克里金插值的约束信息推算全年降水量。结果表明，湿季（5—9月）推算结果明显优于其他时段，且具有较强的稳定性;该时段的推算误差率普遍小于11%，误差率的平均值为-1.24%，可作为推估全年降水量的有效方法。

关键词 时间变差函数;降水量;克里金插值;降水特征;湿季

中图分类号 P 331  文献标识码 A

文章编号 0517-6611（2019）12-0225-06

doi：10.3969/j.issn.0517-6611.2019.12.062

开放科学（资源服务）标识码（OSID）：

Abstract Based on the daily rainfall data at 16 meteorological stations in Shijiazhuang from 1981 to 2010，the semivariogram was analyzed to calculate the variogram parameters （range;sill and nugget），and then the annual precipitation was estimated by Ordinary Kriging（OK） interpolation constrained by these parameters.The results showed that the estimated precipitation of wet season was closer to the actual precipitation compared to others periods and more stability.The calculated error rates of test data in Shijiazhuang were normally less than 11% and the mean error rate was 1.24%.It indicated that time semivariogram was an useful method to estimate annual precipitation under the condition of wet seasonal precipitation was known.

Key words Time semivariogram;Precipitation;Kriging interpolation;Precipitation characteristics;Wet season

連续可靠的降水数据是水文分析、水资源管理和气候研究的重要前提。但是，受人为、经济和气候等因素影响难以获取连续、完整的降水数据，因此通过已有数据推算补全缺失资料就显得十分必要。气候要素本身就是一个随机变量[1]，统计学方法是气候学家们用来研究气候要素常用的方法。利用统计学方法研究降水的文献很多，主要是针对其时空分布[2-4]、频数[5]、降水周期[6]、强度特征[5]和变化趋势[7-8]、预估[9]等。降水资料缺失主要分为2个方面：空间资料缺失，如王刚等[10]通过与该站临近的、相关性高的站点补齐缺失数据;时序资料缺失，张强等[11]用三次样条函数内插方式补齐缺失数据;晏利斌[7]利用多元回归插值方法补齐缺失数据;张志萍等[12]则通过雨量站分区间利用相关站点汛期降水量占全年比率插补全年降水量。通常，气象站点有完整的湿季降水资料而缺少全年完整数据，如张志萍等[12]曾指出大理河流域许多雨量站只记录了5—10月份的降水资料，显然其不能代表年降水量，所以需要通过湿季降水数据扩推全年降水量。然而，对于年内缺失湿季之外降水推算全年降水量的研究很少。变差函数是研究变量空间变化极为有利的工具，是地质统计学的基本内容之一，利用地质统计学研究降水[13-14]、水位[15]等水文参数空间变异特征取得了较好的研究效果。因此，笔者选取石家庄16个气象站点的累年逐日降水资料，分析、获取各站点时间变差函数参数，以这些参数为插值约束进而通过湿季降水量推算全年降水量。

1 资料与方法

1.1 数据来源

中国气象局气象数据中心开放全国国家级站点地面气象资料，覆盖了绝大多数的县级区域。该研究整理和分析了经中国气象局检验修正后（包括16个县站点1981—2010年累年逐日降水量、石家庄站点1955—2015年逐月降水量）的地面气象资料，这些资料具有可靠性高、时间连续性好等特点，得到广泛认可[10]。

1.2时间变差函数 变差函数是用来研究区域变量的空间变化特征的有力工具，根据变差函数的研究理论，把变差函数中的位置和距离替换为时间和时间差，则可以用来分析变量的时间相关性，称其为时间变差函数[16]，时间变差函数γ（t，f）可表达为：

γ（t，f）=12E[Z（t）-Z（t+f）]2（1）

式中，t代表时间，f代表随机过程中的时间间隔，Z（t）为t时刻随机变量的取值，E[]为随机变量的期望。如果假设随机过程是二阶平稳的，时间变差函数表达为：

γ（f）=12E[Z（t）-Z（t+f）]2（2）

根据变差函数理论，当f=0时，γ（f）的取值最小，两者之间的相关性最强;在f趋于变程（α）的过程中γ（f）值逐渐增大，两者之间大的相关性随之减小;当f>α时，γ（f）值不再改变，两者之间没有相关性。

由于受实际取样有限等因素的影响，需先制作试验变差函数，试验变差函数计算公式为：

γ*（f）=12N（f）N（f）i=1E[Z（ti）-Z（ti+f）]2（3）

式中，N（f）是f间隔下样本点对的数目。试验变差函数确定之后，通过一种合适的理论变差函数对其进行拟合，从而形成理论变差函数模型。常用的变差函数理论模型有球状模型、高斯模型、指数模型、线性模型等，由于球状模型可以较好地拟合各站点试验变差函数计算结果，故该研究使用的理论模型为球状模型，其表达式为：

γ（f）=C0f=0

C0+C（32fα-12f3α3）0

C0+Cf>α（4）

式中，C0为块金效应，C为基台值，C0+C为总基台值。

1.3 时间变差函数方法的步骤与原则

①对研究区周边地区站点的累年逐日降水资料进行整理与分析，获得研究区相关站点的时间变差函数各项参数。如果研究区有气象站点，可直接使用该站点通过分析后的变差函数参数;如果没有站点，则可以通过周边站点插值后获取。

②获取累年逐日最大降水量（Pday）以及累年逐月最大降水量（Pmonth）。分析研究区湿季各月降水量，超过Pmonth降水量的月份把超出降水量分配给邻月，优先分配给邻月降水量大的月份;月间降水量分配完毕之后，进行月内降水量分配，月内每天的降水量以Pday优先分配，不能以整数天分配Pday的，Pday分配完后，最后一天分配剩余降水量，把这些分配降水的天放置在该月中间。③利用第一步获得的时间变差函数各项参数作为约束，对湿季的日降水量进行普通克里金插值，最后提取该年每天降水量，求和后即可得到该年的全年降水量。

2 石家庄地区年降水特征

石家庄气象站点在20世纪50年代中期就有连续且完整的年降水资料，其中5—9月多年平均降水量占全年降水量的80%以上，因此该研究中把该时段定义为研究区的湿季。选取石家庄站点的1955—2015年的逐月降水资料，计算每年湿季降水量占该年降水总量的比率（图1）。从图1可以看出，近61年石家庄地区年降水量变化幅度较大，其中1996年降水量最大，高达1 097.1 mm，1972年降水量最小，仅226.0 mm，两者数量悬殊相差甚大。从比率曲线可以看出，石家庄湿季降水量（Pwet）占该年降水量的比率变化范围也较大，在1968年比率最小，仅55.04%，在1988年比率最大，达97.72%，几乎全年降水量都在湿季完成，同时也反映了年内降水量的分配也是极不均匀的。此外，降水量最多（最少）的年份，其占全年降水比率并不是最大（最小），且每年比率各不相同，因此常规的方法难以推算该区年降水量。

3 时间变差函数方法在石家庄地区的应用

3.1 石家庄市气象站点累年逐日降水资料相关性分析

通过对收集到的石家庄16个县站点累年逐月、累年逐日降水数据统计与分析，发现累年逐月、累年逐日数据通过球状模型都可以很好地拟合原始数据。但是，经进一步研究（以行唐为例），发现累年逐月降水数据规律性较差（图2），虽然根据原始数据可以较好地拟合变差函数，但其可能是由于数据点太少而造成的“假象”，如果将其应用到后期的插值过程中，可能会产生较大的误差，从而导致错误的结果，因此该研究仅对累年逐日降水数据的时间变差函数进行分析与利用。

石家庄16个县站点累年逐日降水数据都符合正态分布规律，满足高斯域變差函数基本要求。各站点试验时间变差函数的拟合结果见图3（一维随机变量），通过球模型可以很好地拟合原始数据（R2均在0.94以上），拟合出的各项参数见表1。从图3和表1可以看出，各站点之间的块金效应不同，变化范围较大，为0.54～1.40，相差近2倍;基台值为5.124～7.976，变程为136.4～144.2，它们主要受控于累年日降水量多少及其时间分布。平山站点的块金效应最大（1.40），但是其变程不是最大的;同样，行唐的块金效应最小（0.54），其变程并不是最小的。上述情况可以用变差函数性质解释：变程越大，随机变量的非均质性越小，但块金效应增大会在一定程度上增加随机变量的非均质性。

3.2 石家庄站点各项时间变差函数参数的确定

时间变差函数的各项参数与日降水量及其时间分布有密切的联系，而降水量的空间分布通常可以由插值方法确定，在气象数据稀少等特殊区域，时间变差函数的各项参数也可以用降水量空间插值方法确定。克里金插值、距离反比插值是降水空间插值中最常用的方法。许多学者研究发现，克里金插值比距离反比插值精度更高，结果更可靠[17]。但是在对石家庄16个气象站点的时间变差函数的块金效应、基台值和变程进行变差函数拟合时，发现这些参数空间变异性规律极差（图4），如果使用克里金插值可能会产生较大的误差。另外，庄立伟等[18]在研究东北地区逐日降水空间插值方法时发现，距离反比插值精度比克里金插值精度高，但平滑程度较小，更适合日降水量的空间插值，因此石家庄站点的各项参数由距离反比方法确定，其站点各项参数经插值后结果见表1。

3.3 石家庄站点Pday和Pmonth的确定和分配

石家庄站点时间变差函数各项参数获取后，需要确定该站点的Pday和Pmonth。该站点距离藁城和正定2个站点较近，其大致位于2站的中间位置，因此石家庄站点的Pday和Pmonth由2处的平均值确定（表1）。Pday和Pmonth确定完以后开始进行月间与月内分配。下面以石家庄2016年的Pwet为例（表2）进行分配，具体步骤如下：7月降水量445.6 mm远大于Pmonth（127.0 mm），故7月降水量分配127.0 mm，剩余部分优先分配邻月降水量大的月份，依次分配给5、6、8、9月。分配完成后的各月降水量见表2。

月间降水量分配完毕后，进行月内降水量分配，其中，石家庄站点Pday为9.45 mm。以5月为例优先分配9.45 mm，其中88.4/9.45≈9.35，则连续分配9天Pday，第10天分配剩余降水量3.35 mm，即5月的11—19日日降水量均为9.45 mm，20日为3.35 mm。Pday分配完毕后，把分到降水量的10 d放在该月中间位置，其余日期降水量分配为0mm，另外4个月份做法相同。湿季之外的月份降水量设为未知，通过后期的插值获得，到此，降水量分配完毕。

3.4 年降水量计算结果及分析

收集到了石家庄站点2011—2016年和2015年正定、井陉等6个站点的逐月降水数据进行验证上述方法的可行性。首先，把各个站点按照上述的步骤处理完毕之后导入到ArcGIS 10.2.2软件中，然后，在表1中对应站点时间变差函数参数约束下进行普通克里金插值，处理范围为（0～366，-1～1），最后，提取该年内相应天数的插值数值求和后即可得到该年的年降水值。测验数据的原始相关信息见表3，为证明湿季推算结果可靠性，提供了相邻不同时段的计算结果（表4）。

从表3可以看出，同一站点不同年份以及同一年份不同站点之间Pwet、Pa和月降水量差别明显，这符合石家庄地区比率、年降水量变化幅度较大以及年内降水分布不均的降水特征。

误差率是推算精度的重要指标，从图5可以看出，3—7、6—10和7—11月的计算降水量明显的偏离Pa，它们的推算结果普遍难以接受（表4），不能获得理想的推算效果。虽然4—8月推算结果在特殊情况下产生比5—9月份较优的推算结果，但是4—8月份的推算结果具有很大的不稳定性，最大误差达37.68%，因此，如果采用4—8月降水资料推算Pa，其推算结果可能很大程度偏离Pa，无法获取具有参考价值的推算结果。5—9月份推算结果则相对稳定，虽然可能产生高达16.91%的误差，但是其推算结果一直游离在Pa附近，具有较高的参考价值。

利用5—9月降水资料推算实际降水量（Pa）在石家庄地区获得了较好的应用，其计算误差率普遍在11%以下，其中石家庄2013年计算误差率仅为2.11%，说明该方法可以作为一种全年降水量的推算方法;但石家庄2011年和无极2015年推算误差却分别达16.91%和-12.35%（表4）。误差较大的原因可能是由于Pmonth和Pday不能代表站点自身降水特征以及克里金的圆滑效应引起的。

由于无法收集充足的历史数据确定去掉奇异值后的Pmonth和Pday，对无极和平山借用相邻站点的Pmonth和Pday进行重新计算（表4），这2个站点的推算误差得到了有效的缩小，故在推算前需要对Pmonth和Pday的代表性进行检验，尽可能地体现气象站点自身的降水特征。

克里金插值的圆滑效应也是误差来源之一。在占全年比率大致不变的情况下，Pwet过高或过小都会导致误差率较大;在Pwet大致相近的情况下，占全年比率的过高或过小同样也会导致误差率较大（图6）。前者在比率较大的情况下，由于克里金法的圆滑效应，降水量多或少会导致其周围未分配降水量日期的插值结果多或少，并对误差进行传递，从而造

成较大的误差率。在Pwet大致相同的情况下也是这种原因。

如在2011年石家庄站点，Pwet占全年比率高达89.05%，降水量为600.4 mm，克里金插值的圆滑效应使得该年产生较大的误差率;而2013年石家庄站点比率高达89.26%，但是其Pwet仅为453.7 mm，虽然克里金法仍有圆滑效应存在，但是其引起的绝对误差较小，因而其误差率很小。2015年藁城站点則相反，其Pwet及其比率都较小，从而产生较大的误差。

4 结论

利用石家庄16个气象站点的逐日降水资料，借助变差函数相关理论，通过对湿季降水量推算全年降水量方法的研究，得出以下结论：

（1）不同时段的推算结果有很大差异，其中湿季的推算结果最为可靠。3—7、6—10和7—11月推算结果较大程度地偏离实际降水量（Pa）;4—8月推算结果有较大的不稳定性;而湿季推算结果更可靠，其推算结果和误差的稳定性均较好。

（2）时间变差函数方法在石家庄测验数据中的推算误差普遍小于11%，其中石家庄2011年推算误差最小，仅为2.11%;其误差率的平均值为-1.24%，方差为8.66，具有较高的参考价值，可以作为一种推估全年降水量的有效方法。

（3）时间变差函数方法在石家庄2011年和无极2015年测试数据中误差较大是由于Pmonth和Pday不能代表站点自身降水特征及克里金插值的圆滑作用引起的。其中无极和平山站点借助各自临近站点的Pmonth和Pday使得误差得到一定程度的缩小。

参考文献

[1]江志红，丁裕国，陈威霖.21世纪中国极端降水事件预估[J].气候变化研究进展，2007，3（4）：202-207.

[2] 张秀娟，陈晓光，王尧，等.西北四省区降水的时空变化特征分析[J].安徽农业科学，2012，40（18）：9809-9812.

[3] CUI L F，WANG L C，LAI Z P，et al.Innovative trend analysis of annual and seasonal air temperature and rainfall in the Yangtze River Basin;China during 1960-2015[J].Journal of atmospheric and solarterrestrial physics，2017，164：48-59.

[4] 徐利岗，杜历，姚海娇，等.中亚干旱区降水时空变化特征及趋势分析[J].干旱区资源与环境，2015，29（11）：121-127.

[5] 王志福，钱永甫.中国极端降水事件的频数和强度特征[J].水科学进展，2009，20（1）：1-9.

[6] 刘健，夏军，王明森，等.1961—2015年山东省降水周期变化特征[J].人民黄河，2017，39（4）：6-10.

[7] 晏利斌.1961—2014 年黄土高原气温和降水变化趋势[J].地球环境学报，2015，6（5）：276-282.

[8] 吴慧，吴胜安.近48年海南省极端降水时空变化趋势[J].安徽农业科学，2010，38（19）：10101-10103.

[9] ALIZADEH M J，KAVIANPOUR M R，KISI O，et al.A new approach for simulating and forecasting the rainfallrunoff process within the next two months[J].Journal of hydrology，2017，548：588-597.

[10] 王刚，严登华，张冬冬，等.海河流域1961年-2010年极端气温与降水变化趋势分析[J].南水北调与水利科技，2014，12（1）：1-6，11.

[11] 张强，孙鹏，陈喜，等.1956～2000 年中国地表水资源状况：变化特征、成因及影响[J].地理科学，2011，31（12）：1430-1436.

[12] 张志萍，冉大川，慕志龙.大理河流域降水资料插补方法探讨[J].人民黄河，2006（12）：26-27，78.

[13] 李巍，范文义，毛学刚，等.降雨量空间插值方法比较研究[J].安徽農业科学，2014，42（12）：3667-3669.

[14] DIRKS K N，HAY J E，STOW C D，et al.Highresolution studies of rainfall on Norfolk Island Part II： Interpolation of rainfall data[J].Journal of hydrology，1998，208（3）：187-193.

[15] OHMER M，LIESCH T，GOEPPERT N，et al.On the optimal selection of interpolation methods for groundwater contouring： An example of propagation of uncertainty regarding interaquifer exchange[J].Advances in water resources，2017，109：121-132.

[16] 裴韬，周成虎，李全林，等.基于变差函数分析的地震时间相关性定量估算[J].地震，2002，22（2）：17-21.

[17] PIAZZA A D，CONTI F L，NOTO L V，et al.Comparative analysis of different techniques for spatial interpolation of rainfall data to create a serially complete monthly time series of precipitation for Sicily，Italy[J].International journal of applied earth observations & geoinformation，2011，13（3）：396-408.

[18] 庄立伟，王石立.东北地区逐日气象要素的空间插值方法应用研究[J].应用气象学报，2003，14（5）：605-615.