2013年全国大学生高教杯数学建模C题古塔变形论文2013高教社杯全国大学生数学建模竞赛承诺书我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。

如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。

我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）： C 我们的参赛报名号为（如果赛区设置报名号的话）：所属学校（请填写完整的全名）：长春工业大学参赛队员 (打印并签名) ：1. 武太彬2. 贾光辉3. 牛文正指导教师或指导教师组负责人 (打印并签名)：李纯净（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。

以上内容请仔细核对，提交后将不再允许做任何修改。

如填写错误，论文可能被取消评奖资格。

）日期： 2013 年 9 月_16_日赛区评阅编号（由赛区组委会评阅前进行编号）：古塔的变形摘要本文对古塔的变形问题建立数学模型，它实质上是一个空间解析几何问题。

首先建立空间解析几何模型，并利用这个模型对问题1进行求解，然后对模型进行数据处理和图形分析，精确地给出了中心位置坐标；之后对于问题2，在问题1的基础上我们计算出中心坐标的拟合直线，计算出曲率的值；最后对于问题3，运用AR自回归模型给出古塔的变形趋势。

在问题1中，通过分析这四年古塔的每一层中的8 个离散点，运用MATLAB 数学软件建立观测数据模拟图，依据此图作出每一层相应散点的投影平面，经过计算基本近似正八边形，得出正八边形的中心，并运用空间直线拟合模型，求出古塔的一条轴心拟合直线；并且运用空间平面拟合模型，求出每一层的8个散点拟合的平面，那么直线与每一平面相交的点，即为每一层的中心坐标，以所求第一层中心坐标为例，1986年第一层中心坐标为(566.6616,522.6994,1.9738)，1996年第一层中心坐标为(566.662,522.6992,2.0117)，2009年第一层中心坐标为(566.6588,522.7012,1.7033)，2011年第一层中心坐标(566.6587,522.7012,1.7018)，其它（算上塔尖）14层见正文表5。

在问题2中，首先运用MATLAB软件画出每一年的俯视图，即各年的平面图，可以看出这四年的古塔是逐渐倾斜和弯曲的，且在第五层开始发生了扭曲。

建立曲率数学模型，运用软件求解出相应的曲率值18.27。

在问题3中，首先建立带季节项的AR自回归模型，运用问题2中所求的曲率作为自变量，代入到模型中，运用时间序列分析的统计软件SPSS，得出古塔变形的趋势，即随年代的增长，古塔的倾斜和弯曲将更加严重，且在第五层由于空间中心的改变，将更加扭曲。

在本文最后，对模型的优缺点及改进之处进行分析。

关键词：空间解析几何最小二乘法拟合曲率AR预测模型MATLAB7.0软件1问题的重述与分析由于许多外在原因，古塔会产生各种变形，倾斜，弯曲，扭曲等。

为了了解各种变形量，测绘公司先后于1986年7月，1996年8月，2009年3月和2011年3月对该塔进行了四次观测。

讨论3个问题，问题1，给出确定古塔各层中心位置的通用方法，列表给出各次测量的古塔各层中心坐标。

本文首先对古塔的变形情况进行分析，可以获取位置信息，且只有四次的观测数据信息。

建立适当的坐标系，先研究一层塔的八个点的投影所得的八个点的坐标点，然后再确定各中心的坐标，从而找到各层的中心点，可以拟合成一条直线。

并且运用空间拟合平面模型，求出每一层的8个散点拟合的平面，那么直线与每一平面相交的点，即为每一层的中心坐标[1,2]。

问题2，分析该塔倾斜，弯曲，扭曲等变形情况。

首先运用MATLAB 软件画出每一年的俯视图，即各年的平面图，可以看出四年的古塔是倾斜还是弯曲或是扭曲，建立相应的数学模型。

问题3，分析该塔的变形趋势。

首先建立带季节项的AR 自回归模型，运用问题2中所求的结果作为自变量，代入到模型中，得出古塔变形的趋势。

2 基本假设1. 观测的数据准确；2. 以古塔的地面为水平面；3. 假设人在观测时是围绕塔的中心进行八个不同的方面进行大量的测量；4. 两点的距离作到小数点后百分位即为等长，后面忽略；5. 不考虑异常点；3 符号说明0:x 已知每个点的横坐标； 0:y 已知每个点的纵坐标； 0:z 已知每个点的竖坐标；,,,:a b c d 直线方程的待估系数；,,:m n p 空间直线的对应,,x y z 的方向向量；,a b θθ：分别是边12a a -和边12b b -由其原始位置旋转的角度；W T : 变化率；4 模型的建立与求解4.1 问题1的求解首先根据已知每一年数据中的每一层的8个散点坐标(,,)i j k x y z 运用MATLAB 软件，得到如下每年观测数据的模拟图。

