孙 强，肖 云

(1.宝鸡文理学院 数学系，陕西 宝鸡 721016；2.陕西能源职业技术学院 基础一部，陕西 咸阳 712000)

1 提出问题(2013高教社杯全国大学生数学建模竞赛C题[1])

由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。

某古塔已有上千年历史，是我国重点保护文物。管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。

请你们根据附件1提供的4次观测数据，讨论以下问题：

1.给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。

2.分析该塔倾斜、弯曲、扭曲等变形情况。

3.分析该塔的变形趋势。

2 解决问题

2.1 各层中心坐标

根据1986年第一层8个观测点画出的平面图(见图1)可以看出该古塔为近似的凸正八边形[2]。根据平面解析几何[3]的知识，古塔每层的中心就是每层各顶点坐标的算术平均数。

图1 平面投影图

由于1986年和1996年13层缺少观测点5的数据，所以分别选择观测点4与观测点6的平均值作为观测点5的数据。从而得到各年各层中心点坐标，如表1所示。

表1 中心点坐标表

2.2 关于倾斜问题的计算

以1986年为例，由各层几何中心的横坐标和竖坐标拟合出在xoz面内的直线方程z=95x-53715，即为空间直线在该平面内的投影；由各层几何中心的纵坐标和竖坐标拟合出在yoz面内的直线方程z=-132y+69185，即为空间直线在yoz面内的投影。过直线z=95x-53715作xoz面的垂面，过直线z=-132y+69185作yoz面的垂面，两个垂面的交线即为古塔中心所在的空间直线L1，

利用方向余弦公式可以得到1986年古塔的倾斜角度为：γ1986=44.58170025816806'

同理可得其它各年古塔的倾斜角度(见表2)。

表2 各年古塔倾斜角度表

通过上表数据发现，古塔在逐年倾斜，尤其是1996年到2009年倾斜幅度比较大，可能是由于地震、飓风等因素的影响，加剧了古塔的倾斜。

2.3 关于古塔弯曲的计算

古塔矗立多年，其塔身已经发生细微变化，若通过对各层的中心坐标进行非线性拟合得到空间曲线，其难度较大，因此转换思路，考虑古塔在三维坐标系的两侧面的投影，显然投影为平面曲线，对两投影平面曲线的弯曲程度进行分析，从而确定空间曲线的弯曲程度。

以1986年为例，由各层的几何中心的横坐标和竖坐标拟合出在xoz面内的投影为：

z= -21.832240269999999782157829031348x2

+24850.542571490001137135550379753x

-7071397.5216804295778274536132812

由各层的几何中心的纵坐标和竖坐标拟合出在yoz面内的投影为：

z=-170.93657200000001239459379576147y2

+178493.15410139999585226178169250y

-46595886.313481502234935760498047

平面曲线的曲率计算公式[4]为：

表3 1986年各层中心点弯曲情况

同理可得各年各层中心点的曲率，如表4、表5、表6所示。

表4 1996年各层中心点弯曲情况

续表4

层数xyKxKy4566．8183522．5922Kx=4．242267564159704e-005Ky=7．304663911800054e-0055566．8649522．5563Kx=4．501587261097414e-005Ky=9．113076805948330e-0056566．9118522．521Kx=4．501587261097414e-005Ky=1．152331437252489e-0047566．9506522．5042Kx=5．036472775651844e-005Ky=1．297267560566100e-0048566．9884522．4881Kx=5．299409099587878e-005Ky=1．459892255666876e-0049567．0265522．4714Kx=5．583284127915918e-005Ky=1．658620028790178e-00410567．062522．4572Kx=5．866326239706748e-005Ky=1．856968285815020e-00411567．1102522．4173Kx=6．282091068587297e-005Ky=2．614819568214344e-00412567．1578522．3775Kx=6．732146680390316e-005Ky=3．844106721752882e-00413567．1812522．4504Kx=6．969212705308391e-005Ky=1．963212757186078e-004塔尖567．2543522．2366Kx=7．785485781506534e-005Ky=0．00281763872959

表5 2009年各层中心点弯曲情况

表6 2011年各层中心点弯曲情况

续表6

