唐紀(jì)芳，張子衛(wèi)

（四川職業(yè)技術(shù)學(xué)院,四川 遂寧 629000）

關(guān)于古塔變形研究的數(shù)學(xué)模型

唐紀(jì)芳，張子衛(wèi)

（四川職業(yè)技術(shù)學(xué)院,四川 遂寧 629000）

我國(guó)是一個(gè)文明古國(guó)，各地的古塔建筑非常的多，由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等.為保護(hù)古塔，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，了解各種變形量，以制定必要的保護(hù)措施.本文通過(guò)對(duì)古塔測(cè)量數(shù)據(jù)進(jìn)行分析研究，對(duì)古塔各層的中心、變形的塔身傾斜度、扭曲度、古塔彎度等給出了明確的定義，并運(yùn)用數(shù)學(xué)軟件進(jìn)行了相應(yīng)計(jì)算，對(duì)古塔的各個(gè)測(cè)試時(shí)間的變形程度給出了具體量化結(jié)論，并建立了古塔變形的數(shù)學(xué)模型預(yù)測(cè)其變形趨勢(shì)，給文物管理部門提出了有關(guān)古塔保護(hù)的建議.

古塔變形；各層中心點(diǎn)；傾斜度；彎曲度；扭曲度；變形預(yù)測(cè)

1 問(wèn)題（2013年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽C題）

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

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

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

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

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

（3）分析該塔的變形趨勢(shì).

2 模型假設(shè)

(1)假設(shè)附件1中的數(shù)據(jù)是以同樣的方式對(duì)同樣的8個(gè)點(diǎn)進(jìn)行了4次觀測(cè).

(2)假設(shè)古塔最初建筑時(shí)是正八邊形，而它的中心就是重心.

(3)假設(shè)塔尖就是一個(gè)點(diǎn).

(4)假設(shè)底層不考慮彎曲和扭曲程度.

(5)假設(shè)各層的傾斜、彎曲和扭曲均是以第1層為基礎(chǔ)作比較.

(6)假設(shè)古塔是公元1000年修建.

(7)假設(shè)古塔的變形不受強(qiáng)烈地質(zhì)災(zāi)害影響.

3 模型建立與求解

3.1 問(wèn)題1的求解

（1）分析.根據(jù)附件提供的數(shù)據(jù)信息，發(fā)現(xiàn)所測(cè)數(shù)據(jù)為古塔各層底部角點(diǎn)處的空間坐標(biāo)數(shù)據(jù)，經(jīng)過(guò)試算，1至13層為邊長(zhǎng)相等的空間八邊形，因此我們認(rèn)為，古塔修建時(shí)，1至13層為水平的平面正八邊形，正八邊形的中心就是各頂點(diǎn)的重心，由于時(shí)間久遠(yuǎn)，長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，從而產(chǎn)生變形，各層底面已經(jīng)不再是正多邊形，頂點(diǎn)已經(jīng)變化位置，但中心應(yīng)該在各層底面各頂點(diǎn)的重心附近，因此我們先根據(jù)《圖論》中的知識(shí)得出空間有限個(gè)點(diǎn)的中心的定義，再運(yùn)用matlab軟件編程在各層底部各頂點(diǎn)的重心附近搜尋中心位置.

（2）空間個(gè)n點(diǎn)Ai(i=1,2,…,n)的中心O的定義.若空間一點(diǎn)O滿足min﹛ max﹛ OAi﹜﹜(i=1,2,…,n)，則O叫做空間n個(gè)點(diǎn)Ai(i=1,2,…,n)的中心.即對(duì)空間上的任一點(diǎn)到這n個(gè)點(diǎn)的距離中有一個(gè)最大值，則所有可能的最大值中那個(gè)最小的值所對(duì)應(yīng)的那個(gè)點(diǎn)我們就把它叫做這n個(gè)點(diǎn)Ai(i=1,2,…,n)的中心.

（3）各層底面中心的確定.在研究各層的中心之前，我們首先對(duì)缺失的數(shù)據(jù)進(jìn)行了補(bǔ)充.在1986年和1996年的第十三層的第5個(gè)觀測(cè)點(diǎn)的數(shù)據(jù)都缺失了.我們將2009年和2011年該層的同一個(gè)觀測(cè)點(diǎn)的數(shù)據(jù)差值除以年份差作為每年的變化差值，然后計(jì)算出1986年和1996年的補(bǔ)充數(shù)據(jù).然后，各層底面中心應(yīng)該在各層8個(gè)觀測(cè)點(diǎn)的重心附近，所以我們先求出其重心坐標(biāo)，以其為中心構(gòu)造長(zhǎng)方體型區(qū)域，在這個(gè)區(qū)域內(nèi)編程搜索滿足min﹛ max﹛ OAi﹜﹜(i=1,2,…,n)的點(diǎn)O，其程序見附錄1，計(jì)算結(jié)果如下（表1）：

表1 各年所測(cè)古塔各層中心坐標(biāo)

我們運(yùn)用matlab軟件的三維作圖工具作出古塔在各年所測(cè)的頂點(diǎn)與中心的圖形，見圖一.

圖一

3.2 問(wèn)題2的求解

（1）古塔傾斜角的計(jì)算.我們認(rèn)為，古塔的傾斜分為塔基的傾斜與第i層的傾斜，而傾斜度的刻畫分為與鉛直方向的偏離角度和偏向水平面的具體方位，這個(gè)方位我們用塔身向量在平面上的投影與軸的正方向的夾角來(lái)確定.底層傾斜角的計(jì)算我們采用將底層各頂點(diǎn)擬合的平面的法向量與鉛直方向的偏離角度和偏向水平面的具體方位來(lái)確定，第i層的傾斜程度我們用第1層的中心為起點(diǎn)，第i層的中心為終點(diǎn)的向量與鉛直方向的偏離角度和偏向水平面的具體方位來(lái)刻畫.（見下圖）

