古塔變形分析模型

李 靜，張 媛，顧 忠，賀 欣，王 鵬

(1.鄭州鐵路職業(yè)技術學院，河南 鄭州 450052；2. 鄭州市金水區(qū)教育體育局, 河南 鄭州 450008)

古塔變形分析模型

李 靜1，張 媛1，顧 忠1，賀 欣1，王 鵬2

(1.鄭州鐵路職業(yè)技術學院，河南 鄭州 450052；2. 鄭州市金水區(qū)教育體育局, 河南 鄭州 450008)

根據2013年全國大學生數學建模比賽C題中某古塔4年的觀測數據，給出了確定古塔各層中心位置的通用方法，建立了最優(yōu)化模型，用Lingo軟件求得4次測量的古塔的各層中心坐標。以斜率、曲率、投影、均方差、擬合等知識為基礎，對傾斜度、彎曲度、扭曲度三個指標進行定義，結合使用Excel和Matlab軟件對古塔變形情況進行量化分析，最后根據得到的數據對古塔的變形趨勢進行預測。

傾斜度；彎曲度；扭曲度；均方差；優(yōu)化模型；曲率；擬合

由于古塔長時間受自身或外界因素影響，會發(fā)生各種變形，包括傾斜、彎曲、扭曲。現對古塔4年中個別月份變形情況統(tǒng)計數據進行分析。

1 問題一的模型建立與求解

對于2013年全國大學生數學建模比賽C題(以下簡稱“C題”)中的問題一，從物理結構的穩(wěn)定角度考慮，想象成對8個點結繩拉緊，繩長越短繩子越緊，求塔每層中心位置可以轉化為求某一點到塔每層各個點距離之和最短。我們建立了最優(yōu)化模型，運用Lingo軟件求解出4年古塔每層的中心坐標，整理繪制成Excel表格并用Matlab軟件制圖。

1.1 繪制古塔三維立體圖形

根據C題附件1給出1986年7月古塔x，y，z坐標，運用Matlab軟件對數據進行預處理并且制出1986年古塔三維立體圖(如圖1)，以方便對問題進行思考與研究。

圖1 古塔的基本形狀

1.2 求古塔每層中心坐標

由于古塔常年受自身重量以及自然環(huán)境影響，發(fā)生傾斜、彎曲、扭曲等形變，致使原來在同一平面的8個頂點不在同一平面上。

我們想象成用一根繩子從中心位置處分別和8個點結繩，若要結構最穩(wěn)定，也就是使結繩最緊，即所結繩長之和最短。利用這個思想，可以把求塔每層中心目標轉化成求某一點到塔層各個端點距離之和最短，若求各層中心坐標可轉化成某點到該層8個頂點距離之和最小，求這點坐標分量。

建立最優(yōu)化模型一：

這是個非線性規(guī)劃模型，運用Lingo軟件編程求解，并將4年的各塔層中心目標數據用Excel軟件進行統(tǒng)計制表(如表1)。

表1 各塔層中心坐標 單位：m

續(xù)表

層數 2009年3月xyz 2011年3月xyz1566．7089522．70261．7645566．7091522．70251．76332566．7490522．66667．3091566．7492522．66647．29063566．7885522．631412．7324566．7889522．631812．72704566．8201522．603217．0697566．8206522．602717．05215566．8542522．573021．7095566．8548522．572421．70406566．9384522．516826．2114566．9391522．516126．20507566．9693522．495429．8250566．9701522．494629．81738567．0192522．462133．3403567．0200522．461233．33709567．0688522．428636．8441567．0697522．427636．822610567．1248522．376540．1612567．1258522．375240．144211567．1681522．338944．4328567．1692522．337744．425112567．2112522．301548．7000567．2126522．300248．684013567．2598522．258152．8188567．2611522．256752．8135塔尖567．3360522．214855．0910567．3375522．213555．0870

2 問題二的模型建立與求解

對于C題問題二，我們首先定義了傾斜度、彎曲度、扭曲度3個指標，并給出計算這些指標的公式。接下來分別通過3個指標對古塔變形情況進行分析。

考慮傾斜變形時，我們細化出總偏移量、南北偏移量、東西偏移量、總傾斜度、南北傾斜度、東西傾斜度等指標，利用各層中心坐標，刻畫出古塔4年傾斜程度，得出在2009年傾斜度和偏移量最大。