层数xyKxKy7566．98522．5115Kx=1．391862369344460e-004Ky=1．315042216570195e-0048567．0313522．4788Kx=1．652365186907063e-004Ky=1．588277771021799e-0049567．0825522．4457Kx=1．981423066226998e-004Ky=1．947726404754545e-00410567．1381522．3926Kx=2．445711068920006e-004Ky=2．789927810965960e-00411567．181522．3535Kx=2．908117973953012e-004Ky=3．747825722119750e-00412567．2238522．3147Kx=3．493289592971776e-004Ky=5．185771509623188e-00413567．2725522．2701Kx=4．366927098082264e-004Ky=7．913512996982270e-004塔尖567．3375522．2135Kx=6．052506961316644e-004Ky=0．00150009354147

2.4 古塔扭曲的计算

随着层数的增加，扭曲程度越来越大，故取13层作为研究对象，分析古塔的扭曲情况。

首先考虑1986年到1996年的扭曲情况。选取1986年和1996年第13层的8个观测数据与当年第13层的几何中心构成的两向量的夹角，就可以反映出第13层的扭曲程度。

图2 第13层扭曲角度

其中A1，A2，A3，A4，A5，A6，A7，A8和O为1986年13层的观测点和中心点；A1′，A2′，A3′，A4′，A5′，A6′，A7′，A8′和O′为1996年13层的观测点和中心点；

以此类推，可以算出其余七个夹角，最后求这八个夹角的均值，它揭示了1986年到1996年古塔的扭曲程度。

但是，按原数据计算1996年到2009年发生的扭曲约为89.73o，显然这不符合实际情况。通过作图发现2009年观测点编号的顺序与1996年观测点编号顺序发生紊乱，故对2009年数据编号进行调整(见表7)。

表7 2009年13层数据调整表

用类似的方法可以分别计算出1996年到2009年和2009年到2011年的扭曲角度(见表8)。

2.5 变形趋势的预测

GM(1,1)模型由于其实用性和可操作性，是迄今为止应用最为广泛的灰色模型之一。近年来它被广泛应用于工程、经济、图像处理、环境等许多领域。

表8 古塔的扭曲程度表

建立灰色预测模型[5]时，需要对原始数据进行累加处理，通过累加处理对非负序列找到某种规律，然后建立微分方程。若X(0)={x(0)(1),x(0)(2),…,x(0)(n)|n≥4}是一个原始的非负连续等间距数据序列，GM(1,1)的具体步骤如下所示：

第一步：原始数据X(1)={x(1)(1),x(1)(2),…,x(1)(n)}序列通过依次累加生成(AGO)。这里

一阶差分方程的离散化得到：

x(0)(k)+az(1)(k)=b

这里背景值被定义为：

z(1)(k)=0.5×x(1)(k)+0.5×x(1)(k-1)

第三步：通过最小二乘法估计发展系数a和灰输入b。

YN=[x(0)(2),x(0)(3),…,x(0)(n)]T

第四步：对于预测模型确立以下一阶差分的白化表达式：

第五步：求解等式并从AGO中获得序列X(1)的预测值如下所示：

对上面等式通过累减生成运算(IAGO)，

x(0)(k)=x(1)(k)-x(1)(k-1)就能够得到原始数据序列的预测值，即：

k=2,3,…,

用GM(1,1)模型，可以得到古塔未来几年的变化趋势(见表9)。

表9 GM(1,1)预测古塔未来中心变化趋势值

通过对古塔未来变化的预测，可以很好地揭示古塔的变形情况及发展趋势。通过表中数据可以看出，其变化趋势和我们上面所做的倾斜、弯曲、扭曲情况基本一致，需要对古塔及时维护与加固处理，减少人民生命与财产的损失。但模型还有不足之处，需要再进一步的研究，如能够绘制出空间曲线，那么对该问题的理解就更加清晰明确了。

[1] 2013年高教社杯全国大学生数学建模竞赛赛题[EB/OL]. http://www.mcm.edu.cn/problem/2013/2013.html .

[2] 胡志晓.古塔倾斜观测和数据分析[J].江苏建筑，2011(6):34-35，44.

[3] 杨文茂，李全英.空间解析几何[M].武汉：武汉大学出版社，2006:30-32.

[4] 华东师范大学数学系.数学分析[M].北京：高等教育出版社，2011:250-251.

[5] 卓金武. MATLAB在数学建模中的应用[M].北京：北京航空航天大学出版社，2011:30-32.