1）第一層底面（塔基）的傾斜度

我們先求出各頂點(diǎn)的擬合平面，設(shè)擬合平面方程為Axil+Byil+Czil=1(C＞0)，代入八個(gè)頂點(diǎn)坐標(biāo)得到八個(gè)關(guān)于A,B,C的一次方程，然后用matlab軟件解超定方程的矩陣左除法得到擬合平面的法向量（程序見附錄3），m=﹛A,B,﹜C得到它與鉛直方向向量{0，0,1}的夾角度數(shù)

在水平面上偏離ox軸的正方向的夾角為

2）古塔第1層至第j層的傾斜度的計(jì)算

表2 各年豎直傾斜度α

表3 各年水平偏轉(zhuǎn)度θ

我們?nèi)「鲗拥钠骄鶅A斜度為古塔的傾斜度，結(jié)果見上表.

（2）古塔扭曲度的計(jì)算.我們?cè)诿繉铀臏y(cè)量數(shù)據(jù)中選定第一個(gè)觀測(cè)點(diǎn)與該層的中心點(diǎn)連線成一個(gè)空間向量，然后將古塔第j-1層與第j層兩層的這兩個(gè)向量的夾角作為古塔第j-1層的局部扭曲度，運(yùn)用表2數(shù)據(jù)我們得到如下結(jié)果，見表4.

表4 扭曲度

1986年最大扭曲度在第十層，扭曲度為1.467943度，1996年最大扭曲度在第十層，扭曲度為1.445559度，2009年最大扭曲度在第五層，扭曲度為0.917536度，2011年最大扭曲在第五層，扭曲度為0.93176度.

（3）古塔彎曲度的計(jì)算.古塔在第j層的彎曲度，我們認(rèn)為是由Oj-1,Oj,Oj+1這三個(gè)點(diǎn)確定。這三個(gè)點(diǎn)可確定一個(gè)圓，把這個(gè)圓的半徑的倒數(shù)作為古塔在第j層的彎曲度.

接下來(lái)我們計(jì)算彎曲度，我們以第1,2,3層為例：

表5 各年彎曲度對(duì)比表

1986年、1996年、2009年、2011年最大扭曲分別發(fā)生在第11層、第10層、第9層、第8層，彎曲度為分別為：

3.3 問(wèn)題3的求解

古塔變形趨勢(shì)的預(yù)測(cè)我們假設(shè)古塔建于公元1000年，其傾斜度，扭曲度、彎曲度均為0，我們以時(shí)間t為變量，建立第t年古塔最大傾斜度，最大扭曲度、最大彎曲度的函數(shù)f(t),g(t),h(t).

首先，取古塔每年最大豎直及水平偏轉(zhuǎn)度，建立了表6，再用這5個(gè)點(diǎn)進(jìn)行回歸，，通過(guò)回歸方法求出傾斜函數(shù)f(t)，然后根據(jù)此函數(shù)計(jì)算，得出2011——2026年預(yù)測(cè)數(shù)據(jù)見表7，作出相應(yīng)的圖形，見圖二、圖三.

表6

傾斜函數(shù)：

表7

圖二

圖三

取古塔每年最大扭曲度建立了表8，利用這5個(gè)點(diǎn)進(jìn)行回歸分析，通過(guò)回歸方法求出扭曲函數(shù)，即最大扭曲度函數(shù)：，然后根據(jù)函數(shù)計(jì)算，得出2012——2026年預(yù)測(cè)數(shù)據(jù)見表9，作出相應(yīng)的圖形，見圖四.

表8

表9

圖四

表10

表11

圖五

4 模型評(píng)價(jià)與改進(jìn)

本文在實(shí)際問(wèn)題的基礎(chǔ)上，通過(guò)合理的假設(shè)，預(yù)測(cè)出古塔的最大可能的傾倒方向，可以協(xié)助有關(guān)部門對(duì)古跡的維護(hù).

充分利用excle來(lái)對(duì)大量的數(shù)據(jù)進(jìn)行計(jì)算，完成表格的輸入、統(tǒng)計(jì)、分析等多項(xiàng)工作，可生成精美直觀的表格、圖表.

Matlab是一種科學(xué)計(jì)算軟件，專門以矩陣的形式處理數(shù)據(jù)，利用matlab計(jì)算可以使我們所計(jì)算的數(shù)據(jù)更加精確.

[1]陳國(guó)華,韋程?hào)|.數(shù)學(xué)建模與數(shù)學(xué)建模方法[M].天津：南開大學(xué)出版社，2012.

[2]陳光亭,裘哲勇.數(shù)學(xué)建模[M].北京：高等教育出版社，2012.

[3]梁?？?古塔變形測(cè)量方法探討[J].城市勘測(cè)，2011，（3）.

[4]吳望名.圖論及其應(yīng)用[M].北京：科學(xué)出版社，1984.

責(zé)任編輯：張隆輝

O 29

A

1672-2094(2013)06-0155-05

2013-10-28

唐紀(jì)芳（1982-），女，四川遂寧人，四川職業(yè)技術(shù)學(xué)院應(yīng)用數(shù)學(xué)與經(jīng)濟(jì)系講師，碩士。

張子衛(wèi)（1963-），男，四川遂寧人，四川職業(yè)技術(shù)學(xué)院應(yīng)用數(shù)學(xué)與經(jīng)濟(jì)系副教授。