降雨空间插值分析.doc

第五章 降雨空间插值分析降雨空间插值分析是系统的中间件，其主要任务是把流域内175个雨量站的资料利用空间插值方法合理地插值到分布式水文模型所应用的空间网格上，以便于利用历史和实时自动测报雨量进行模型的率定和模拟验证，其输出结果以数据库或数据文本方式储存。该层次的功能主要包括以下三个方面：（1）把175站雨量信息合理地插值到计算网格；（2）雷达降雨与分布式水文模型耦合接口；（2）暴雨数值预报与分布式水文模型耦合接口。5.1 概述降雨空间插值分析是本系统的关键技术之一。该部分的功能实现途径如下：（1）建立统一的基础空间数据库，包括统一的网格、单元、区域、子流域划分及编码，实现对同一区域对象的地理、水文、气象综合描述。系统的基本分辨率规定为空间1km1km；（2）多源降雨信息的同化及整合。无论是自动测报实时雨量，还是历史数据（包括月、日、时等时段），通过该软件都可以生成网格上的空间分布数据。（3）数值天气预报产品转化为1km网格的空间数据；（4）雷达信息转换为分布式水文模型所用网格的空间数据。5.2 空间插值方法空间插值方法的主要思想是：由分布的流域上的各个测站（xi, yi, zi）（x, y为坐标值，z为雨量值），拟合出该时段降雨量在流域上的分布函数f (x, y)，进而求得在该函数在计算网格上的积分：5-1则网格上的面平均雨量为：5-2在实际操作时，分布函数的拟合是采用加权的最小二乘拟合得出，但是对于复杂的空间分布函数，其求解并不是简单的问题。一般情况下多选用多项式函数来作为数学表达式，另外还要求解上的可行性和便利性，目前趋势面的求解均采用最小二乘法，一般来说只有线性表达式以及可转化为线性的表达式方可求解。目前流行较多的方法有：算术平均、距离反比加权平均、最短距离法、空间函数拟合插值等。算术平均方法比较简单，如果网格内有雨量站点，则该网格内的平均雨量为网格内站点雨量的平均值，但是小花间网格要4万多个，而雨量站点165个，该方法不能适用。以下重点介绍距离反比加权平均、最短距离法、克里格法和空间函数拟合插值方法。5.2.1. 距离反比加权插值距离反比加权插值（FIDW，Inverse Distance Weighted Interpolation）是根据网格中心点附近的雨量站点资料插值网格中心点雨量（图5.1），以此代表网格上的平均雨量，计算公式如下：5-3式中：Zp是相邻点的高程，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的点Xh处的随机量Z(X)与Z(X+h)具有某种程度的自相关，而且这种自相关性依赖于两点间的距离h与变量特征。在某种意义上说这就是区域化变量的结构性特征。（2）协方差函数。协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。在概率理论中，随机向量X与Y的协方差被定义为：5-4区域化变量5-5在空间点x和x+h处的两个随机变量Z(x)和Z(x+h)的二阶混合中心矩定义为Z(x)的自协方差函数，即5-6区域化变量Z(x)的自协方差函数也简称为协方差函数。一般来说，它是一个依赖于空间点x和向量h的函数。设Z(x)为区域化随机变量，并满足二阶平稳假设，即随机函数Z(x)的空间分布规律不因位移而改变，h为两样本点空间分隔距离或距离滞后，Z(xi)为Z(x)在空间位置xi处的实测值，Z(xi+h)是Z(x)在xi处距离偏离h的实测值i=1,2,N(h)，根据协方差函数的定义公式，可得到协方差函数的计算公式为：5-7在上面的公式中，N(h)是分隔距离为h时的样本点对的总数，和分别为和的样本平均数，即，5-8在公式中N为样本单元数。一般情况下，特殊情况下可以认为近似相等。若（常数），协方差函数可改写为如下：5-9式中：m为样本平均数，可由一般算术平均数公式求得，即5-10（3）变异函数。变异函数又称变差函数、变异矩，是地统计分析所特有的基本工具。在一维条件下变异函数定义为，当空间点x 在一维x 轴上变化时，区域化变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差的一半为区域化变量Z(x)在x轴方向上的变异函数，记为(h)，即：5-11在二阶平稳假设条件下，对任意的h有，因此上式可以改写为：5-12从上式可知，变异函数依赖于两个自变量x和h，当变异函数(x,h)仅仅依赖于距离h 而与位置x无关时，可改写成(h)，即5-13设Z(x)是系统某属性 Z在空间位置x处的值，Z(x)为一区域化随机变量，并满足二阶平稳假设，h为两样本点空间分隔距离，Z(xi)和Z(xi+h)分别是区域化变量在空间位置xi和xi+h处的实测值i=1,2,.,N(h)，那么根据上式的定义，变异函数(h)的离散公式为：5-14变异函数揭示了在整个尺度上的空间变异格局，而且变异函数只有在最大间隔距离1/2处才有意义。（4）克里格估计量。假设x是所研究区域内任一点，Z(x)是该点的测量值，在所研究的区域内总共有n个实测点，即x1,x2,.,xn，那么，对于任意待估点或待估块段V的实测值Zv(x)，其估计值是通过该待估点或待估块段影响范围内的n个有效样本值的线性组合来表示，即5-15式中：i为权重系数，是各已知样本在Z(xi)在估计值时影响大小的系数，而估计值的好坏主要取决于怎样计算或选择权重系数i在求取权重系数时必须满足两个条件，一是使估计值的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值和实际值Zv(x)之差的平方和最小，在数学上，这两个条件可表示为5-16（5）普通克里格分析方法。设Z(x)为区域化变量，满足二阶平稳和本征假设，其数学期望为m，协方差函数c(h)及变异函数(h)存在。即5-17对于中心位于x0 的块段为V，其平均值为Zv(x0)的估计值以5-18进行估计。在待估区段V的邻域内，有一组n个已知样本v(xi)(i=1,2,n)，其实测值为Z(xi)(i=1,2,n)。克里格方法的目标是求一组权重系数i(i=1,2,n)，使得加权平均值：5-19成为待估值段V的平均值Zv(x0)的线性、无偏最优估计量，即克里格估计量。为此，要满足以下两个条件：1、无偏性。要使成为Zv(x) 的无偏估计量，即，当 时，则有：。这时，是的无偏估计量。2、最优性。在满足无偏性条件下，估计方差为5-20由方差估计可知5-21为使估计方差最小，根据拉格朗日乘数原理，令估计方差的公式为：5-22求以上公式对i和的偏导数，并令其为0，得克里格方程组5-23整理后得：5-24解上述n+1阶线性方程组，求出权重系数i和拉格朗日乘数，并带入公式，经过计算可得克里格估计方差，即：5-25以上三个公式都是用协方差函数表示的普通克里格方程组和普通克里格方差。5.2.4 空间函数拟合插值方法可用作空间插值的分析方法有：趋势面分析和曲面插值分析。趋势面分析拟合的曲面不要求所拟合的曲面通过抽样点，但要求所拟合的曲面在整体上要充分逼近实际的空间曲面。通过插值方法拟合的曲面要求所拟合的曲面必须要通过所有抽样的数据点。趋势面分析必然存在两方面的问题需要考虑：一是数学曲面类型的确定，二是拟合精度。一般情况下多选用多项式函数(表5-1)来作为数学表达式，另外还要求解上的可行性和便利性，目前趋势面的求解均采用最小二乘法，一般来说只有线性表达式以及可转化为线性的表达式方可求解。图5.3为趋势面分析例子。表5-1 用于插值分析的通用多项式独立项项次表面性质项数Z=a00平面1A1X+a2Y1线性2a3XY+a4X2+a5Y22二次抛物面3a6X3+a7Y3+a8XY2+a9X2Y3三次曲面4a10X4+a11Y4+a12X3Y+a13X2Y2+a14XY34四次曲面5图5.3 趋势面分析示例图5.3中，a是原始数据抽样方案，b、c、d是根据a所给出的抽样点数据计算出的N1，2，3时多项式趋势面的等值线表示。一次趋势面反映了空间分布由北到南的渐增趋势，二次趋势面反映了这种趋势在中部开始分化为向东南和西南渐增，三次趋势面是对二次趋势面的进一步修饰，反映出渐增趋势的具体走向，从图中可以看出，三次趋势面已经基本接近于实际分布，但舍去了原等值线的细小弯曲。曲面插值分析可分为全局插值和分块插值两种。全局插值就是用一个统一点数学曲面来描述待插值曲面，这个统一的曲面必须在全部抽样点上与抽样数据吻合。从数学上讲，要找一个简单曲面并使其通过全部抽样点是不可能的。为了解决这个问题，Harly提出了用多个简单曲面叠加进行全局插值的多面函数插值法，但这种方法也仅适用于小范围插值。分块插值是指将整个分析区域，根据数据格网剖分为若干个小单元，每个单元上独立地拟合一个小曲面片，若干小曲面片连接起来构成整个空间曲面。一般地，要求相邻曲面片在边界上连续，即给出相同的数值，特殊情况还要求由若干小曲面片所连接起来的空间曲面是一个光滑曲面，以便更逼真地描述曲面形态。基于正方形网格的分块插值方法有：双线性多项式曲面插值和双三次多项式曲面插值。双线性多项式曲面插值公式为：5-26此模型具有4个待定系数，可以通过正方形网格的4个顶点的抽样数据值确定。对于单元网格，则有以下计算公式5-27式中：下标p为插值点，下标1、2、3、4分别为（1，0）、（1，1）、（0，1）和（0，0）4个点（图2.16）。双线性多项式曲面插值方法拟合的曲面是连续的，但并不光滑（图5.