基于地质统计学理论的海拉尔河流域降水时空变异性研究

基金项目:水利部公益性行业可研专项经费(201201050)；国家自然基金青年基金(51109115).

作者简介:刘晓民（1981），男(蒙古族)，内蒙古赤峰人，博士，主要从事3S技术与水资源管理方面研究。Email：13204717007@163.com

通讯作者:刘海燕（1978），女，内蒙古赤峰人，高级工程师，主要从事水利信息化方面研究。Email：lhy_lj@sina.comDOI：10.13476/j.cnki.nsbdqk.2014.04.004

第12卷总第73期·南水北调与水利科技·2014年第4期刘晓民等·基于地质统计学理论的

海拉尔河流域降水时空变异性研究

2,3

（1.内蒙古农业大学 水利与土木工程建筑学院,呼和浩特 010018；

2.清华大学 水沙科学与水利水电工程国家重点实验室,北京 100084；3.北京信息职业技术学院 计算机工程系,北京 100018；）

摘要：将地质统计学理论引入到水文水资源系统领域，并结合随机函数理论，利用少数测站的年降水量资料，确定出海拉尔河流域年降水量的实验变差函数。然后，进行不同方向实验变差函数的筛选、归并、拟合及其套合，形成水文变量几何异向性的结构函数，应用多时段各向异性泛克立格空间估计理论，实现了任意空间点、不同大小流域平均年降水量的分析计算。最后，通过与算术平均法、距离反比法、泰森多边形法等方法计算结果的对比分析与论证，证实了该理论与方法的准确性。

关键词：海拉尔河流域；地质统计学；多时段；泛克立格

中图分类号：P333.1文献标志码：A文章编号：16721683（2014）04001605

Spatial and Temporal Variation of Rainfall in the Hailar River Basin Based on Geostatistics Theory

LIU Xiaomin1,WAN Zheng1,LIU Haiyan2,3

(1.College of Water Conservancy and Civil Engineering,Inner Mongolia Agricultural University,Inner Mongolia,010018,

2.State Key Lab oratory of Hydroscience and Engineering,Tsinghua University Beijing 100084,China;

3.Department of Computer Engineering,Beijing Information Technology College,Beijng,100018,China)

Abstract:The geostatistics theory was introduced into the research field of hydrology and water resources and combined with random function theory to determine the experimental variation function of annual precipitation in the Hailar River Basin based on the measured annual precipitation data at a few stations.Then,the structural function of geometric anisotropy of hydrological variables was obtained by selecting,merging,and fitting the experimental variation functions in different directions. The theory of spatial estimation with anisotropy and universal kriging of multiple time periods was applied to achieve the analysis and calculation of average annual precipitation at an arbitrary point in space and in a basin with various scales.Finally,the advancement,convenience,and accuracy of the theory and method used in the paper were confirmed through the comparison of the results determined by the method of arithmetic average,inverse distance method,and Thiessen Polygon method.

Key words:Hailar River Basin;geostatistics theory;multiple time;universal kriging

流域平均降水量的计算方法主要有泰森多边形法、算术平均法、距离反比法或距离平方反比法（美国称客观运行法）、等降水量线法，以及上世纪70年代发展起来的地质统计学方法，，也简称克立格法克立格法对于平稳和非平稳地质信息的处理，是较为有效而精确的方法［13］，它包括简单克立格法、泛克立格法、泛克立格法、协同克立格法4种方法。对于水文变量而言，大多数学者多采用泛克立格法进行计算，但泛克立格法方法存在缺陷：其一，由于我国水文测站在大江大河上分布较多，而在中小河流上则稀疏分布，为了更好地研究局部地区和流域的降水量或其他水文变量的空间变异性，必须通过周边大量水文信息插补研究区域的数据，工作量较大；其二，对于小流域特别是水文测站较少或无水文测站的小流域，该方法在利用邻域水文资料时，由于各变量成因类型不同，破坏了变量一致性假设。

