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

现代测量数据处理作业20基于matlab的开采沉陷数据处理姓名：戴超本文是基于Matlab的开采沉降数据处理，对数据的处理方法采用线性回归和多项式拟合两种方法，主要是利用最小二乘原理进行线性最小二乘拟合，通过曲线拟合的方法反映观测量之间的规律，建立拟合函数，从而确定最佳估计参数。Matlab进行曲线拟合主要有两种方法:回归法拟合和多项式拟合，下文将结合开采沉降实例来比较这两种方法各自的优点。1线性回归模型在Matlab统计工具箱中使用命令regress()实现多元线性回归，调用格式为：b=regress(y，x)或[b，bint，r，rint，stats]=regess(y，x，alpha)其中：1）y表示一个1n的因变量数据矩阵。2）x是pn矩阵，自变量x和一列具有相同行数，值是1的矩阵的组合。如：对含常数项的一元回归模型，可将x变为2n矩阵，其中第一列全为1。3）alpha为显著性水平(缺省时设定为0.05)，输出向量b回归细数估计值（并且第一值表示常数，第二个值表示回归系数）。4）bint为b的置信区间。5）r、rint为残差及其置信区间。6）stats是用于检验回归模型的统计量，有四个数值，第一个是2r，其中R是相关系数；第二个是F统计量值；第三个是与统计量F对应的概率P，当alphap时拒绝0H，回归模型成立；第四个2S是误差方差估计值。说明：相关系数2r越接近1，说明回归方程越显著；1,1pnpFFalpha时拒绝0H，F越大，说明回归方程越显著；与F对应的概率alphap时拒绝0H，回归模型成立。2多项式拟合模型在Matlab统计工具箱中使用命令polyfit(x,y,n)实现多项式拟合，调用格式为：[p,S,mu]=polyfit(x,y,n)现代测量数据处理作业21或[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构建范德蒙矩阵V，然后进行QR分解，得到的上三角矩阵。8）df：自由度，df=length(y)-(n+1)。df0时，为超定方程组的求解，即拟合点数比未知数(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序号11121314151617x2.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,0.77,0.11,0.80,0.45,0.46,0.95,0.22,0.45,0.66,0.59,0.80,0.现代测量数据处理作业2287];X=[ones(length(y),1),x'];Y=y';[b,bint,r,rint,stats]=regress(Y,X)Matlab输出结果：b=0.10780.4220bint=-0.00580.22130.34080.5032r=[0.0516-0.0264-0.10350.0303-0.0203-0.0129-0.15390.10990.0046-0.0065-0.1072-0.09880.04260.1091-0.00310.08460.0997]'rint=[-0.12360.2268；-0.19510.1422；-0.25420.0472；-0.14610.2068；-0.19620.1556；-0.18870.1628；-0.2876-0.0202；-0.05640.2761；-0.16820.1775；-0.18020.1673；-0.25360.0391；-0.25360.0561；-0.12650.2117；-0.05610.2744；-0.18040.1743；-0.08610.2553；-0.06720.2666]stats=[0.8910122.64520.00000.0069]Matlab运算的相关结果见表2。表2线性回归相关结果回归系数回归系数估计值回归系数置信区间0ˆb0.1077-0.00580.22131ˆb0.41230.34080.50320.89102r22.64521F0000.0P0069.02S验证模型的有效性：（1）由结果知，残差r均未超过rint中各残差相应的置信区间，运算结果有效。（2）0ˆb的置信区间含零点，说明结果无效。（3）0.90412r，说明显著性较好。（4）0ˆb、1ˆb置信区间不大，说明有效性较好。在Matlab中，输入rcoplot(r,rint)，进行残差图分析。得到残差图如图1所示：现代测量数据处理作业23图1可以看出第七个点的置信区间不包含零点，认为这个数据异常，将其取出再次计算，得到最终结果：b=[0.23930.3406]'bint=[0.13140.3473；0.26760.4135]r=[0.0186-0.0001-0.05840.0168-0.0241-0.01430.0907-0.0618-0.0688-0.0556-0.03120.0630-0.04100.07020.0959]'rint=[-0.10850.1456；-0.12200.1219；-0.16720.0504；-0.11120.1448；-0.15100.1029；-0.14160.1131；-0.02470.2060；-0.17530.0517；-0.18220.0445；-0.16160.0503；-0.14650.0842；-0.05620.1823；-0.16530.0833；-0.05050.1910；-0.01710.2090]stats=[0.8867101.76070.00000.0035]改进后的回归模型的系数、系数置信区间与统计量，如表3。表3线性回归相关结果回归系数回归系数估计值回归系数置信区间0ˆb0.23930.13140.34731ˆb0.34060.26760.41358867.02r7607.011F0000.0P0035.02S验证模型的有效性：（1）残差r均未超过rint中各残差相应的置信区间，运算结果有效。（2）0ˆb、1ˆb的置信区间均不含零点，说明结果有效。（3）8867.02r，较为接近1，说明显著性较好。（4）0ˆb、1ˆb置信区间不大，说明有效性有效。（5）101.7607F，数值较大，说明显著性较好。现代测量数据处理作业24（6）概率0000.0P，05.0P，拒绝0H，回归模型成立。所以回归模型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）多项式拟合模型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]';现代测量数据处理作业25y=[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.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的CommandWindow窗口中，得到结果：p1=0.42200.1078s1=R:[2x2double]df:15normr:0.3217根据得到的拟合方程系P数列出拟合方程：1078.04220.0xy为了更直观的得到线性多项式一次拟合图像，输入下列命令代码：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.45,0.66,0.59,0.80,0.87]';[p2,s2]=polyfit(x,y,2)现代测量数据处理作业26将上述指令代码输入到Matlab的CommandWindow窗口中，得到结果：p1=-0.21630.9899-0.2044s1=R:[3x3double]df:14normr:0.1649根据得到的拟合方程系P数列出拟合方程：2044.09899.0-0.2163xy为了更直观的得到线性多项式二次拟合图像，输入下列命令代码：x1=0:0.01:2.5;y2=polyval(p2,x1);plot(x,y,'*r',x1,y2)得到多项式二次拟合图像，如图5。图5多项式二次拟合对该组观测值进行三次多项式拟合，拟合所需的命令代码：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的CommandWindow窗口中，得到结果：p1=0.0805-0.53711.3658-0.3266s1=R:[4x4doub