降雨空间插值分析

第五章降雨空间插值分析降雨空间插值分析是系统的中间件，其主要任务是把流域内175个雨量站的资料利用空间插值方法合理地插值到分布式水文模型所应用的空间网格上，以便于利用历史和实时自动测报雨量进行模型的率定和模拟验证，其输出结果以数据库或数据文本方式储存。

该层次的功能主要包括以下三个方面：（1）把175站雨量信息合理地插值到计算网格；（2）雷达降雨与分布式水文模型耦合接口；（2）暴雨数值预报与分布式水文模型耦合接口。

5.1 概述降雨空间插值分析是本系统的关键技术之一。

该部分的功能实现途径如下：（1）建立统一的基础空间数据库，包括统一的网格、单元、区域、子流域划分及编码，实现对同一区域对象的地理、水文、气象综合描述。

系统的基本分辨率规定为空间1km×1km；（2）多源降雨信息的同化及整合。

无论是自动测报实时雨量，还是历史数据（包括月、日、时等时段），通过该软件都可以生成网格上的空间分布数据。

（3）数值天气预报产品转化为1km网格的空间数据；（4）雷达信息转换为分布式水文模型所用网格的空间数据。

5.2 空间插值方法空间插值方法的主要思想是：由分布的流域上的各个测站（xi, yi, zi ）（x, y 为坐标值，z 为雨量值），拟合出该时段降雨量在流域上的分布函数f (x, y)，进而求得在该函数在计算网格上的积分：()⎰⎰=dA y x f P , 5-1则网格上的面平均雨量为：A P P = 5-2在实际操作时，分布函数的拟合是采用加权的最小二乘拟合得出，但是对于复杂的空间分布函数，其求解并不是简单的问题。

一般情况下多选用多项式函数来作为数学表达式，另外还要求解上的可行性和便利性，目前趋势面的求解均采用最小二乘法，一般来说只有线性表达式以及可转化为线性的表达式方可求解。

目前流行较多的方法有：算术平均、距离反比加权平均、最短距离法、空间函数拟合插值等。

算术平均方法比较简单，如果网格内有雨量站点，则该网格内的平均雨量为网格内站点雨量的平均值，但是小花间网格要4万多个，而雨量站点165个，该方法不能适用。

以下重点介绍距离反比加权平均、最短距离法、克里格法和空间函数拟合插值方法。

5.2.1. 距离反比加权插值距离反比加权插值（FIDW ，Inverse Distance Weighted Interpolation ）是根据网格中心点附近的雨量站点资料插值网格中心点雨量（图5.1），以此代表网格上的平均雨量，计算公式如下：∑∑=-=-=r1p n p r 1p n p p y x,d d ZZ 5-3式中：Z p 是相邻点的高程，d 是插值点到p 点的距离；n 是参数，范围从1.0到6.0，通常用的值是2.0。

-n 表示越靠近被插值点越重要。

图5.1 FIDW 方法计算网格平均雨量示意图网格跨过边界的插值生成的插值面与实际不一致；插值面可能在主要的类型、区域和分类内局部有效。

FIDW 方法的一种假象就是帐篷支柱影响。

也就是局部的最大值和最小值都位于测量点的位置。

当用雨量计测量时，给人的印象是在雨量站雨的强度最大，这显然是不合理的。

图5.2示例了用IDW 方法通过雨量站的点雨量插值降雨等值面，这个面是通过距离平方反比方法用所有的数据点插值得到的。

等值面显示了降雨主要发生在雨量站的周围。

图5.2 距离反比加权插值方法结果示意图5.2.2 最短距离法最短距离法是比较简单的一种方法，即用与网格中心点最近的雨量站点资料代表网格平均雨量。

当网格趋近于无穷小时，该方法与水文学中常用的泰森多边形相同。

该方法的缺陷是雨量空间分布在两个雨量站控制分界线处不连续。

5.2.3 克里格法（Kriging）克里格法（Kriging）是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法。

克里格法的适用条件是区域化变量存在空间相关性。

应用克里格法首先要明确三个重要的概念。

一是区域化变量；二是协方差函数，三是变异函数。

（1）区域化变量。

当一个变量呈空间分布时，就称之为区域化变量。

这种变量反映了空间某种属性的分布特征。

矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属性。

区域化变量具有双重性，在观测前区域化变量Z(X)是一个随机场，观测后是一个确定的空间点函数值。

区域化变量具有两个重要的特征。

一是区域化变量Z(X)是一个随机函数，它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点X 与偏离空间距离为h 的点X ＋h 处的随机量Z(X)与Z(X+h)具有某种程度的自相关，而且这种自相关性依赖于两点间的距离h 与变量特征。

在某种意义上说这就是区域化变量的结构性特征。

（2）协方差函数。

协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。

在概率理论中，随机向量X 与Y 的协方差被定义为：()()()[]EY Y EX X E Y X Cov --=, 5-4区域化变量()()w v u x x x Z X Z ,,= 5-5在空间点x 和x+h 处的两个随机变量Z(x)和Z(x+h)的二阶混合中心矩定义为Z(x)的自协方差函数，即()()[]()()[]()[]()[]h X Z E X Z E h X Z X Z E h x Z x Z Cov +-+=+, 5-6 区域化变量Z(x)的自协方差函数也简称为协方差函数。

