计算方法实验报告

关于最小二乘拟合曲线在重力异常资料处理中的应用的试验报告
常用的重力异常资料处理方法有以下类型：圆滑、平均法、趋势分析、空间延拓、导数换算等。

其中圆滑法又分为平均圆滑和多项式圆滑。

平均圆滑能够有效地消除随机误差和干扰，但圆滑后有可能改变异常形态特征，给一些利用异常性态特征进行异常解释的工作带来困难。

多项式圆滑是假设局部区域内有效异常形态能用多项式函数拟合，将偏离这个函数的部分视为误差或干扰。

(1) 剖面数据多项式圆滑
假设在一个局部区域内异常圆滑值可用多项式
表示，则有
在x 附近的局部区域[x-M Δx, x+M Δx]内圆滑异常用该函数拟合。

若取 N 阶多项式拟合异常，有
令
这样，共有N+1个未知变量，2M+1个方程，写成矩阵形式为
显然，多项式阶次越高，对异常数据拟合越好。

一般情况下，希望圆滑值不要过于逼近原始值，否则，不能取得圆滑效果。

另一方面，为了保证方程组可解，通常要求N≤ 2M ，即保持正定或超定。

实际中，位场异常数据的多项式圆滑阶次很少超过3次，即N≤3。

以五点二次（即M=2, N=2）圆滑为例，如此所构成的是一个超定方程组，可用最小二乘方法求解。

+++=2210)(x a x a a x f +++=∆2210)(x a x a a x g []M M M M i x i x a x i x a x i x a a x i x g N N ),1(,,2,1,0,1,),1(,)()()()(2210-----=∆⋅+++∆⋅++∆⋅++=∆⋅+∆ ∑
+=-=∆=+=∆-+-=111,)(12,,2,1,)1(N j j i ij i i i x u x g g M i x i M x x ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡==⋅++++++++)1)(12(1)12()1(22221)1(11211)1(10)12(1221,,N M M N N N N M M u u u u u u u u a a a g g g U X G G X U
根据最小二乘法，有
考虑到数据点的对称分布，上式改写成
有
实际计算时，采用窗口内的异常点进行圆滑，若将坐标作一下替换
则有
显然，每移动一次计算窗口，仅x=0时的圆滑值是所需计算，其它坐标上圆滑不需要考虑。

从上式可知，当x=0时，只需解出a0即可得到窗口中心点的圆滑值，即
五点二次圆滑计算公式为
同理，可以导出三点二次、七点二次、七点三次九点三次圆滑计算公式。

无论几点圆滑公式，事实上是计算点周边点上异常值的加权之和。

2525105222210221
21101)()()(x a x a a x g x a x a a x g x a x a a x g ++=∆++=∆++=∆ min )]([5
122210→∆-++=Φ∑=i i i i x g x a x a a min ][22
22210→-++=Φ∑-=i i i i g x a x a a ⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=++=++⇒⎪⎪⎪⎪⎭
⎪⎪⎪⎪⎬⎫=-++=∂Φ∂=-++=∂Φ∂=-++=∂Φ∂∑
∑∑∑∑∑∑∑∑∑∑∑∑∑-=-=-=-=-=-=-=-=-=-=-=-=-=-=2222242223122202222322221220222222221022222102222210122221000][20][20][2i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i x g x a x a x a x g x a x a x a g x a x a a x g x a x a a a x g x a x a a a g x a x a a a 2,1,0,1,2++--=x ⎪⎩⎪⎨⎧+++--=+++++--=++++++=++++--++--++--)44(34010)22(0100)(10021012202101212101220g g g g g a a g g g g g a g g g g g a a )](3)(1217[351)0(221100+-+-+-++==∆g g g g g a g ))]}2()2([(3))()([12)(17{351)(x x g x x g x x g x x g x g x g ∆+∆+∆-∆-∆+∆+∆-∆+∆=∆
(2) 平面多项式圆滑
与剖面圆滑类似，用二维多项式模拟异常，即
用同样的方法，可以导出平面圆滑计算公式，如九点二次圆滑公式如下：
同理，可以导出二十五点二次、四十九点二次、二十五点三次、四十九点三次圆滑计算公式。

现就一具体问题应用matlab 的polyfit 函数来进行试验：
eg1. 用二次多项式拟合数据：
x=[0.1 0.4 0.5 0.7 0.7 0.9];
y=[0.61 0.92 0.99 1.52 1.47 2.03];
并作出数据点和拟合曲线图。

解：用polyfit 求二次多项式的系数，作出曲线图，程序如下：
clear,clf,clc
x=[0.1 0.4 0.5 0.7 0.7 0.9];
y=[0.61 0.92 0.99 1.52 1.47 2.03];
cc=polyfit(x,y,2) %求出A 与B 的系数
二次项的系数为：
xx=x(1):0.1:x(length(x));
yy=polyval(cc,xx);
plot(xx,yy,'--')
hold on
plot(x,y,'x')
axis([0,1,0,3])
xlabel('x')
ylabel('y')
作图如下：
25423210),(y a xy a x a y a x a a y x g +++++=∆)]},(),((),(),([()],(),(),(),([2),(5{91),(y y x x g y y x x g y y x x g y y x x g y y x g y y x g y x x g y x x g y x g y x g ∆+∆-∆+∆+∆+∆+∆-∆+∆+∆-∆-∆-∆+∆+∆-∆+∆+∆+∆-∆+∆=∆919291929592
919291----
eg2.给定一组数据点如下：
t=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 10.00 10.20 10.32 10.42 10.50 10.55 10.58 10.60]*10.^(-3);
求其最小二乘拟合曲线。

解：作出散点图，观察散点图易知变量t 、y 不成线性关系，但是近似成y=βe t
α关系，或成t
b a y +=1关系。

对于近似成y=βe t
α关系，令：
Y=log(y), B=log(β), T=1/t, A=α.
则有：Y=AT+B
程序：
clear,clf,clc
t=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];
y=[4.00 6.40 8.00 8.80 9.22 9.50 9.70 9.86 10.00 10.20 10.32 10.42 10.50 10.55 10.58 10.60]*10.^(-3);
plot(t,y,'.')
C=polyfit(1./t,log(y),1) %求出A 与B 的系数
t1=1:0.1:16;
y1=exp(-4.4807)*exp(-1.0567*1./t1);
hold on
plot(t1,y1,'--') %绘出y与t的近似函数关系曲线图
grid on
作图如下：
可以看到拟合曲线效果非常好。

从上面第二个例子还可以看出，在重力异常资料处理方法中，有些很接近指数形式的拟合曲线可以先化成指数形式，然后再化成多项式，利用matlab的polyfit函数求出系数，最后画出拟合曲线，拟合效果会更佳。