实验报告五学院名称：理学院 专业年级： 姓 名： 学 号：课 程：数学模型与数学建模 报告日期：2015年12月8日一、实验题目例2.2.1 水库库容量与高程设一水库将河道分为上、下游两个河段，降雨的开始时刻为8时，这是水位的高程为168m ，水库容量为38109.21m ⨯，预测上游的流量()()s m t Q /3，d 取值如表2.2.1所示。

表2.2.1 上有流量()t Q 的预测已知水库中水的容量()3810mV 与水位高程H （m ）的数值关系为表2.2.2表2.2.2 水库库容量与水位高程的关系如果当日从8时开始，水一直保持s m /10003的泄流量，根据所给数据，预报从降雨时刻到56h 以内每小时整点时刻水库中水的库容量与水位高程。

例2.2.2 地下含沙量某地区有优质细沙埋在地下，某公司拟在此处采沙，已得到该地区钻探资料图的一角如下表，在每个格点上有三个数字列，都是相对于选定基点的高度（m ），最上面的数字是覆盖表面的标高，中间的数字是沙层顶部的标高最下面的数字是沙层底部的标高，每个格子都是正方形，边长50m 。

画星号处，即沼泽表层地带，没有钻探数据。

试估计整个矩形区域内的含沙量。

二、实验目的插值模型是数据挖掘的另一类模型，插值（Interpolation ）的目的是根据能够获得的观测数据推测缺损的数据，此时观测数据(){}ni i i y x 1,=被视为精确的基准数据，寻找一个至少满足条件的函数()x y y =，使得()n i x y y i i ,,2,1,Λ==，在本节我们强调的是插值模型的应用，而不是插值方法的构造。

三、问题陈述2.2.1 一维插值例2.2.1 水库库容量与高程 2.2.2 二维插值例2.2.2 地下含沙量 2.2.3 泛克里金插值四、模型及求解结果2.2.1 一维插值一元函数差值公式为()()∑==ni i i x y x y 1λ其中()x i λ是满足条件()iji x δ=λ的函数，依据插值的公式，如最近邻差值，线性插值、分段三次Hermite插值等，分别取阶梯函数、线性函数、三次多项式函数等，相应的数学表达式可以查阅本科生数值计算教材。

下面先通过简单的Matlab一维插值令interpret1了解相应的计算结果。

例：2.2.1 水库库容量与高程为了给出每小时的报告，需要补充每小时整点时刻上有流量的数据，以及相应不同库容量的水位高程。

假设（a）已知数据准确（b）相邻两个时刻之间的流量变化是现行的（c）相邻两个水位高程之间的高程对水的库容量的变化也是线性的首先，利用Matlab线性插值令，确定每小时的上游流量q（t）.由程序在Matlab中运行的结果为：q =1.0e+004 *Columns 1 through 90.3600 0.4050 0.4500 0.4950 0.5400 0.6000 0.6600 0.7200 0.7800Columns 10 through 180.7975 0.8150 0.8325 0.8500 0.8675 0.8850 0.9025 0.9200 0.9350Columns 19 through 270.9500 0.9650 0.9800 0.9950 1.0100 0.9629 0.9157 0.86860.8214Columns 28 through 360.7743 0.7271 0.6800 0.6329 0.5857 0.5386 0.4914 0.4443 0.3971Columns 37 through 450.3500 0.3000 0.2500 0.2410 0.2320 0.2230 0.2140 0.2050 0.1960Columns 46 through 490.1870 0.1780 0.1690 0.1600然后确定每个时刻t 的水库容量()t v ，因为，水库容量=原库存量+流入量-泄流量()s m /1038，即：()()()81036001036009.21588-⨯-⨯+=--⎰t ds s q t v t这里我们遇到数值积分，被积函数()s q 没有解析表达式，只有一个数列表示，iq 表示在i整点时刻的流量，利用Matlab 逐点积分指令()y x cumptrapz ,，可以得到水库容量()t v v =在每一刻[]56,8∈t 的值。

最后确定每时刻t 水库的水位高程h （t ），因为最大水库容量已经远远超出了已知数据范围，需要利用外插方法补充数据，确定水库高程对水库容量的依赖关系h=H （v ）。

