数学建模项目实践

杨晶晶　欧冰

摘 要： 为及时了解古塔变形的情况，制定必要的保护措施，作者首先通过四年的观测数据建立数学模型，拟合出古塔各个塔层的中心坐标、中心轴线，再依据中心轴线的变化分析古塔的倾斜、弯曲、扭曲等变化情况和变化趋势。

关键词： 数学模型 数据拟合 斜率 曲率

一、问题的提出与分析

由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，制定必要的保护措施。某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。文物部门希望通过这4次观测数据确定古塔各层中心位置和古塔各层中心坐标，同时分析该塔倾斜、弯曲、扭曲等变形情况。

为解决此问题，我们需要建立相关数学模型，通过模型的解答与分析判断古塔的心位置，古塔各层中心坐标，以及该塔倾斜、弯曲、扭曲等变形情况。为方便建立数学模型，我们假设古塔各个塔层和塔尖的高度都相同；假设古塔每层有8个观测点且8个点在各层内离层底相对高度相同；假设每个年份对每一层的观测点都固定；假设每层观测点在古塔的中心轴周围。通过把每一层的相关观测数据进行算术平均化处理，通过对每一层的不同年份的数据进行算术平均，得出13个塔层4个年份的三维坐标，最后运用数据拟合的方法得到古塔各个塔层和塔尖的中心坐标方程及中心轴线方程。然后分析塔的倾斜情况。我们对已知的多点观测数据进行算术平均处理，把处理后的数据进行多项式拟合后得出拟合方程组，构建空间直角坐标系，通过勾股定理并且进行数角转换便可得出倾斜角度。最后分析塔的弯曲情况。通过对已知数据进行2次多项式拟合，选取低、中、高三个层次利用曲率比较分析古塔4个年份的弯曲程度。曲率越大，表示古塔的弯曲程度越大。

二、古塔中心坐标、中心轴线的求解

我们把各测量数据进行算术平均化处理，通过数据拟合得出各年份的各层拟合中心位置坐标，以第一层为例，1986年的中心点坐标为（566.657142，522.70516，1.787375），1996年的中心点坐标为（566.66，493.1，1.78），2009年的中心点坐标为（566.67，522.74，1.76），2011年的中心点坐标为（566.67，522.74，1.76）。其余各层可用类似的方法得到中心点坐标。

为了减少误差，使结果更精确，我们对得出的4个年份的中心点坐标再进行一次算术平均处理和数据拟合处理，即可得出第一层的中心点坐标为M1（566.67，522.71，2.13），第二层的中心点坐标为M2（566.72，522.68，6.38）。其余各层的中心点坐标可类似得出。通过各层的中心点坐标可以直观看出各个塔层的水平方向的中心坐标没有发生较大偏差，说明塔层越高偏移相对性大且呈现出倾斜的趋势。

在上述模型的求解过程中，我们使用的是算术平均，经过实践，我们发现运用加权平均的处理方法可以减少塔中心坐标的误差，使模型更加贴近真实情况。我们对年份进行一个加权平均处理，随着年份的增加，所占权重越小；我们把权重设为1986年占40%，1996占30%，2009年占20%，2011年占10%。经过类似计算处理可得到各塔层中心坐标为：

M1（566.7，522.7，2.1），M2（566.7，522.7，6.4），M3（566.8，522.7，10.6），

M4（566.8，522.6，14.9），

M5（566.9，522.6，19.1），

M6（566.9，522.6，23.4），

M7（566.9，522.5，27.7），

M8（566.99，522.5，31.9），

M9（567.1，522.5，36.2），

M10（567.1，522.5，40.4），

M11（567.1，522.4，44.7），

M12（567.2，522.4，48.9），

M13（567.2，522.4，53.2），M14（567.2，522.3，55.1）。

如上图所示：从图形上不但可以直观地看出古塔空间中轴线、各个塔层段的中心点，更可以通过塔层中心点的连线判断出古塔的变形趋势。

三、古塔的变形

我们利用上述所求得的古塔中心轴线，建立空间三维坐标系，再通过数据拟合出四条空间直线，即可求得斜率。

对于古塔弯曲的分析要从古塔的曲率入手。首先我们用Matlab软件分别对1986年7月、1996年8月、2009年3月和2011年3月这四次观测的数据进行拟合，拟合成四个方程组。我们以1986年的多项式为例，进行空间曲线的曲率计算。

1986年：

y■=-0.574605251x■■+650.946992376x■-183835.161349422z■=-34.4235380x■■+39129.9111344x■-11119834.7987881

首先讨论其在P点（x1，y1，15）的曲率，将z■=15代入方程，解得x■=566.791589879，y■=522.62741530，再把曲线看做是质点在空间中的运动轨迹，则质点的速度向量及加速度向量可以由曲线方程中的一阶及二阶倒数给出：速度向量v=（1，-0.4159，107.9675），加速度向量a=（0，-1.1492，-68.8471）。加速度可以分解为两部分：一个是与速度方向平行的切向加速度，另一个是与速度方向垂直的法向加速度。由圆周运动中速度与加速度的关系就可以确定圆周运动的半径，也就是曲率半径，从而求出曲率。圆周运动关系式为：|a|=|v|■/2。其中r是圆周运动的半径。切向加速度大小就是加速度向量在速度向量上的投影值，由向量的运算规则可以求出切向加速度大小为：（v*a）/|v|，则法向加速度大小为：■，曲率半径为：r=|v|^3/■。曲率k■=■/|v|^3。再把速度向量、加速度向量代入，解得k■=0.000444。同理可解得1996年、2009年、2011年的曲率分别为k■=0.000481，k■=0.000141，k■=0.000142411。

类似地，可求得各年份在塔高30米、45米的曲率分别为：

K■=0.000423，K■=0.000459，K■=0.000142，K■=0.000144076.（30米处）

K■′=0.000405，K■′=0.000441，K■′=0.000145，K■′=0.000147087.（45米处）

由下图可分析出从1986年至1996年曲率大幅降低，1996年至2009年、2009年至2011年都呈增长趋势，可能是随着观测年份的增加，古塔的年龄增长，故弯曲的程度会随之提高。这充分证明随着年份的增加，塔身离1986年的基准越大，且倾斜的幅度使先升高后减小，因此若塔继续保持这种势态塔就不会倒塌。

四、模型评价

在建模过程中我们运用算术平均和加权平均，使多项数据简化，方便于模型求解。在数据处理上，运用Matlab的线性拟合、二次多项式拟合、三次多项式拟合等方法求得空间直线方程。并给出了斜率公式、曲率公式计算古塔的倾斜、弯曲和扭曲等变形情况，还通过画图、列表等方法对倾斜、弯曲、扭曲这三种情况进行了数据分析。不足的地方是，由于时间关系未能对其他模型进行优化处理。

参考文献：

[1]史峰.MATLAB函数速查手册.北京：中国铁道出版社，2011.

[2]薛毅.数学建模基础.北京：北京工业大学出版社，2005.

[3]张德丰.MATLAB数值分析与应用.北京：国防工业出版社，2009.

[4]叶其孝.大学生数学建模竞赛辅导材料[M].长沙：湖南教育出版社，2002.

[5]张序.测量学.南京：东南大学出版社，2007.

基金项目：中山职业技术学院教研教改课题《与电子类专业教学相结合的〈电路数学〉探究型项目式教学的研究》（编号：JYC1110），主持人：欧冰。