前面讲述的都是一维插值

前面讲述的都是一维插值,其特点是节点为一维实数组,插值函数是一元函数(曲线)。若节点是二维数组,则插值函数就是二元函数,即曲面。例如在某区域测量了若干点(节点)的高程(节点值),为了画出该区域较精确的等高线图,就要先插入更多的点(插值点),计算这些点的高程(插值),这就是一个二维插值问题。

二维插值问题的数学描述为:已知二元函数g(x, y)在某矩形区域R = [a,b]X[c,d]上互异节点(xi,yi)的函数值zij,如何求出在R 上任一点(x, y)处的函数值g(x, y)的近似值。

二维插值方法的基本思想是:根据已知数据点(xi,yi,zij),构造一个具有一定光滑性、简单易于计算的函数z = f (x, y)(已知曲面)作为z = g(x, y)(未知曲面)的近似,使曲面z = f (x, y)通过已知数据点,即f(xi,yi)=zij。然后计算点(x, y)的函数值f (x, y)作为g(x, y)的近似值。

常见的二维插值可分为网格节点插值和散乱节点插值。其中网格节点插值适用于节点比较规范的情况,即在包含所给节点的矩形区域内,节点由两族平行于坐标轴的直线的交点所组成。

散乱节点插值适用于一般的节点,多用于节点不太规范(即节点为两族平行于坐标轴的直线的部分交点)的情况。

1、网格节点插值法

网格节点插值问题可简述为:已知数据点(xi,yi,zij),其中i =1,2,…,m, j =1,2,…,

n,a=x1

网格节点插值法主要有以下几种形式:

(1)分片线性插值

分片线性插值即是在以(xi,yi),(xi+1,yi),(xi+1,yi+1),(xi,yi+1)为顶点的小矩形Rij上作线性插值。其分片插值函数为:

显然,分片线性插值函数在R上连续。

(2) 分片双线性插值

分片双线性插值曲面由一片一片空间二次曲面构成,其分片函数表达式为:

f (x, y) ??(Ax ??B)(Cx ??D), (x, y)?R?ij,∈

其中四个待定系数可由R ij四个顶点值唯一确定。但分片双线性插值曲面在R上不连续。

(3) 分片双三次样条插值

分片双三次样条插值函数在每一片(即小矩形R ij)上的表达式为:

F( x, y) =(A1 + A2 x + A3 x2 + A4 x3)(B1 + B2 x1 + B3 x2 + B4 x3)

其中待定系数Ai,Bi (i ??1,2,3,4)可由R ij四个顶点的值及插值函数f (x, y)在x和y 方向的光

' \' \

滑性(即f' x、fx、fxx、fy、fyy连续)和相应的边界条件所唯一确定。

2、散乱节点插值

散乱节点插值问题可简述为:已知函数z ??f (x, y)在矩形域R上散乱分布的n个互异节点(xi,yi)处的函数值zij,求R 上任一插值点(x, y)≠(xi,yi)的插值z。

常用的方法是反距离加权平均法,或称Shepard方法。其基本思想是,在点(x, y)≠(xi,yi),定义其插值函数的函数值为节点处函数值按(x, y)与节点距离的某种形式反比作为权重的加权平均,即某一点的函数值受周围各点的影响,较近的点影响较大,较远的点影响较小,其影响权数与距离平方成反比。例如若记

rij?则插值函数(曲面)可定义为

(x?xi)2?(y?yj)2

zij??f(x,y)??Wij(x,y)zij??i,j?rij?0rij?0rij 其中

1Wij(x,y)?2rij?ri,j12ij

这样定义的插值曲面是全局相关的,对曲面上任一点作数值计算都要涉及到全体已知数据,因而在已知数据量大的情况下计算量相当大,而且插值曲面在节点(xi,yi)处不光滑。虽然如此,由于该方法思想简单,故有种种改进方法。

另外一种常用的散乱数据插值方法是Kriging插值.它建立在地统计学理论基础上,利用区域化变量的原始数据和半方差函数的结构特征,对位采样点的区域化变量的取值进行线性最佳无偏估计.

Matlab中有一些计算二维插值的程序。如

z=interp2(x0,y0,z0,x,y,'method')

其中x0,y0分别为m 维和n维向量,表示节点,z0为nXm维矩阵,表示节点值,x,y为一维数组,表示插值点,x与y应是方向不同的向量,即一个是行向量,另一个是列向量,z为矩阵,表示得到的插值,'method'的用法同上面的一维插值。

如果是三次样条插值,可以使用命令

pp=csape({x0,y0},z0,conds,valconds),y=ppval(pp,{x,y})

其中x0,y0分别为m 维和n维向量,z0为mXn 维矩阵,具体使用方法同一维插值。

例2 在一丘陵地带测量高程,x和y 方向每隔100米测一个点,得高程如下表,试拟合一曲面,确定合适的模型,并由此找出最高点和该点的高程。

y x 100 200 300 400

clear,clc x=100:100:500; y=100:100:400;

100 636 698 680 662 200 697 712 674 626 300 624 630 598 552 400 478 478 412 334 500 450 420 400 310

z=[636 697 624 478 450 698 712 630 478 420 680 674 598 412 400 662 626 552 334 310]; pp=csape({x,y},z');

cz1=fnval(pp,{100:10:500,100:10:400});

cz2=interp2(x,y,z,100:10:500,[100:10:400]','spline'); figure;

mesh(100:10:500,100:10:400,cz2);

zi = griddata(x,y,z,100:10:500,[100:10:400]','v4'); figure;

mesh(100:10:500,100:10:400,zi);

[x y]=meshgrid(100:100:500,100:100:400); [xi,yi] = meshgrid(100:10:500,100:10:400);

[Vint]=gIDW(x(:),y(:),z(:),xi,yi,-2,'n',length(x(:))); figure;

mesh(100:10:500,100:10:400,Vint);

(2)kriging 插值

kriging 插值作为地统计学中的一种插值方法由南非采矿工程师

D.G.Krige于1951年首次提出,是一种求最优、线形、无偏的空间内插方法。在充分考虑观测资料之间的相互关系后,对每一个观测资料赋予一定的权重系数,加权平均得到估计值。

这里介绍普通Kriging插值方法的基本步骤:

1.该方法中衡量各点之间空间相关程度的测度是半方差,其计算公式为:

联系客服:779662525#qq.com(#替换为@) 苏ICP备20003344号-4