基于灰色模型的古塔變形分析預(yù)測

孫 強，肖 云

(1.寶雞文理學(xué)院 數(shù)學(xué)系，陜西 寶雞 721016；2.陜西能源職業(yè)技術(shù)學(xué)院 基礎(chǔ)一部，陜西 咸陽 712000)

1 提出問題(2013高教社杯全國大學(xué)生數(shù)學(xué)建模競賽C題[1])

由于長時間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護(hù)古塔，文物部門需適時對古塔進(jìn)行觀測，了解各種變形量，以制定必要的保護(hù)措施。

某古塔已有上千年歷史，是我國重點保護(hù)文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進(jìn)行了4次觀測。

請你們根據(jù)附件1提供的4次觀測數(shù)據(jù)，討論以下問題：

1.給出確定古塔各層中心位置的通用方法，并列表給出各次測量的古塔各層中心坐標(biāo)。

2.分析該塔傾斜、彎曲、扭曲等變形情況。

3.分析該塔的變形趨勢。

2 解決問題

2.1 各層中心坐標(biāo)

根據(jù)1986年第一層8個觀測點畫出的平面圖(見圖1)可以看出該古塔為近似的凸正八邊形[2]。根據(jù)平面解析幾何[3]的知識，古塔每層的中心就是每層各頂點坐標(biāo)的算術(shù)平均數(shù)。

圖1 平面投影圖

由于1986年和1996年13層缺少觀測點5的數(shù)據(jù)，所以分別選擇觀測點4與觀測點6的平均值作為觀測點5的數(shù)據(jù)。從而得到各年各層中心點坐標(biāo)，如表1所示。

表1 中心點坐標(biāo)表

2.2 關(guān)于傾斜問題的計算

以1986年為例，由各層幾何中心的橫坐標(biāo)和豎坐標(biāo)擬合出在xoz面內(nèi)的直線方程z=95x-53715，即為空間直線在該平面內(nèi)的投影；由各層幾何中心的縱坐標(biāo)和豎坐標(biāo)擬合出在yoz面內(nèi)的直線方程z=-132y+69185，即為空間直線在yoz面內(nèi)的投影。過直線z=95x-53715作xoz面的垂面，過直線z=-132y+69185作yoz面的垂面，兩個垂面的交線即為古塔中心所在的空間直線L1，

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

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

表2 各年古塔傾斜角度表

通過上表數(shù)據(jù)發(fā)現(xiàn)，古塔在逐年傾斜，尤其是1996年到2009年傾斜幅度比較大，可能是由于地震、颶風(fēng)等因素的影響，加劇了古塔的傾斜。

2.3 關(guān)于古塔彎曲的計算

古塔矗立多年，其塔身已經(jīng)發(fā)生細(xì)微變化，若通過對各層的中心坐標(biāo)進(jìn)行非線性擬合得到空間曲線，其難度較大，因此轉(zhuǎn)換思路，考慮古塔在三維坐標(biāo)系的兩側(cè)面的投影，顯然投影為平面曲線，對兩投影平面曲線的彎曲程度進(jìn)行分析，從而確定空間曲線的彎曲程度。

以1986年為例，由各層的幾何中心的橫坐標(biāo)和豎坐標(biāo)擬合出在xoz面內(nèi)的投影為：

z= -21.832240269999999782157829031348x2

+24850.542571490001137135550379753x

-7071397.5216804295778274536132812

由各層的幾何中心的縱坐標(biāo)和豎坐標(biāo)擬合出在yoz面內(nèi)的投影為：

z=-170.93657200000001239459379576147y2

+178493.15410139999585226178169250y

-46595886.313481502234935760498047

平面曲線的曲率計算公式[4]為：

表3 1986年各層中心點彎曲情況

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

表4 1996年各層中心點彎曲情況

續(xù)表4

層數(shù)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年各層中心點彎曲情況

續(xù)表6

層數(shù)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 古塔扭曲的計算

隨著層數(shù)的增加，扭曲程度越來越大，故取13層作為研究對象，分析古塔的扭曲情況。

首先考慮1986年到1996年的扭曲情況。選取1986年和1996年第13層的8個觀測數(shù)據(jù)與當(dāng)年第13層的幾何中心構(gòu)成的兩向量的夾角，就可以反映出第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年古塔的扭曲程度。

但是，按原數(shù)據(jù)計算1996年到2009年發(fā)生的扭曲約為89.73o，顯然這不符合實際情況。通過作圖發(fā)現(xiàn)2009年觀測點編號的順序與1996年觀測點編號順序發(fā)生紊亂，故對2009年數(shù)據(jù)編號進(jìn)行調(diào)整(見表7)。

表7 2009年13層數(shù)據(jù)調(diào)整表

用類似的方法可以分別計算出1996年到2009年和2009年到2011年的扭曲角度(見表8)。

2.5 變形趨勢的預(yù)測

GM(1,1)模型由于其實用性和可操作性，是迄今為止應(yīng)用最為廣泛的灰色模型之一。近年來它被廣泛應(yīng)用于工程、經(jīng)濟(jì)、圖像處理、環(huán)境等許多領(lǐng)域。

表8 古塔的扭曲程度表

建立灰色預(yù)測模型[5]時，需要對原始數(shù)據(jù)進(jìn)行累加處理，通過累加處理對非負(fù)序列找到某種規(guī)律，然后建立微分方程。若X(0)={x(0)(1),x(0)(2),…,x(0)(n)|n≥4}是一個原始的非負(fù)連續(xù)等間距數(shù)據(jù)序列，GM(1,1)的具體步驟如下所示：

第一步：原始數(shù)據(jù)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)

第三步：通過最小二乘法估計發(fā)展系數(shù)a和灰輸入b。

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

第四步：對于預(yù)測模型確立以下一階差分的白化表達(dá)式：

第五步：求解等式并從AGO中獲得序列X(1)的預(yù)測值如下所示：

對上面等式通過累減生成運算(IAGO)，

x(0)(k)=x(1)(k)-x(1)(k-1)就能夠得到原始數(shù)據(jù)序列的預(yù)測值，即：

k=2,3,…,

用GM(1,1)模型，可以得到古塔未來幾年的變化趨勢(見表9)。

表9 GM(1,1)預(yù)測古塔未來中心變化趨勢值

通過對古塔未來變化的預(yù)測，可以很好地揭示古塔的變形情況及發(fā)展趨勢。通過表中數(shù)據(jù)可以看出，其變化趨勢和我們上面所做的傾斜、彎曲、扭曲情況基本一致，需要對古塔及時維護(hù)與加固處理，減少人民生命與財產(chǎn)的損失。但模型還有不足之處，需要再進(jìn)一步的研究，如能夠繪制出空間曲線，那么對該問題的理解就更加清晰明確了。

[1] 2013年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽賽題[EB/OL]. http://www.mcm.edu.cn/problem/2013/2013.html .

[2] 胡志曉.古塔傾斜觀測和數(shù)據(jù)分析[J].江蘇建筑，2011(6):34-35，44.

[3] 楊文茂，李全英.空間解析幾何[M].武漢：武漢大學(xué)出版社，2006:30-32.

[4] 華東師范大學(xué)數(shù)學(xué)系.數(shù)學(xué)分析[M].北京：高等教育出版社，2011:250-251.

[5] 卓金武. MATLAB在數(shù)學(xué)建模中的應(yīng)用[M].北京：北京航空航天大學(xué)出版社，2011:30-32.