考慮彎曲變形時，借用均方差的思想，利用各層中心點坐標數據，將古塔14層分為13段，求出古塔每段斜率，用各段斜率與上下中心連線的斜率的差的平方和開方這一偏離程度作為彎曲度，偏離程度越大，彎曲度越大，得出在2009年彎曲度最大。

考慮扭曲變形時，把每年古塔各角在每層的觀測點投影到平面xOy上，擬合出8條曲線，利用這8條曲線的平均曲率的總和來刻畫古塔每年的扭曲度，得出在1996年扭曲度最大。

2.1 古塔傾斜情況分析

結合問題一得出的各塔層中心坐標，建立空間直角坐標系(如圖2)，這里塔尖中心點P投影到塔底中心點O所在平面z=zO上，記為點Q，QM ⊥CM，QN ⊥CN，根據相應的空間幾何知識可得出總偏移量、南北偏移量、東西偏移量、總偏移度、南北偏移度、東西偏移度。

圖2 塔頂塔底中心坐標示意圖

將4年的塔底和塔頂中心坐標代入以下各傾斜指標:

利用MATLAB軟件對以上表達式進行編程求解，得出傾斜指標具體數據，同時根據中心坐標的移動得出古塔有南偏、東偏趨勢，將變形數據匯總整理得到表2。

表2 古塔傾斜指標匯總

年份 總偏移量m南偏移量m東偏移量m總傾斜度°南傾斜度°東傾斜度°高度m19860．75060．46810．58680．80620．50280．630355．128119960．76080．47480．59440．81710．51000．638455．124020090．79450．48780．62710．85360．52410．673755．091020110．79620．48900．62840．85550．52540．675255．0870

為了清晰觀察4年傾斜指標的變化，根據表2數據，運用Excel軟件分別繪制古塔方向偏移量折線圖(如圖3)和古塔方向傾斜度折線圖(如圖4)。

圖3 古塔4年偏移量走勢圖

圖4 古塔4年傾斜度走勢圖

由表2整理數據結合圖3、圖4分析出：古塔在1986年、1996年、2009年發(fā)生傾斜形變的大小，隨著古塔建立年代的久遠，偏移量、傾斜度變化呈逐步上升趨勢，在不采取對古塔變形預防或補救的情況下，由于古塔材質變形或損壞并且長時間承受自重、氣溫、風力等因素影響，古塔傾斜量只會逐年增加，滿足理論依據；總偏移量和傾斜角度變化最大為1996—2009年，其余年份偏移變化差值均勻增加，另外，1996—2009年塔的高度減少也較大。

古塔每年承受自重、氣溫、風力等影響是自然的(忽略地震、颶風等自然災害)，每年古塔發(fā)生形變量應大致呈現均勻增加趨勢，而1996—2009年差值較大可能因為2008年四川地震，此次大型自然災害使其發(fā)生傾斜程度較大，符合理論推斷。

2.2 古塔彎曲情況分析

具體分析古塔彎曲情況，即將古塔彎曲程度進行量化，用具體指標反映古塔4年中的彎曲形變大小，自定義變量彎曲度β。

這里以塔的各層中心點坐標作為參考數據，將古塔分為13段，塔尖—13層、13層—12層、…、3層—2層、2層—1層，如圖5。

圖5 中心坐標曲線各段斜率圖

結合圖6分析得出：1986年第11層、12層、13層斜率有下降趨勢，意味著古塔11層以上有彎曲，1996年第10層以上斜率減小，即古塔10層以上發(fā)生彎曲，而1986年、1996年前10層古塔斜率基本趨于穩(wěn)定，說明古塔前10層沒有發(fā)生彎曲；對于2009年、2011年斜率走勢圖基本吻合,在第5層、第7層、第9層斜率減小，表明古塔在這3層有不同程度的彎曲變形，其他各點沒有太大波動，10層以上繼續(xù)保持彎曲。

我們借用均方差思想，用各段斜率k(i)與k的偏差的平方和開方刻畫彎曲度

圖6 平均斜率走勢圖

用Maltab軟件計算斜率，運用Excel軟件對相應數據進行計算，同樣方法處理每年每塔層中心坐標得出所需的彎曲度(如表3)。