为了克服克立格法在水文领域应用的局限性，本文在计算海拉尔河流域上各流域平均降水量时，将多时段泛克立格空间分析理论［46］进一步延伸，采用几何异向性的变差函数来模拟将多时段水文变量的结构函数，从而考虑了不同方向上水文数据的分布特点与变异程度，以解决小流域或局部区域内少站点情况下降水量的最优估计问题。

1多时段泛克立格空间估计理论

首先假定一时空随机函数Z（x，t），式中x代表空间坐标，t为时间维。对流域内的任一测站xi，Z（xi，t）为一随机过程，而Z（x，t）在各个时刻tj（j=1，2，…，T）的实现则组成一个时间序列Z（xi，tj）；多个测站、在各时刻，Z（x，t）的观测值则组成一个随机时间序列簇。在有限邻域内空间上任意两站的时间序列都是相关的，测站距离越大，其相关程度逐渐减弱。

1.1基本假设

① x，y随机过程Z（x，t）是一阶平稳的，即：

E［Z（x，t）］=mx（1）

②x、y，随机过程Z（x，t）、Z（y，t）的互方差函数C（x，y，t）是平稳的，即：

E［Z（x，t）Z（y，t）］=mxmy+C（x，y）（2）

③ t，根据随机过程Z（x，t）的性质，确定其增量［Z（x，t）-Z（y，t）］具有非平稳的方差函数和非平稳的数学期望［m（x，t）-m（y，t）］，且它们都具有一阶平稳的特性。对x、y、t：

E［Z（x，t）-Z（y，t）］=m（x，y）-m（y，t）=mx-my（3）

12E［Z（x，t）-Z（y，t）］2=r（r，y，t）=r（x，y）（4）

式中：r（x，y）为变差函数。

④设Z（x，t）可以分解为飘移mx和剩余R（x，t）两部分，即Z（x，t）=mx+R（x，t），其中，mx代表长历时现象；R（x，t）代表短周期内变化。

1.2变差函数的计算

在确定上述假设的前提下，结合各水文测站的时间序列{Zi（tj）；i=1，2，…，N；j=1，2，…，T}，推算出实验变差函数，其公式如下：

mi=1T∑Tj=1Zi(tj)（5）

Ri(tj)=Zi(tj)-mi（6）

空间点i(xi)处剩余Ri(t)的一个实现用Ri(tj)代替，则不同两点i、k在时段tj内的剩余Ri(tj)、Rk(tj)可形成一组数据对，如时段为tj(j=1,2,…，T)，则可获得T个数据对Ri(tj)、Rk(tj)，即可估计i、k两点剩余的变差函数：

r(i,k)=12T∑Tj=1［Ri(tj)-Rk(tj)］2（7）

如果结构函数r(i,k)在各向都是同性的，则r(i,k)与i、k两点的方向没有关系，与两点的距离存在一定关系，则上式变形为：

r(i,k)=r(hik)=T2T∑11T［Ri（tj）-Rk（tj）］2（8）

式中：r（h）代表空间上每两站的年降水量时间序列。如果有N个测站，即可获得［N（N-1）/2］个r（h）。以各水文测站之间的距离作h为横坐标，以其对应的r（h）为纵坐标，绘制实验变差函数图。

根据对不同水文测站所获取的资料分析可知，受地形、地貌、地质成因条件等影响，大多数水文变量都呈现出较为明显的异向性。根据流域或地区内水文测站的分布情况，可按照前文所列的方法，做出各向的变差函数图，再将各向变差函数图套合，最终形成套合变差函数。

1.3多时段泛克立格方程组

前文方法求得的变差函数在各个时段内是相同的，但可根据站年降水量序列{Zi(tj);i=1,2,…，N；J=1，2，…，T}，利用泛克立格空间估计理论，估计任意时段t0、空间任一点x0处的Z（x0，t0）值。其估计式为：

Z*（x0，t0）=∑Nd=1λt0dZd（t0）（9）

式中：λt0d为t0时段第个测站观测值Zd（t0）的估计权重。

