基于灰色模型的古塔變形分析預測

孫 強，肖 云

(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.