基于matlab的开采沉陷数据处理

1、基于matlab 的开采沉陷数据处理 姓名：戴超本文是基于Matlab的开采沉降数据处理，对数据的处理方法采用线性回归和多项式拟合两种方法，主要是利用最小二乘原理进行线性最小二乘拟合，通过曲线拟合的方法反映观测量之间的规律，建立拟合函数，从而确定最佳估计参数。Matlab 进行曲线拟合主要有两种方法:回归法拟合和多项式拟合，下文将结合开采沉降实例来比较这两种方法各自的优点。1线性回归模型在Matlab统计工具箱中使用命令regress()实现多元线性回归，调用格式为：b=regress(y，x)或b，bint，r，rint，stats = regess(y，x，alpha)其中：1）表示一个的

2、因变量数据矩阵。2）是矩阵，自变量和一列具有相同行数，值是1的矩阵的组合。 如：对含常数项的一元回归模型，可将变为矩阵，其中第一列全为1。3）为显著性水平(缺省时设定为0.05)，输出向量回归细数估计值（并且第一值表示常数，第二个值表示回归系数）。4）为的置信区间。5）、为残差及其置信区间。6）是用于检验回归模型的统计量，有四个数值，第一个是，其中是相关系数；第二个是统计量值；第三个是与统计量对应的概率，当时拒绝，回归模型成立；第四个是误差方差估计值。说明：相关系数越接近1，说明回归方程越显著； 时拒绝，越大，说明回归方程越显著；与对应的概率时拒绝，回归模型成立。2多项式拟合模型 在Matla

3、b统计工具箱中使用命令polyfit(x,y,n)实现多项式拟合，调用格式为：p,S,mu=polyfit(x,y,n)或p,S=polyfit(x,y,n)或p=polyfit(x,y,n)其中：1）x，y为已知数据点向量，分别表示横、纵坐标。 2）n为拟合多项式的次数。3）向量p是返回n次拟合多项式系数，从高次到低次。4）矩阵S用于生成预测值的误差估计。5）mu=mean(x); std(x)，mean(x)求x每一列的均值，std(x)求x的标准差。6）S是一个结构体数组(struct)，用来估计预测误差，包含了R，df和normr。7）R：polyfit函数中，先根据输入的x构建范德蒙