假设区域化变量Z（x，t0）的飘移式为：

m（x，t0）=kl=0alfl（x）（10）

式中：fl（x）为第l个飘移基函数；al为其系数。

通常，飘移式多采用线性飘移式或二次飘移式。在二维空间下，也就是：

线性式：m（x1，x2，t0）=a0+a1x1+a2x2（11）

二次式：m（x1，x2，t0）=a0+a1x1+a2x2+a3x21+a4x22+a5x1x2（12）

式中：（x1，x2）为空间点x的坐标；x1、x2、x21、x22、x1x2均为基函数；a0为飘移常数因子。

可依据地区实际情况，初步确定飘移式形式。然后根据数据的空间变化特点或趋势面分析法进一步确定。确定的形式要通过模型检验，即给定多个飘移式形式，分别代入泛克立格方程组：

∑Nβ=1λt0βr（xd，xβ）+∑Kl=1ulfl（xd）+u0=r（xd，x0）；

D=1，2，…，N

∑Nd=1λt0d=1

∑Nd=1λt0dfl（xd）=fl（x0）l=1，2，…，K（13）

其估计方差为：σ2ukt0=∑Nd=1λt0dr（xd，x0）+∑Kl=1ulfl（x0）+u0（14）

式中：u0、u1为拉格朗日乘子；r（xd，xβ）为变量在空间点xd，xβ之间的变差函数值。

上述泛克立格估计方程式可用于：

①估计任一时段t0，第个降水量测站的Zi（t0）值。据此可进行模型检验和各测站缺测资料插补。

②估计任一时段t0，空间任一点x0处的Z（t0，t0）。

③估计任一时段t0，空间任一区域V（x0）处的Z（x，t）的平均值，即估计：

1V（x0）∫V（x0）Z（x，t0）dx（15）

2研究实例

本文选取海拉尔河流域作为分析范围，共收集流域内及其邻域共34个站点的1956年-2006年51年同步期降水量资料，论证本方法及其他常规计算方法的差别及精度。

2.1年降水量的平稳性检验

通过傅立叶周期分析，海拉尔河流域及其邻域内各站年降水量都具有7 a、18 a和44 a左右的周期，对海拉尔河流域及其邻域内的17个气象站、8个降水量站和9个水文站的1956年-2006年降水量资料统一进行44 a滑动平均，降水量站、水文站点结果见图1（由于篇幅限制气象站不予列出，

图1各降水量站、水文站年降水量的44 a滑动平均值变化

Fig.1Variation of sliding average value of annual precipitation

in 44 Years at different rainfall stations and hydrologic stations

下同）。可以看出，海拉尔河流域及其邻域内各站年降水量资料基本是平稳的，满足对系列要求的平稳性假设条件。

根据海拉尔河流域及其邻域内各站年降水量8个44 a大周期滑动平均值，分别计算了各站点年降水量滑动平均值变化百分率。发现除海拉尔区气象站、海拉尔水文站、免渡河站、红花尔基水文站外，其余各站均呈波动变化减少趋势；而海拉尔区气象站、海拉尔水文站变化平稳；免渡河站、红花尔基水文站呈增加趋势；全流域的年降水量变化率表现出从上游到下游逐渐增加的规律，符合该地区上游主要为高山和林地、下游为草原的地貌特征；河谷地区年降水量变化率较小，均在0~1.5%之间。

2.2空间飘移式的确定

根据各个测站1956年-2006年年降水量资料估计它们的多年平均降水量mi（i=1，2，…，40），以mi为因变量，各站所在位置的坐标（xi，yi）（为相对坐标）为自变量进行了趋势面分析，结果见表1。

通过趋势面分析发现，年降水量存在线性飘移关系和非线性飘移关系，即：

m（x1，x2）=a0+a1x1+a2x2

m（x1，x2）=a+a1x1+a2x2+a3x21+a4x22+a5x1（16）表1海拉尔河流域及其邻域内各站1956年-2006年年降水量的趋势面分析