5625645665685705205251986年观测数据模拟图图1 1986年观测数据模拟图5625645665685705205251996年数据观测数据模拟图图2 1996年观测数据模拟图5625645665685705725185205225245262009年观测数据模拟图图3 2009年观测数据模拟图5625645665685705725185205225245262011年观测数据模拟图图4 2011年观测数据模拟图由图1到图4的模拟图可以直观地看到，四次测量古塔的显著变化，相同坐标上塔的上部越来越小，把图2和图3叠加到一个图。

5655705205255300204060古塔扭曲变化对比图图5 古塔扭曲变化对比图图5中，红线代表1996年的塔图，黄线代表2009年的塔图。

两个塔图叠加在一起对比更加显明，而在第五层出现扭曲，从第五层往上越来越小，到2011年成了一个点的塔尖。

经过以上的图形分析，本文以空间解析几何为背景，运用最小二乘拟合理论给出空间直线与平面的拟合。

下面建立空间内曲线的最小二乘问题，以空间直线为例研究曲线拟合的方法。

已知空间直线的标准方程000x x y y z z X Y Z ---== 整理得直线射影式方程0000z-z +x =az+b y z-z +y =cz+d X x ZY Z⎧=⎪⎪⎨⎪=⎪⎩（）（），其中a X Z =;00b x -z X Z =c Y Z =;00d y -z Y Z =； 这样直线可以看作是用这 2 个方程表示的平面相交的直线，所以可以分别对 2 个方程进行数据拟合。

设i i ˆˆx=az +b 表示按拟合方程x=az+b 求得的近似值。

一般地，它不同于实测值i x 两者之差m 2i i i=1Q [x -(az )]X b =+∑，同理可得m2i i i=1Q [-(cz )]X y d =+∑当Q 取最小值时,,,a b c d 的值即为方程的系数，即满足下列方程时Q 值最小Q Q Q Q 0,0,0,0y y X Xa b c d∂∂∂∂====∂∂∂∂ 有111122111111,m m m m i i i i i i i i m m m m m mi i i i i i i i i i i i i i bN a z x dN c z y b z a z x z d z c z y z ==========⎧⎧+=+=⎪⎪⎪⎪⎨⎨⎪⎪+=+=⎪⎪⎩⎩∑∑∑∑∑∑∑∑∑∑ 令 12111m z z z F ⎡⎤=⎢⎥⎣⎦(1) 方程组(1)可写成'',FF A FX FF B FY ==其中''''[3.4]11[];[x ...x ];[];[...].m m A a b X B c d Y y y ==== 根据m 组数据点解方程组就可以求得,,,a b c d 的值。

566.8567567.2522.4522.6522.81996年古塔各层中心点拟合直线分析图图6 1996年古塔各层中心点拟合直线分析图566.8567567.2522.4522.6522.82009年古塔各层中心点拟合函数分析图图7 2009年古塔各层中心点拟合函数分析图从图6，图7中可以得出：由1996年古塔各层中心点拟合直线与2009年古塔各层中心点拟合直线，方程为，由一次方程的系数k 看出1996年为0.0105,2009年为0.0112，可以得出古塔的中心发生了一次倾斜，且2009年比1986年倾斜角度变大。

运用附件中MATLAB 的程序经过调试，得到如下结果：1986年0.0104566.64110.0066522.7125x z y z =+⎧⎨=-+⎩1996年0.0105566.64110.0067522.7124x z y z =+⎧⎨=-+⎩2009年0.0112566.66430.0084522.7436x z y z =+⎧⎨=-+⎩2011年0.0112566.66420.0084522.7436x z y z =+⎧⎨=-+⎩接下来运用最小二乘法拟合求空间平面方程： 设空间平面方程为z ax by c =++，其中,,a b c 为待估参数。

设古塔的每层8个点，设点坐标{(,,),1,2,,8}i i i x y z i =，考虑到数据在,,x y z 三个方向均存在误差，得矩阵形式：()L A E x L E ∧+=+即1122,3881122,3881111,11n x y x y n x y v v x y v v x y A E v x y v v ∧⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦。

12,1812,13,1,,n z z L n n z v z a v z x b L E c z v ⎡⎤⎡⎤⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎣⎦通常采用矩阵奇异值分解解算待定参数的整体最小二乘解。

12118(1)[]0T T m m A L U U V U V +-+∑⎡⎤⎡⎤==∑⎢⎥⎢⎥⎣⎦⎣⎦其中1112111112121222111[],,[]0Tmm m m V V V U U U V V ∑⎡⎤==∑=⎢⎥∑⎣⎦则参数的整体最小二乘估计为：11222TLS x V V -=-残差矩阵为：1221222[][]TT L E E U V V ∧=∑[5,6]。