高教社杯全国大学生数学建模竞赛C题论文

2012高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。我们参赛选择的题号是（从A/B/C/D中选择一项填写）：C 我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：临沧师范高等专科学校参赛队员(打印并签名)：1.饶兵2.田兴平3.熊生指导教师或指导教师组负责人(打印并签名)：穆志勇日期：2013年9赛区评阅编号（由赛区组委会评阅前进行编号）：2012高教社杯全国大学生数学建模竞赛编号专用页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：2013高教社杯全国大学生数学建模竞赛C题论文摘要：文物是人类的重要历史文化遗产，它不仅反映过去，还可以借鉴未来，文物保护是一项重要的工作，不仅要有高度的文物保护意识，还应该有强有力的措施和先进的科学技术手段。本题给出了某古塔在1986年7月、1996年8月、2009年3月和2011年3月进行的4次观测数据，数据包含了古塔每一层以及塔尖上的采样点的三维坐标值（x,y,z）。要求确定古塔各层的中心坐标，分析古塔的倾斜和变形情况，最后再分析古塔的变形趋势，以便于加强防护措施，保护古塔。本文根据题目附件给出的三维坐标数据，做出了各次测量数据采样点的三维散点图，画出各条相连边的连线，求出该空间多边形的中心点坐标值。算出了各楼层的中心坐标值之后，可以画出4次测量的各楼层中心连线图，再连接上塔的顶点，就刻画出古塔的外形。根据该连线图，我们可以分析古塔的倾斜、扭曲、弯曲等变形情况。再将4次测量画出的中心连线图作比较，就可以知道古塔的变形趋势。关键词：文物保护、古塔、垂径定理、中心、线性回归一、问题重述古塔的变形由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。请你们根据附件1提供的4次观测数据，讨论以下问题：1.给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。2.分析该塔倾斜、弯曲、扭曲等变形情况。3.分析该塔的变形趋势。二、问题分析根据数据的散点图初步判断古塔是至少具有八个侧面的多边形棱锥体（或者棱台体），也有可能是圆锥体（或者圆台体）。每层楼具有八个采样点，将八个点投影在水平面（x,y平面）上，采样点的在水平面的坐标值为（x,y），将相连的两个采样点用直线连接，连线就是外接圆的玄，做出八条边（亦及八条外接圆的玄）的中垂线，中垂线的交点（或者多个交点的平均值）就是圆心的水平坐标（xo,yo）,而（外接圆的圆心）是该楼层的中心在水平面上的投影，也就说楼层中心的水平坐标值也是（xo,yo）。再根据该楼层八个采样点的垂直坐标（z轴）的平均值，来近似的估算中心点的垂直坐标zo的值。最后就得出了该楼层的中心三维坐标值（xo,yo,zo）。三、模型的假设3.1、假设古塔每一层的水平横截都是中心对称图形。（保证中心点坐标的平均值居中）3.2、古塔每层的重心就是它的几何中心。四、定义与符号说明直线AB的斜率直线AB的中垂线的斜率AB与BC的中垂线的交点xy平面内直线和y轴的夹角θ空间直线和xy水平面的夹角五、模型建立与求解5.1画散点图，对塔的形体做初步判断在第u次测量中，根据数据附件中给出的每次测量的数据，画出了各次采样点的散点图，并将古塔同一楼层的点用线连起来，例如2011年的测量数据采样点连线图如下：根据图像初步判断塔体是至少具有八个侧面的多边形棱锥体（或者棱台体），也有可能是圆锥体（或者圆台体）。并且假设其每个水平截面都是中心对称图形。5.2求出古塔每层楼的中心坐标值并画出连线图在第u次测量过程中，第j楼的八个观测点的坐标值分别为：（xj,yj,zj）,楼层序数j=1,2,...8。提取八个采样点在水平面（x,y）上的投影坐标值：（xj,yj）,j=1,2,...8，连接相连的两个点，画得一个八边形，如图所示：做出第一条变AB的中垂线，中垂线与AB交于PAB点，坐标值为（x12,y12）,根据中点坐标公式可得：，。直线AB的斜率为：，那么直线AB的中垂线的斜率为：根据直线的点斜式方程，已知直线上的一个点，又已知直线的斜率，则AB的中垂线的方程可以确定。方程为：，带入的值得：。同理可算得BC的中点PCB的坐标值和BC边中垂线的斜率。PCB（x23,y23）中，,。也可以确定了BC的中垂线的方程。AB和BC直线的中垂线交于点，已知两直线的方程，联立两条中垂线的方程，可以解得它们的交点坐标值：。以此类推，可以算得所有边（八条边）的中垂线的方程，然后出它们两两相交的交点坐标。如果八边形的外接圆存在，那么所有中垂线的交点都应该重合于圆心，但是由于测量误差和形变，这些中垂线的交点并不重合，因此，取所有交点坐标的平均值，作为外接圆的圆心值，也就是八边形的中心坐标值。实现的方法是用组合，先算出第1条中垂线（即AB边的中垂线）与其它7条边的中垂线的交点，再算第2条中垂线（即BC边的中垂线）与第3条到第8条中垂线的交点，依次类推，算出所有中垂线之间的交点坐标。如果是八条边则产生个交点，如果是L条边则产生个交点。总的交点的坐标值为：,其中i,j为边的序号，并且要求另外是等价的，只能出现一次。对于中心点的垂直坐标，用,其中L是边的条数，也是顶点的个数。根据上述方法，算得每层楼的中心点坐标，然后用直线连接，如图所示：5.3数据和分析从图中可以发现，第1、2次（1986年和1996年）测量的数据非常接近，图像几乎重合；第3、4次（2009年和2011年）测量的数据非常接近，图像几乎重合。而在第2次（1996年）测量值和第3次（2009年）测量值之间，存在着较大差异。第1次1986年测量各楼层中心坐标值列表楼层坐标轴12345678910111213楼顶x566.68566.73566.77566.81566.85566.89566.93566.96567567.04567.11567.17567.21567.25y522.7522.66522.63522.59522.56522.53522.5522.47522.44522.42522.39522.36522.32522.24z1.78747.320212.75517.07821.72126.23529.83733.35136.85540.17244.44148.71252.83455.123第2次1996年测量各楼层中心坐标值列表楼层坐标轴12345678910111213楼顶x566.68566.73566.77566.81566.85566.9566.93566.97567.01567.04567.11567.18567.22567.25y522.7522.66522.62522.59522.56522.52522.5522.47522.44522.42522.39522.35522.31522.24z1.7837.314612.75117.07521.71626.2329.83233.34536.84840.16844.43548.70752.8355.12第3次2009年测量各楼层中心坐标值列表楼层坐标轴12345678910111213楼顶x566.76566.81566.84566.88566.91567.02567.06567.12567.19567.24567.28567.32567.37567.34y522.69522.65522.61522.58522.54522.52522.51522.5522.49522.43522.38522.33522.28522.21z1.76457.30912.73217.0721.70926.21129.82533.3436.84440.16144.43348.752.81855.091第4次2011年测量各楼层中心坐标值列表楼层坐标轴12345678910111213楼顶x566.77566.81566.85566.88566.91567.02567.06567.12567.19567.24567.28567.32567.37567.34y522.69522.65522.61522.58522.54522.52522.51522.5522.49522.43522.38522.33522.28522.21z1.76327.290512.72717.05221.70426.20529.81733.33736.82240.14444.42548.68452.81355.0875.3.1静态分析由于第1、2次（1986年和1996年）测量的数据非常接近，所以选取第2次（1996年）测量的数据为例进行分析。根据1996年数据画得古塔各楼层中心连线图，如图红色所示，为了便于观察，相对前一个图而言，将图像旋转了大约180度。在xy水平面（z=0)、xz侧面（y=552.8时的平面）和yz侧面（x=566.5时的平面）的投影都接近直线，因此可以用线性拟合的方法，将数据拟合为空间直线。因此对连线做了空间直线拟合，拟合后的图像如图所示。另外因为楼顶坐标（上图中的五角星）不是某一楼层的中心，因此我们不放在拟合直线的参考范围之内。算出了该拟合直线的空间直线方程为x=Az+C1,y=Bz+C2,据此我们可以算出直线在xy水平面内的投影方程为，根据此方程可以算得古塔中心连线的倾斜方向与x轴和y轴的夹角。与y轴的夹角，=-54.439°；那么与x轴的夹角等于144.439°。为了算得空间直线与xy水平面的夹角，我们在空间直线上任取两点（x1,y1,z1）,(x2,y2,z2),课算得空间直线和它在xy平面内的夹角，可证明这个夹角就是空间直线和xy平面的夹角。设夹角为：θ，则，就得，任取两点带入，可算得89.276°，也就说该直线偏离了垂直（z轴）方向0.72353°。第3、4次（2009年和2011年）测量的数据非常接近，因此取第3次（2009年）测量的数据为例进行分析。根据数据画得如下的图像：从图中可以发现，各个楼层的中心连线图是扭曲的，在xy水平面（z=0)、xz侧面（y=552.8时的平面）和yz侧面（x=566.5时的平面）的投影都是弯曲的，说明各楼层的中心分布不是线性的，不能简单的用直线来代替。但我们还是先用线性拟合直线的算法来算一个平均值，后面再具体分析扭曲的情况。用前面的算法，可以算得xy水平投影与y轴夹角为-60.120°，与x轴夹角为：150.12°根据拟合直线的方程，可算得拟合直线和水平面的夹角为89.180°，也就说该直线偏离了垂直方向（z轴方向）0.81989°。根据中心连线及其投影图，我们将它分为三段来进行计算：1-5楼、5-9楼、9-13楼。分别用直线代替了这三条线段，如图所示：根据拟合直线的方程算得：夹角名称楼层段水平投影和y轴的夹角（度）水平投影和x轴的夹角（度）所在直线和xy水平面的夹角（度）偏离垂直线的角度（与z轴的夹角）1-5层-44.0669134.066989.39860.60135-9层-79.6511169.651188.97191.02809-13层-40.2861130.286189.02820.9717根据图表数据分析，中间段5-9层的角度与上下两层的角度相差较大，偏离垂直线的角度也是最大，倾斜严重，最有可能在此段发生坍塌。5.3.2动态分析。为了比较1996年和2009年的数据，从总体拟合直线上来比较两组数据，列表如下：数据名测量序数水平投影和y轴的夹角（度）水平投影和x轴的夹角（度）所在直线和xy水平面的夹角（度）偏离垂直线的角度（与z轴的夹角）2009年（第3次）-60.120150.1289.1800.819891996年（第2次）-54.439144.4389.2760.72353差值-5.681-5.690-0.0960.09636可以发现总体上它们的倾斜的角度的差值并不明显。将1996年和2009年的数据拟合图的起点重合后，两次拟合连线的对图比如下：将2009年（第3次）分为三段拟合直线的夹角，分别和1996年（第2次）的拟合直线的夹角值进行比较。古塔2009年测量所得角度和1996年测量所得角度分层比较表：角度量分段数水平投影和y轴的夹角之差水平投影和x轴的夹角之差所在直线和xy水平面的夹角之差偏离垂直线（z轴）的角度差1-5层10.372-10.3630.1226-0.12225-9层-25.21225.221-0.30410.30449-13层14.153-14.143-0.24780.2481由表格可以看出，5-9层在1996年至于2009年期间，发生了较大的形变，角度变化最大，致使整个塔身发生了扭曲。六、对模型的评价6.1模型优点简化了塔的模型，确定了每层塔的中心（重心），以连线代替塔体。可以充分展示塔身的倾斜、变形、弯曲和扭曲情况。6.2模型缺点实际的塔体情况可能很复杂，其切面不一定是中心对称图形。只能通过平均值来尽量接近中心值。6.3改进如果能增加测量点的数据个数是最好的，在限制现有数据个数的情况下，可以进一步用数据筛选方法，淘汰偏差较大的数据。七、参考文献[1]李柏年.MATLAB数据分析方法[M].北京:机械工业出版社,2012:12-19.[2]王能超.数值分析简明教程[M].北京：高等教育出版社，2003：36-37.[3]于晓秋，李欣，张宏礼.MATLAB数值实验[M].沈阳：辽宁科学技术出版社，2004：34-40.八、附录附录一，主程序（保存为文件C2013.m）：gcsj1=xlsread('D:\gcsj.xls','sheet1','C4:E111');%从电子表格中读取对应测量序数的数据gcsj2=xlsread('D:\gcsj.xls','sheet1','I4:K111');gcsj3=xlsread('D:\gcsj.xls','sheet1','O4:Q111');gcsj4=xlsread('D:\gcsj.xls','sheet1','U4:W111');symsYoXoZoyoxoXYyo2xo2%画散点连线图foru=1:4figure;fork=0:12x=gcsj((k*8+1):(8*k+8),1);%读取第k+1楼层的坐标数据值y=gcsj((k*8+1):(8*k+8),2);z=gcsj((k*8+1):(8*k+8),3);x=x(~isnan(x));%去除空白数据NaNy=y(~isnan(y));z=z(~isnan(z));r1=size(x);r=r1(1);plot3(x,y,z,'r-*');%画采样坐标点连线图gridonholdonplot3([x(1),x(r)],[y(1),y(r)],[z(1),z(r)],'r-*');%画第1点和最后1个的连线图xlabel('x轴');ylabel('y轴');zlabel('z轴');endend%计算古塔每个楼层的中心坐标值figure;foru=1:4switchu%读取第u次测量数据case1gcsj=gcsj1;case2gcsj=gcsj2;case3gcsj=gcsj3;case4gcsj=gcsj4;endYo=0;Xo=0;fork=0:12xq=gcsj((k*8+1):(8*k+8),1);%读取第k+1楼层坐标数据yq=gcsj((k*8+1):(8*k+8),2);zq=gcsj((k*8+1):(8*k+8),3);Yo=0;Xo=0;ll=0;x=xq(~isnan(xq));%去除空白数据NaNy=yq(~isnan(yq));z=zq(~isnan(zq));r1=size(x);r=r1(1);forj=1:r-2forn=j+1:r-1K1=(x(j)-x(j+1))/(y(j+1)-y(j));%算中垂线斜率K2=(x(j+2)-x(j+1))/(y(j+1)-y(j+2));equation1=sym('Y-(y(j)+y(j+1))/2-K1*(X-(x(j+1)+x(j))/2)=0');%解中垂线交点方程equation2=sym('Y-(y(j+2)+y(j+1))/2-K2*(X-(x(j+1)+x(j+2))/2)=0');[xo,yo]=solve(equation1,equation2,'X','Y');yo2=eval(yo);xo2=eval(xo);Yo=Yo+yo2;%累加中心点坐标值Xo=Xo+xo2;ll=ll+1;%累加中心点个数endK1=(x(r-1)-x(r))/(y(r)-y(r-1));K2=(x(1)-x(r))/(y(r)-y(1));equation1=sym('Y-(y(r-1)+y(r))/2-K1*(X-(x(r)+x(r-1))/2)=0');equation2=sym('Y-(y(1)+y(r))/2-K2*(X-(x(1)+x(r))/2)=0');[xo,yo]=solve(equation1,equation2,'X','Y');yo2=eval(yo);xo2=eval(xo);Yo=Yo+yo2;Xo=Xo+xo2;ll=ll+1;%累加中心点个数endYok(u,k+1)=Yo/ll;%计算中心点坐标平均值Xok(u,k+1)=Xo/ll;Zok(u,k+1)=sum(z)/r;ends=size(gcsj);%开始计算顶点坐标xd=gcsj(105:s(1),1);yd=gcsj(105:s(1),2);zd=gcsj(105:s(1),3);xd2=xd(~isnan(xd));%去除空白数据NaNyd2=yd(~isnan(yd));zd2=zd(~isnan(zd));p1=size(xd2);p=p1(1);xt(u)=sum(xd2)/p;%算得顶点的坐标值（xt,yt,zt）yt(u)=sum(yd2)/p;zt(u)=sum(zd2)/p;huatu(u,Xok,Yok,Zok,xt,yt,zt);%画各楼层中心连线图end%1996年数据分析figure;title('1996年测量数据楼层中心连线分析图');plot3(Xok(2,:),Yok(2,:),Zok(2,:),'r-o')%画中心坐标连线图holdon;gridon;xlabel('x轴');ylabel('y轴');zlabel('z轴');plot3([Xok(2,13),xt(2)],[Yok(2,13),yt(2)],[Zok(2,13),zt(2)],'r-pentagram');%最后和顶点连线plot(Xok(2,:),Yok(2,:),'k-o')%中心坐标连线在xy平面上的投影图plot([Xok(2,13),xt(2)],[Yok(2,13),yt(2)],'k-pentagram');%和顶点的连线在xy平面上的投影plot3(566.5*ones(1,13),Yok(2,:),Zok(2,:),'g-o')%画连线在zy平面上的投影plot3([566.5,566.5],[Yok(2,13),yt(2)],[Zok(2,13),zt(2)],'g-pentagram');%最后和顶点连线plot3(Xok(2,:),522.8*ones(1,13),Zok(2,:),'b-o')%画连线在zx平面上的投影plot3([Xok(2,13),xt(2)],[522.8,522.8],[Zok(2,13),zt(2)],'b-pentagram');%最后和顶点连线%1996年数据分析拟合figure;nihe(Xok(2,:),Yok(2,:),Zok(2,:));%调用拟合函数title('1996年测量数据楼层中心连线分析图');%2009年数据分析打印figure;plot3(Xok(3,:),Yok(3,:),Zok(3,:),'r-o')%画中心坐标连线图holdon;gridon;xlabel('x轴');ylabel('y轴');zlabel('z轴');plot3([Xok(3,13),xt(3)],[Yok(3,13),yt(3)],[Zok(3,13),zt(3)],'r-pentagram');%最后和顶点连线plot(Xok(3,:),Yok(3,:),'k-o')%中心坐标连线在xy平面上的投影图plot([Xok(3,13),xt(3)],[Yok(3,13),yt(3)],'k-pentagram');%和顶点的连线在xy平面上的投影plot3(566.5*ones(1,13),Yok(3,:),Zok(3,:),'g-o')%画连线在zy平面上的投影plot3([566.5,566.5],[Yok(3,13),yt(3)],[Zok(3,13),zt(3)],'g-pentagram');%最后和顶点连线plot3(Xok(3,:),522.8*ones(1,13),Zok(2,:),'b-o')%画连线在zx平面上的投影plot3([Xok(3,13),xt(3)],[522.8,522.8],[Zok(3,13),zt(3)],'b-pentagram');%最后和顶点连线title('2009年测量数据楼层中心连线分析图');figure;nihe(Xok(3,:),Yok(3,:),Zok(3,:));%总体拟合直线，并算出夹角title('2009年测量数据楼层中心连线拟合图');legend('中心连线图','拟合后的直线图');%2009年分段讨论：figure;holdon;gridon;xlabel('x轴');ylabel('y轴');zlabel('z轴');%1-5楼数据拟合nihe(Xok(3,1:5),Yok(3,1:5),Zok(3,1:5));%5-9楼数据拟合nihe(Xok(3,5:9),Yok(3,5:9),Zok(3,5:9));%9-13楼数据拟合nihe(Xok(3,9:13),Yok(3,9:13),Zok(3,9:13));title('2009年测量数据楼层中心连线分段拟合图');Xokc=Xok(3,:)-Xok(2,:);%算两次测量坐标的差值Zokc=Zok(3,:)-Zok(2,:);Yokc=Yok(3,:)-Yok(2,:);%1996年与2009年拟合直线比较画图figure;XOKK=Xok(2,:)+Xok(3,1)-Xok(2,1);%将坐标起点重合以便于比较YOKK=Yok(2,:)+Yok(3,1)-Yok(2,1);nihe2(XOKK,YOKK,Zok(2,:));%1-5楼数据拟合nihe2(Xok(3,1:5),Y