Table 1Trend surface analysis of annual precipitation in the Hailar River Basin and its neighborhood stations from 1956 to 2006

测站类型站点相对x

坐标相对y

坐标相对坐标的二次项1980年-2006

年平均值线性趋势

模拟值非线性趋

势模拟值x/kmy/kmx^2/万km2y^2/万km2 xy/万km2P/mmP^/mmP*/mmPP^

/mmPP*

/mm降水量站水文站免渡河357.65 442.712.82 962.3194.7460.3417.9414.442.446黑山头257.55 572.46.63 105.2143.5325.4350.8342.5-25.5-17.2兴安402.15 410.216.22 927217.6443.8446.3441.7-2.52.1维纳河3205 371.410.22 885.2171.9424.5399.5406.925.117.6托河5065 563.225.63 094.9281.5510.9499.751111.2-0.1银阿451.75 613.520.43 151.1253.6486.6464.5490.122.1-3.5塔尔气364.45 318.313.32 828.4193.8454.9429443.22611.7拉布达林301.35 598.59.13 134.3168.7393.6375.6377.818.115.8乌尔其汗381.95 49414.63 018.4209.8451.5429.5426.42225.1牙克石330.45 469.410.92 991.4180.7408.1400.2395.97.912.2红花尔基276.15 354.47.62 867147.9410.5374.2389.236.321.3伊敏牧场264.25 395.772 911.4142.6332.6364.8369.9-32.2-37.3辉河口259.95 437.96.82 957.1141.3323.7359.8356.9-36.1-33.2海拉尔263.25 457.16.92 977.9143.6357.2360.7355.3-3.51.9头站258.65 490.76.73 014.7142283.2356.1347.4-72.9-64.1完工202.55 4604.12 981.1110.5305324.3314.6-19.3-9.6嵯岗144.85 472.72.12 995.179.2267.9289.1270.7-21.3-2.82.3稳健实验变差函数的获取与变差函数的最优

拟合针对海拉尔河流域，分别将34个测站1956年-2006年降水量资料代入公式（5）、（6）、（7）计算出了海拉尔河流域不同空间点之间年降水量的实验变差函数。为了获得较为稳健的实验变差函数，删除了特异邻差值，采用稳健统计估计中的切尾估值法，消除了个别误差较大的资料对实验变差函数的影响。在此基础上，分0°、90°、45°、-45°四个方向，并以方向搜索变化角度±22.5°为标准，筛选出不同方向上的实验变差函数点，组成不同方向上的实验变差函数，见图2。

由图2可以看出，不同方向上的实验变差函数点群较居中，趋势明显，可以初步认为海拉尔河流域年降水量的空间结构特征是异向性的。

根据四个方向年降水量实验变差函数变化趋势，选用球状理论变差函数模型（17），用VPLOT程序对其进行了最优拟合。结果见表2。

r（h）=0，

c0+c132（ha）-12（ha）3

c0+c1，

h=0

0

h≥a(km)（17）

由表2可以看出，不同方向上的理论变差函数其基台值相同，只是变程各不相同。通过绘制年降水量的方向变程图2，发现本区年降水量具有几何异向性。经优选模拟，求得方向变程图的半长轴为555 km、半短轴为315 km、方向角为58°、各向异性比为1.762。

图2（a）海拉尔河流域0°、90°方向年降水量的变差函数

Fig.2（a）Variogram of annual precipitation in 0°

and 90° directions in Hailar River Basin

图2(b)海拉尔河流域45°、-45°方向年降水量的变差函数

Fig.2(b)Variogram of annual precipitation in 45° and

-45° directions in Hailar River Basin

表2年降水量不同方向理论变差函数的拟合参数

Table 2Fitting parameters of theoretical variation function of

annual precipitation in different directions

方向参数C0 /mm2C/mm2a/km0°505 200367.27590°505 200457.65345°505 200511.732-45°505 200291.6852.4理论变差函数和飘移关系式的检验