表3 4年各塔層彎曲度

根據表3數據明顯可以看出：隨著古塔建立時間的增加，彎曲度也明顯呈上升趨勢，與古塔傾斜偏移量變化趨勢大致一致。1996—2009年之間年彎曲度差值較大，2009—2011年時間較短，所以彎曲度沒有發(fā)生太大變化；整體走勢還是趨于穩(wěn)定型增長，符合天氣以及隨著年代的久遠古塔自身材質的改變與損壞規(guī)律。

2.3 古塔扭曲情況分析

針對古塔扭曲情況的分析，忽略彎曲等因素的影響。

把古塔各角每年的觀測點投影到平面xOy上，擬合出8條曲線，利用這8條曲線的平均曲率的和來刻畫古塔每年的扭曲度。

把擬合出的曲線y=ax3+bx2+cx+d的13個點(或12個點)代入曲率公式

計算出每條曲線的平均曲率，則每一年古塔的扭曲度定義為

γ越大，則扭曲程度越大；γ越小，扭曲程度越小。

把1986年7月13層的1號觀察點，投影在xOy平面上，運用Matlab軟件進行擬合，得出三次方程式

y=-8.3x3+1.4×104x2-7.9×106x+1.5×109， 然后將擬合方程代入曲率表達式，得出1號觀察點的平均曲率為0.073 073 846×10-10；同時可以求出其他7個觀察點的平均曲率，將這8個曲率相加即為扭曲度

γ=387.222 249 4×10-10。

同理將其他年份數據進行處理求得相應的扭曲度(如表4)。

表4 古塔4年扭曲度量化表

注：由于所得數據偏小，我們目的是分析古塔4年扭曲度變化，故將表中所得數據乘以 10-10。

我們得到了古塔4年的總扭曲度分別為：

387.222 249 4×10-10，68 503.476 02×10-10，93 806.860 49×10-10，93 821.802 53×10-10。

由表4分析得：扭曲度越來越大，古塔扭曲情況逐年嚴重，同樣與傾斜、彎曲呈現的變化趨勢相一致，而且從低層到高層扭曲度逐漸增加。

還可以看出：1986—1996年總體扭曲度變化最為明顯，1996年、2009年、2011年這些年間扭曲度變化相應較??；在1986年和1996年測量中，4號觀測點和7號觀測點各自對應的兩條棱的扭曲程度最顯著，其中，4號觀測點的棱的扭曲度變化幅度高于7號觀測點；在2009年和2011年測量中，2號觀測點和4號觀測點各自對應的2條棱的扭曲程度最為顯著，而7號觀測點對應的棱在這兩年的扭曲度反而有所下降。

[1]黃強. 古塔變形監(jiān)測的探討[J]. 測繪與空間地理信息，2013,36(6)：217-220.

[2]韓中庚. 數學建模實用教程[M]. 北京：高等教育出版社，2012.

[3](美)Mark M.Meerschaert.數學建模方法與分析[M].北京：機械工業(yè)出版社，2015.

[4]楊為民. 建筑物傾斜的計算方法的探討[J]. 施工技術與測量技術，2010,30(5)：175-176.

[責任編輯：趙 偉]

The Deformation Analysis Model of the Ancient Pagoda

LI Jing1, ZHANG Yuan1,GU Zhong1,HE Xin1, WANG Peng2

(1. Zhengzhou Railway Vocational and Technical College, Zhengzhou 450052,China;2. Zhengzhou Jinshui District Education Sports Bureau,Zhengzhou 450008,China)

According to the 2013 national college students' mathematical modeling competition problem C a pagoda four years of observation data, gives the general method of the ancient towers layers center position, the optimization model is established, using Lingo software, obtained four times while hiking in the measurement of each layer center coordinates; And by slope, curvature, projection, the mean square error, fitting, such as knowledge, slope, bending and torsion degrees "define" the three indexes and combined with using Excel and Matlab software, the quantitative analysis of the deformation of the tower, and according to the data to predict the deformation trend of tower.

gradient; bending; distortion; mean square error; optimization model;curvature; fitting

2016 - 03 - 07

李靜(1976— )，女，河北棗強人，鄭州鐵路職業(yè)技術學院公共教學部講師,研究方向為數學教學。

O29

A

1008-6811(2016)-04-0015-06