4、矩阵V，然后进行QR分解，得到的上三角矩阵。8）df：自由度， df=length(y)-(n+1)。df>0时，为超定方程组的求解，即拟合点数比未知数(p(1)p(n+1)多。9）normr：标准偏差、残差范数，normr=norm(y-V*p)，此处的p为求解之后的数值。3实例分析某矿区测得17个观测站成果如表1，除表所列数据外，这些观测站的其他地质采矿条件均相同。表1 某矿区各观测站实测下沉值序号12345678910x1.211.942.171.451.571.600.371.380.800.85y0.670.900.920.750.750.770.110.800.450.46序

5、号11121314151617 x2.250.50.711.051.151.441.57y0.950.220.450.660.590.800.87Matlab进行曲线拟合主要有两种方法：回归法拟合和多项式拟合,下面分别对两种方法进行分析，并比较这两种方法各自的优点。1）线性回归模型以这17期观测数据为基准，进行Matab线性回归计算，Matlab进行线性回归程序代码为：x=1.21,1.94,2.17,1.45,1.57,1.60,0.37,1.38,0.80,0.85,2.25,0.50,0.71,1.05,1.15,1.44,1.57;y=0.67,0.90,0.92,0.75,0.75,

6、0.77,0.11,0.80,0.45,0.46,0.95,0.22,0.45,0.66,0.59,0.80,0.87;X=ones(length(y),1),x'Y=y' b , bint , r , rint , stats =regress ( Y , X)Matlab输出结果：b = 0.1078 0.4220bint = -0.0058 0.2213 0.3408 0.5032r =0.0516 -0.0264 -0.1035 0.0303 -0.0203 -0.0129 -0.1539 0.1099 0.0046 -0.0065 -0.1072 -0.0988 0.

7、0426 0.1091 -0.0031 0.0846 0.0997'rint =-0.1236 0.2268；-0.1951 0.1422；-0.2542 0.0472；-0.1461 0.2068；-0.1962 0.1556；-0.1887 0.1628；-0.2876 -0.0202；-0.0564 0.2761；-0.1682 0.1775；-0.1802 0.1673；-0.2536 0.0391；-0.2536 0.0561；-0.1265 0.2117；-0.0561 0.2744；-0.1804 0.1743；-0.0861 0.2553；-0.0672 0.2666s

8、tats =0.8910 122.6452 0.0000 0.0069Matlab运算的相关结果见表2。表2 线性回归相关结果回归系数回归系数估计值回归系数置信区间 -0.0058 0.22130.3408 0.5032 验证模型的有效性：（1）由结果知，残差r均未超过rint中各残差相应的置信区间，运算结果有效。（2）的置信区间含零点，说明结果无效。（3），说明显著性较好。（4）、置信区间不大，说明有效性较好。在Matlab中，输入rcoplot(r,rint)，进行残差图分析。得到残差图如图1所示：图1 可以看出第七个点的置信区间不包含零点，认为这个数据异常，将其取出再次计算，得到最终结果

9、：b = 0.2393 0.3406 'bint =0.1314 0.3473；0.2676 0.4135r =0.0186 -0.0001 -0.0584 0.0168 -0.0241 -0.0143 0.0907 -0.0618 -0.0688 -0.0556 -0.0312 0.0630 -0.0410 0.0702 0.0959 'rint =-0.1085 0.1456；-0.1220 0.1219；-0.1672 0.0504；-0.1112 0.1448；-0.1510 0.1029；-0.1416 0.1131；-0.0247 0.2060；-0.1753 0.

10、0517；-0.1822 0.0445；-0.1616 0.0503；-0.1465 0.0842；-0.0562 0.1823；-0.1653 0.0833；-0.0505 0.1910；-0.0171 0.2090stats = 0.8867 101.7607 0.0000 0.0035改进后的回归模型的系数、系数置信区间与统计量，如表3。表3 线性回归相关结果回归系数回归系数估计值回归系数置信区间0.1314 0.34730.2676 0.4135 验证模型的有效性：（1）残差r均未超过rint中各残差相应的置信区间，运算结果有效。（2）、的置信区间均不含零点，说明结果有效。（3），较为

11、接近1，说明显著性较好。（4）、置信区间不大，说明有效性有效。（5），数值较大，说明显著性较好。（6）概率，拒绝，回归模型成立。所以回归模型Y=0.2393+0.3406x成立。此时，在Matlab中，再次输入rcoplot(r,rint)，进行残差图分析。得到残差图如图2所示：图2可以看出各点的置信区间均包含零点，认为所有数据正常。在Matlab命令窗口中输入：z=b(1)+b(2)*xplot(x,Y, 'k+',x,z, 'r')可得到离散点与线性回归图形，如图3。图3由图3可看出，除个别点偏离曲线较远，多元线性回归曲线基本可以满足离散点分布。2）多项式拟

12、合模型Matlab进行多项式拟合，可直接调用polyfit函数。首先，将表2中的观测数据写成：x=1.21,1.94,2.17,1.45,1.57,1.60,0.37,1.38,0.80,0.85,2.25,0.50,0.71,1.05,1.15,1.44,1.57'y=0.67,0.90,0.92,0.75,0.75,0.77,0.11,0.80, 0.45,0.46,0.95,0.22,0.45,0.66,0.59,0.80,0.87'然后，参考Matlab中polyfit调用格式，以多项式一次拟合为例，写出进行多项式拟合所需的命令代码：x=1.21,1.94,2.17,1

13、.45,1.57,1.60,0.37,1.38,0.80,0.85,2.25,0.50,0.71,1.05,1.15,1.44,1.57'y=0.67,0.90,0.92,0.75,0.75,0.77,0.11,0.80, 0.45,0.46,0.95,0.22,0.45,0.66,0.59,0.80,0.87' p1,s1= polyfit(x,y,1)将上述指令代码输入到Matlab的Command Window窗口中，得到结果：p1 = 0.4220 0.1078s1 = R: 2x2 double df: 15 normr: 0.3217根据得到的拟合方程系P数列出拟合

14、方程：为了更直观的得到线性多项式一次拟合图像，输入下列命令代码：x1=0:0.01:2.5;y1=polyval(p1,x1);plot(x,y,'*r',x1,y1)得到多项式一次拟合图像，如图4。图4 多项式一次拟合对该组观测值进行二次多项式拟合，拟合所需的命令代码：x=1.21,1.94,2.17,1.45,1.57,1.60,0.37,1.38,0.80,0.85,2.25,0.50,0.71,1.05,1.15,1.44,1.57'y=0.67,0.90,0.92,0.75,0.75,0.77,0.11,0.80, 0.45,0.46,0.95,0.22,0.

15、45,0.66,0.59,0.80,0.87' p2,s2= polyfit(x,y,2)将上述指令代码输入到Matlab的Command Window窗口中，得到结果：p1 = -0.2163 0.9899 -0.2044s1 = R: 3x3 double df: 14 normr: 0.1649根据得到的拟合方程系P数列出拟合方程：为了更直观的得到线性多项式二次拟合图像，输入下列命令代码：x1=0:0.01:2.5;y2=polyval(p2,x1);plot(x,y,'*r',x1,y2)得到多项式二次拟合图像，如图5。图5 多项式二次拟合对该组观测值进行三次多

16、项式拟合，拟合所需的命令代码：x=1.21,1.94,2.17,1.45,1.57,1.60,0.37,1.38,0.80,0.85,2.25,0.50,0.71,1.05,1.15,1.44,1.57'y=0.67,0.90,0.92,0.75,0.75,0.77,0.11,0.80, 0.45,0.46,0.95,0.22,0.45,0.66,0.59,0.80,0.87' p3,s3= polyfit(x,y,3)将上述指令代码输入到Matlab的Command Window窗口中，得到结果：p1 = 0.0805 -0.5371 1.3658 -0.3266s1 = R

17、: 4x4 double df: 13 normr: 0.1577根据得到的拟合方程系P数列出拟合方程：为了更直观的得到线性多项式三次拟合图像，输入下列命令代码：x1=0:0.01:2.5;y3=polyval(p3,x1);plot(x,y,'*r',x1,y3)得到多项式三次拟合图像，如图6。图6 多项式三次拟合然后依据下面命令代码，将三条拟合曲线绘于同一图上：x1=0:0.1:2.5; % 0到2.5的100个点(每隔0.01画一个点)的x坐标 y1=polyval(p1,x1);y2=polyval(p2,x1);y3=polyval(p3,x1);plot(x,y,&

18、#39; *r',x1,y1, ':r' ,x1,y2, ' -r',x1,y3, ' -.r' ) % 同时画三个函数通过Matlab计算，得到多项式拟合图像如图7。图7 多项式拟合图像由图7可见，多项式一次、二次、三次拟合越来越接近离散点的真实分布。综合这三次拟合的结果，可以看出，由于未知系数的增多，多项式一次、二次、三次拟合的自由度依次减小，标准偏差越来越小，拟合的准确度越来越高，详情见表2。表2一次拟合二次拟合三次拟合自由度151413标准偏差0.32170.16490.15773）分析对比Matlab进行线性回归，得到的初始的拟

19、合方程： (1)经过异常点剔除得到拟合方程： (2)Matlab进行多项式拟合，得到的相应拟合方程：一次拟合方程： (3)二次拟合方程： (4)三次拟合方程： (5)由方程(1)、(3)知多项式拟合一次曲线和线性回归曲线是一样的，但式(2)经过异常点剔除，所以其方程的有效性比(1)、(3)要高。由图7知，多项式拟合中，三次拟合最为接近离散点，且由于多项式拟合一次曲线和剔除异常点后的线性回归曲线类似，所以这5个方程中，(5)式最为接近离散点，显著性最强。由于三次多项式拟合最接近原始离散点，对三次多项式拟合利用polyconf函数拟合结果进行评价，利用polyconf函数求polyfit所得的回归

20、多项式在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y±DELTA，alpha缺省时为0.05。 polyconf()函数的调用格式为：Y=polyconf(p,x,s)Y,DELTA=polyconf(p,x,s,alpha)说明：Y=polyconf(p,x,s)使用polyfit函数的选项输出s给出Y的95%置信区间YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数，1-alpha为置信度。Matlab命令代码为：alpha=0.05;Y,DELTA=polyconf(p3,x,s3,alpha)Matlab运行结果如下：Y=0.6822 0.8893 0.9304 0.7699 0.8053 0.8134 0.1093 0.7469 0.4635 0.4957 0.9441 0.2321 0.4012 0.6085 0.6562 0.7667 0.805