套合理论变差函数模型和飘移关系式确定好后，要经过检验方可用于降水量的估计。检验过程如下：对于任一时段（年）tj和测站xi，利用本时段tj内其它各站的年降水量资料代入公式（10）求出测站i的年降水量估计值Z*i（tj）；求Zi（tj）-Z*i（tj）。如果对所有的测站i（i=1，2，…，N）及时段tj（j=1，2，…，T）满足：① 误差的平均值Me近似于零；② 误差的方差σ2e接近于平均克立格估计方差S*2；③ 误差的直方图看上去似一个正态分布；④ 误差的绝对标准偏差直方图看上似一个正态分布，且误差的绝对值在1.96倍的平均克立格标准差内的数据个数大于95%，即可确认该理论变差函数和飘移关系式。

图3和表3给出了海拉尔河流域内部分水文站历年降水量实测值与克立格法估计值的对比。可以看出：乌尔其汗水文站、红花尔基水文站、陈巴尔虎旗气象站，估计值比实测值偏小；而伊敏牧场水文站、辉河口水文站、头站水文站、嵯岗水文站、鄂温克族自治旗气象站，估计值比实测值偏大；牙克石水文站、海拉尔水文站、完工水文站、海拉尔区气象站、牙克石市气象站、兴安降水量站、维纳河降水量站、免渡河降水量站，估计值与实测值相差不大。而海拉尔河流域整体来看，实测值与估计值仅相差0.9 mm，充分证明克立格法的计算精度很高。

表3年降水量的检验统计结果

Table 3Statistics results of annual precipitation

mm

项目水文站乌尔

其汗牙克

石红花

尔基伊敏

牧场辉河

口海拉

尔头站完工嵯岗实测值429.5385.7382.8322.9314.7347.3286.3299.4275.5估计值409.5380.4349.5344.7328.2339.2336.9300.3286.9偏差205.333.321.813.48.150.60.911.4项目气象站降水量站海拉

尔区牙克

石市陈巴

尔虎

旗鄂温克

族自治

旗兴安维纳

河免渡

河海拉尔河

全流域实测值347.6381.5321.5325467.7404.4422.7357.2估计值340.4389.3299.4337.6459.3409.4417.2358.0偏差7.27.82212.68.455.50.9检验结果的几个主要统计参数：

Me=-0.4(mm)，σ2e=4231，S*2=4229，P{|Z-Z*|≥196S*}=66/1683=3.98%。

海拉尔河流域年降水量估计误差直方图和绝对标准化偏差直方图分别见图4和图5。从流域检验结果的统计参数可以看出，上述几个检验条件均已满足，说明理论变差函数和飘移关系式选择准确、拟合较优。本模型可用于海拉尔河流域年降水量时间和空间上的内插估值及不同大小流域平均年降水量的计算。

2.5流域平均降水量的检验分析

对每一个小流域而言，以2 km×2 km的梅花形网格剖分流域，以网格交点作为估值点，利用模型式（10）、式（11）对各点历年年降水量进行最优估计，则某年流域内所有估值点估计值的算术平均值即为该年流域平均降水量。

为了检验克立格法计算流域平均年降水量的精度与有效性，分别利用算术平均法和泰森多边形法计算了各流域1956年2006年流域平均年降水量：

图4年降水量实测值与克立格法估计值的对比

Fig.4Comparison on measuredannual precipitation and

estimated annual precipitation using universal kriging

图5年降水量估计误差直方

Fig.5Estimation error histogram of annual precipitation

图6年降水量绝对标准化偏差直方

Fig.6Absolute normalized deviation

histogram of annual precipitation

A.算术平均值法

针对各类型降水量站点，流域内部各站点取权重1，流域外部各站点取权重0.5，然后加权平均计算出11个流域平均降水量。

B.泰森多边形法

以流域内各类型降水量站点所控制的多边形面积与流域面积的比值为权重系数，然后加权平均计算出11个流域平均降水量。

将上述各水文站以上流域三种方法计算的流域平均年降水量进行对比，见图6。