最后利用函数复合得到水位高程()())(t v H t h = 确定每小时的上游流量q （t ）利用函数复合得到水位高程()())(t vH t h=2.2.2 二维插值二维插值大致可以分为两种，规则点插值和散乱点插值，前者即利用在方形网格点),(iiyx给定的数据mj n i z ij ,,2,1,,,2,1,ΛΛ==，在加密网格点上补充相应函数值，插值公式为∑∑===ni mj ij ij y x l z y x z 11),(),(其中ijl 是满足jsik s k ij y x l δδ=),(的二元函数，如双线性函数、双三次多项式函数、样条函数等，一句不同的插值方式选取，相应的表达式可以从本科生数值计算教材中查阅。

这里先通过Matlab 二维插值指令interp2了解相应的计算结果散乱点插值是要利用在不规则排列的观测点),(j j y x 上的数据nj z j ,,2,1,Λ=确定其他点上的函数值，插值公式形式为ijj j i nj j j y x l y x l z y x z δ==∑=),(,),(),(1常用到的插值方法，如距离加权反比插值，Kriging 插值，需要查阅专业书籍，下面先通过Matlab 二维插值指令griddata 了解相应的计算结果例2.2.2 地下含沙量假设地下沙层连续变化，于是可以利用积分运算矩形区域内的含沙量。

首先通过插值补充缺失的数据，采用一维样条插值。

为了提高梯形积分精度，利用二维样条插值加密数据，得到边长为0.5m的长方形网格上的数据。

最后，将二重积分化为累次积分，利用梯形公式得到积分的近似值。

v =6.4290e+005如果在一维插值补齐缺失数据中用线性插值代替上面的样条插值，将得到沙层体积6345253m ，和以上结果有些差别。

对于加密二维网格，应该加到多密？注意，从理论上讲，网格越密计算精度越高，但是从计算角度看，网格越密计算误差累计也越大，所以需要通过逐步加密网格，确定精度最优的积分值。

如果直接利用二维插值，例如采用Matlab 的散乱点插值指令，仍用二次梯形公式，会得到相近的结果。

v =6.3572e+0052.2.3 泛克里金插值考虑到沙子是一种特殊的地址，下面试用地质学常用的泛克里金插值方法计算。

先介绍克里金插值方法，这本身就是运用统计学方法建模的一个有趣的例子。

1951年南非地质学家Krige 讲()221,R x x x ∈=处矿藏的贮藏量()x f 看成是随机函数()x F 的一个实现，提出了依据观测值()n j y x j j ...,3,2,1,,=，寻找基函数(){}n j x j ...,2,1,=λ，以获得()x F 的具有形式 ()()()∑=jjjx F x x F λ*的最小方差无偏估计，取()x F *的条件期望()()()()n j f x F x F E x fj j ....,2.1,|**===给出未知点x 估值的插值方法，即考虑形如()()()()()x R x P x R x M x F K+=+=∑=1αααφ的随机函数，其中()x M 在地质学中称为飘移，在一些随机分析的场合也成为趋势，αP 是均值为ααp P E =)(的未知的随机变量，}{αφ是多项式基函数，例如这里取1)(1=x φ，12)(x x =φ，23)(x x =φ，214)(x x =φ，225)(x x =φ，216)(x x x =φ，()x R 是满足0))((=x R E ，)())()((h h x R x R E σ=+的随机函数)(h σ是给定的且满足当;0)(→∞→h h σ时，有当1)(0→→h h σ时，有的核函数。

这里按通常取法，取高斯核函数()2ex p )(h h -=σ根据无偏估计要求，()()()()∑===1*j jjx EF x x EF x EF λ即：()()()∑∑∑====K nj JjK X P x x P 111ααααααφλφ 从而得到无偏条件，对任意的α有：()()()∑=jj j x x x ααφλφ在这个约束条件下极小化方差 ()()()∑-jjjx F x x F E 2)(λ=()()()()()()∑∑∑∑---ααααααφλφjjjjjx P x F x x x P x F E 2))()((=()()()∑-jjjx R x x R E 2)(λ=()()()()()()()()()∑∑∑+-jj j jkk j k j x R E x R x R E x R x R E x 22λλλ =：)0()()(2)()(σ+Λ-ΛΛx D x x A x TT其中n n k j x x A ⨯-=))((σ，1))(()(⨯-=n j x x x D σ，1))(()(⨯=Λn j x x λ，1)0(=σ。

引入拉格朗日乘子()K γγγ...,21Γ，由拉格朗日函数()()()()()()∑∑==⎪⎪⎭⎫⎝⎛--+Λ-ΛΛ=ΓΛKn j j j TTx x x x D x x A x L 11)()(212ααααφλφγ：， 的极小值点必为其驻点的结论，得到关系式 ()()()⎪⎪⎭⎫ ⎝⎛=⎪⎪⎭⎫ ⎝⎛ΓΛ⎪⎪⎭⎫⎝⎛ΦΦx x D x AT φ0 其中n n k j x x A ⨯-=))((σ，,))(()(K n j x x ⨯=Φαφ1))(()(⨯=Λn j x x λ，,)(1⨯=ΓK αγ1))(()(⨯-=n j x x x D σ，,))(()(1⨯=K x x αφφ于是，在未知点x 上的随机函数值()x F 可以用它的最小方差无偏估计的条件期望()()()()n j f x F x F E x f j j ....,2.1,|**===近似，从而得到泛克里金插值函数： ()()()∑⎪⎪⎭⎫⎝⎛ΓΛ=Λ==j T T T j j f x f x f x x f 0)()(*λ =()⎪⎪⎭⎫⎝⎛⎪⎪⎭⎫ ⎝⎛ΦΦ-00~)()(1f A x x D TTT φ =()g x x D TT )()(φ 其中()1,1)0()0(~,⨯⨯⨯===K K K n jO O f f 练习：证明上面构造的泛克里金函数是连续的，且通过已知点，即插值函数满足： n j f x f j j ,...,2,1,)(==*用泛克里金插值得到的边长为0.5m 的加密网格上的沙层厚度，通过二次梯形公式，得到沙层体积近似值为6373203m五、程序代码2.2.1 一维插值>> x=0:4:20;>> y=[37 51 45 74 83 88]; >> xx=0:1:20;>> y1=interp1(x,y,xx,'nearest'); %最近邻插值，间断函数 >> y2=interp1(x,y,xx); %线性插值，连续函数>> y3=interp1(x,y,xx,'cubic'); %分段三次Hermite 插值，一阶连续函数 >> y4=interp1(x,y,xx,'splinet'); %样条插值，二阶倒数连续的函数 >> subplot(2,2,1),plot(x,y,'kd',xx,y1),title('nearest') >> subplot(2,2,2),plot(x,y,'kd',xx,y2),title('linear') >> subplot(2,2,3),plot(x,y,'kd',xx,y3),title('cubic') >> subplot(2,2,4),plot(x,y,'kd',xx,y4),title('spline')例2.2.1 水库库容量与高程利用Matlab线性插值令，确定每小时的上游流量q（t）>> T=[8,12,16,24,30,44,46,56];>> Q=[3600,5400,7800,9200,10100,3500,2500,1600];>> t=8:56;>> q=interp1(T,Q,t,'linear')%得到每小时上游流量>> plot(T,Q,'kd',t,q),title('time-flow')水库容量()t vv=在每一刻[]56,8∈t的值>> v=21.9+36*10^(-6)*cumtrapz(t,q-10^3*ones(size(t))); %每小时水库容量>> vmax=max(v); %得到最大的水库容量为30.7524（10^8立方米）利用函数复合得到水位高程()())(t vH t h=>> V=[21.9 23.93 24.06 24.12 24.33 24.47 24.3 24.75]; >> H=[168 168.75 168.8 168.85 168.9 168.95 169 169.05]; >> h=interp1(V,H,v,'linear','extrap'); %得到每小时水位高程>> hmax=max(h);%最大水位高程171.0508米>> plot(V,H,'kd',v,h),title('volume-atitude')2.2.2 二维插值二维插值指令>> x=0:4:16;y=0:4:16;>> z=[620 730 800 850 870;760 880 970 720 1050880 1080 630 1250 1280;980 1180 1320 1450 14201060 1230 1390 1500 1500];>> [X,Y]=meshgrid(0:16,0:16);>> Z1=interp2(x,y,z,X,Y,'nearest');Z2=interp2(x,y,z,X,Y); >> Z3=interp2(x,y,z,X,Y,'cubic');Z4=interp2(x,y,z,X,Y,'spline'); >> subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title('spline')散乱点插值>> x=[1 2 3 4 5];>> y=[4 1 3 2 5];>> z=[4 1 6 2 5];>> [X,Y]=meshgrid(0:0.2:5);>> Z1=griddata(x,y,z,X,Y,'nearest');>> Z2=griddata(x,y,z,X,Y);>> Z3=griddata(x,y,z,X,Y,'cubic');>> Z4=griddata(x,y,z,X,Y,'v4');>> subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title('v4')例2.2.2 地下含沙量通过一维插值补充缺失的数据>> D=[23.0,23.1,23.2,23.4,23.5,24.0,24.0,24.0;23.1,23.3,23.4,23.4,23.5,24.2,24.1,24.1];>> E=[19.9,20.0,20.0,19.8,19.9,20.0,19.8,19.6;19.8,19.7,19.4,20.0,20.1,20.3,20.3,20.5];>> F=[6.0,3.2,1.6,1.0,1.1,1.0,0.8,0.9;2.2,1.4,0.6,0.5,0.3,-0.2,-0.1,0.0];>> P=[1,6,7,8];K=1:8;>> A=[22.4,22.5,23.0,23.2];>> A1=interp1(P,A,K,'spline');a=[A1;D]>> B=[20.0,18.4,17.8,18.0];>> B1=interp1(P,B,K,'spline');b=[B1;E]>> C=[5.8,0.5,0.4,0.4];>> C1=interp1(P,C,K,'spline');c=[C1;F]二维样条插值加密数据，得到边长为0.5m的长方形网格上的数据>> [M,N]=meshgrid(1:8,0:2);>> delta=0.01;[X,Y]=meshgrid(1:delta:8,0:delta:2);>> a1=interp2(M,N,a,X,Y,'spline');>> b1=interp2(M,N,b,X,Y,'spline');>> c1=interp2(M,N,c,X,Y,'spline');>> mesh(X,Y,a1),hold on,>> mesh(X,Y,b1),hold on,mesh(X,Y,c1)最后，将二重积分化为累次积分，利用梯形公式得到积分的近似值。