一般来说，它是一个依赖于空间点x 和向量h 的函数。

设Z(x)为区域化随机变量，并满足二阶平稳假设，即随机函数Z(x)的空间分布规律不因位移而改变，h 为两样本点空间分隔距离或距离滞后，Z(x i )为Z(x)在空间位置x i 处的实测值，Z(xi+h)是Z(x)在x i 处距离偏离h 的实测值[i=1,2,…,N(h)]，根据协方差函数的定义公式，可得到协方差函数的计算公式为：()()()()[]()()[]()∑=+-+-='h N i i i i i h x Z h x Z x Z x Z h N h c 11 5-7在上面的公式中，N(h)是分隔距离为h 时的样本点对的总数，()i x Z 和()h x Z i +分别为()i x Z 和()h x Z i +的样本平均数，即∑==N i i i x Z N x Z 1)(1)(，∑=+=+N i i i h xZ N h x Z 1)(1)( 5-8在公式中N 为样本单元数。

一般情况下)()(h x Z x Z i i +≠，特殊情况下可以认为近似相等。

若m h x Z x Z i i =+=)()(（常数），协方差函数可改写为如下：[]∑=-+=)(12^)()()(1)(h N i i i m h x Z x Z h N h c 5-9式中：m 为样本平均数，可由一般算术平均数公式求得，即∑==n i ix Z N m 1)(1 5-10（3）变异函数。

变异函数又称变差函数、变异矩，是地统计分析所特有的基本工具。

在一维条件下变异函数定义为，当空间点x 在一维x 轴上变化时，区域化变量Z(x)在点x 和x+h 处的值Z(x)与Z(x+h)差的方差的一半为区域化变量Z(x)在x 轴方向上的变异函数，记为γ(h)，即：[])()(21),(h x Z x Z V h x ar +-=γ 5-11[][][]{}22)()(21)()(21h x Z E x Z E h x Z x Z E +--+-=在二阶平稳假设条件下，对任意的h 有，()[]()[]h X Z E X Z E +=，因此上式可以改写为：[]2)()(21),(h x Z x Z E h x +-=γ 5-12从上式可知，变异函数依赖于两个自变量x 和h ，当变异函数γ(x,h)仅仅依赖于距离h 而与位置x 无关时，可改写成γ(h)，即[]2)()(21)(h x Z x Z E h +-=γ 5-13设Z(x)是系统某属性 Z 在空间位置x 处的值，Z(x)为一区域化随机变量，并满足二阶平稳假设，h 为两样本点空间分隔距离，Z(xi)和Z(xi+h)分别是区域化变量在空间位置xi 和xi+h 处的实测值[i=1,2,...,N(h)]，那么根据上式的定义，变异函数γ(h)的离散公式为：[]2)(1^)()()(21)(∑=+-=h N i i i h x Z x Z h N h γ 5-14变异函数揭示了在整个尺度上的空间变异格局，而且变异函数只有在最大间隔距离1/2处才有意义。

（4）克里格估计量。

假设x 是所研究区域内任一点，Z(x)是该点的测量值，在所研究的区域内总共有n 个实测点，即x 1,x 2,...,x n ，那么，对于任意待估点或待估块段V 的实测值Z v (x)，其估计值是通过该待估点或待估块段影响范围内的n 个有效样本值()()n i x Z i v ,...,2,1=的线性组合来表示，即()()∑=*=n i i i vx Z x Z 1λ 5-15式中：λi 为权重系数，是各已知样本在Z(xi)在估计值时影响大小的系数，而估计值的好坏主要取决于怎样计算或选择权重系数λi在求取权重系数时必须满足两个条件，一是使估计值的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值和实际值Z v (x)之差的平方和最小，在数学上，这两个条件可表示为()()[]0=-*x z x Z E v v ()()[]()()[]min 2→-=-**x z x Z E x z x Z Var v v v v 5-16（5）普通克里格分析方法。

设Z(x)为区域化变量，满足二阶平稳和本征假设，其数学期望为m ，协方差函数c(h)及变异函数λ(h)存在。

即()[]m x Z E =()()()[]2m h x Z x Z E h c -+= ()()()[]221h x Z x Z E h +-=γ 5-17对于中心位于x0 的块段为V ，其平均值为Z v (x 0)的估计值以()()⎰=V V dx x Z V h Z 1 5-18进行估计。

在待估区段V 的邻域内，有一组n 个已知样本v(x i )(i=1,2,…n)，其实测值为Z(x i )(i=1,2,…n)。

克里格方法的目标是求一组权重系数λi (i=1,2,…n)，使得加权平均值：()()∑=*=n i i i Vx Z x Z 1λ 5-19成为待估值段V 的平均值Z v (x 0)的线性、无偏最优估计量，即克里格估计量。

为此，要满足以下两个条件：1、无偏性。

要使()x Z V*成为Z v (x) 的无偏估计量，即[][]V V Z E Z E =*，当[]m Z E V =*时，则有：11=∑=n i i λ。