a.对于乌尔其汗水文站以上流域，算术平均值法、泰森多边形法计算的流域平均降水量相差不大，但两种方法结果均高于克立格法计算值，这是因为计算时取用图里河气象站年降水量，加大了高降水区域的权重。这与图里河气象站降水量高于流域平均降水量的情况是一致的。

b.对于大桥屯水文站以上流域，克立格法计算值高于算术平均法与泰森多边形法计算的流域平均年降水量，这是因为计算时取用了牙克石市气象站、乌尔其汗水文站、维纳河降水量站的年降水量，加大了低降水区域的权重。这与三个测站降水量低于流域平均降水量的情况是一致的。

图6海拉尔河部分水文站以上流域

三种方法计算的流域平均降水量对比

Fig.6Comparison of average annual precipitation at some hydrological

stations in the Hailar River Basin obtained from three methods

c.对于牙克石水文站、红花尔基水文站、伊敏牧场水文站、坝后水文站、海拉尔水文站以上流域，克立格法计算值均高于算术平均法、泰森多边形法计算的流域平均降水量，这是因为参与计算的降水量站点多集中在流域中下游，而海拉尔河流域年降水量的总体变化趋势是自流域上游向下游逐渐减少，因而导致算术平均法、泰森多边形法计算值比实际偏低。（其他流域限于篇幅，略）

总体上克立格法与泰森多边形的计算值相差不大，而两种方法的计算结果与算术平均法相差较大，主要是因为海拉尔河流域各类型降水量站点尽管较多，但分布不均。因此，克立格法或泰森多边形法均可作为计算流域平均年降水量的选择，但克立格法计算出没有降水量站点分布的小流域平均降水量。

3结语

本文将多时段各向异性泛克立格空间估计理论引入流域平均年降水量的分析计算中，结果表明，多时段泛克立格模型充分考虑了水文数据空间上缺乏而时间上富足的特点，利用几站的水文时间序列资料即可建立起反映水文变量各时段空间变异性的变差函数，从而实现小流域或局部区域水文变量的最优估计问题。与国内常用的其他方法计算结果进行比较可知，多时段各向异性泛克立格法也可得出较理想的结果，而且适用于没有降水量站点分布的小流域，因此可以作为流域降水量计算的方法。

参考文献（References）：

［1］李恩羊.地质统计学在作物需水量研究中的应用［J］.武汉大学学报：工学版,1991（1）:311.(LI Enyang.On the Use of Geostatistics for Study of Crop Water Requirements［J］.Engineering Journal of WuHan University,1991（1):311..(in Chinese))

［2］袁新,李恩羊.参照作物需水量的空间变异性［J］.水利学报,1990(2):3539.(YUAN Xin,LI Enyang.Spatial Variability of Crop Water Requirements ［J］.Journal of Hydraulic Engineering,2009(2):810.(in Chinese))

［3］王仁铎,胡光道.线性地质统计学［M］.地质出版社,1988.(WANG Renduo,HU Guangdao.Linear Geostatistics ［M］.Beijing:Geology Publishing House,1988.(in Chinese))

［4］候景儒,黄竞先.地质统计学的理论和方法［M］.地质出版社,1990.(HOU Jingru,HUANG Jingxian.Theory and Method of Geostatistics［M］.Beijing:Geology Publishing House,1990.(in Chinese))

［5］雷志栋,杨诗秀.土壤特性空间变异性初步研究［J］.水利学报,2009(2)810.(LEI Zhidong,YANG Shixiu.Preliminary Investigation of the Spatial Variability of Soil Properties［J］.2009,2:810.(in Chinese))

［6］刘廷玺,朝伦巴根.多时段泛克立格空间估计理论及其在水文领域中的应用［J］.水利学报,1995（2）:7683.(LIU Tingxi,CHAO Lun bagen.The Theory of Universal Kriging Spatial Estimation with Multiple Time Period and its Application in Hydrology［J］.Journal of Hydraulic Engineering,1995（2):7683.(in Chinese))

推荐访问:海拉尔 变异性 降水 河流 地质