數學建模項目實踐

楊晶晶　歐冰

摘 要： 為及時了解古塔變形的情況，制定必要的保護措施，作者首先通過四年的觀測數據建立數學模型，擬合出古塔各個塔層的中心坐標、中心軸線，再依據中心軸線的變化分析古塔的傾斜、彎曲、扭曲等變化情況和變化趨勢。

關鍵詞： 數學模型 數據擬合 斜率 曲率

一、問題的提出與分析

由于長時間承受自重、氣溫、風力等各種作用，偶然還要受地震、颶風的影響，古塔會產生各種變形，諸如傾斜、彎曲、扭曲等。為保護古塔，文物部門需適時對古塔進行觀測，了解各種變形量，制定必要的保護措施。某古塔已有上千年歷史，是我國重點保護文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進行了4次觀測。文物部門希望通過這4次觀測數據確定古塔各層中心位置和古塔各層中心坐標，同時分析該塔傾斜、彎曲、扭曲等變形情況。

為解決此問題，我們需要建立相關數學模型，通過模型的解答與分析判斷古塔的心位置，古塔各層中心坐標，以及該塔傾斜、彎曲、扭曲等變形情況。為方便建立數學模型，我們假設古塔各個塔層和塔尖的高度都相同；假設古塔每層有8個觀測點且8個點在各層內離層底相對高度相同；假設每個年份對每一層的觀測點都固定；假設每層觀測點在古塔的中心軸周圍。通過把每一層的相關觀測數據進行算術平均化處理，通過對每一層的不同年份的數據進行算術平均，得出13個塔層4個年份的三維坐標，最后運用數據擬合的方法得到古塔各個塔層和塔尖的中心坐標方程及中心軸線方程。然后分析塔的傾斜情況。我們對已知的多點觀測數據進行算術平均處理，把處理后的數據進行多項式擬合后得出擬合方程組，構建空間直角坐標系，通過勾股定理并且進行數角轉換便可得出傾斜角度。最后分析塔的彎曲情況。通過對已知數據進行2次多項式擬合，選取低、中、高三個層次利用曲率比較分析古塔4個年份的彎曲程度。曲率越大，表示古塔的彎曲程度越大。

二、古塔中心坐標、中心軸線的求解

我們把各測量數據進行算術平均化處理，通過數據擬合得出各年份的各層擬合中心位置坐標，以第一層為例，1986年的中心點坐標為（566.657142，522.70516，1.787375），1996年的中心點坐標為（566.66，493.1，1.78），2009年的中心點坐標為（566.67，522.74，1.76），2011年的中心點坐標為（566.67，522.74，1.76）。其余各層可用類似的方法得到中心點坐標。

為了減少誤差，使結果更精確，我們對得出的4個年份的中心點坐標再進行一次算術平均處理和數據擬合處理，即可得出第一層的中心點坐標為M1（566.67，522.71，2.13），第二層的中心點坐標為M2（566.72，522.68，6.38）。其余各層的中心點坐標可類似得出。通過各層的中心點坐標可以直觀看出各個塔層的水平方向的中心坐標沒有發生較大偏差，說明塔層越高偏移相對性大且呈現出傾斜的趨勢。

在上述模型的求解過程中，我們使用的是算術平均，經過實踐，我們發現運用加權平均的處理方法可以減少塔中心坐標的誤差，使模型更加貼近真實情況。我們對年份進行一個加權平均處理，隨著年份的增加，所占權重越?。晃覀儼褭嘀卦O為1986年占40%，1996占30%，2009年占20%，2011年占10%。經過類似計算處理可得到各塔層中心坐標為：

M1（566.7，522.7，2.1），M2（566.7，522.7，6.4），M3（566.8，522.7，10.6），

M4（566.8，522.6，14.9），

M5（566.9，522.6，19.1），

M6（566.9，522.6，23.4），

M7（566.9，522.5，27.7），

M8（566.99，522.5，31.9），

M9（567.1，522.5，36.2），

M10（567.1，522.5，40.4），

M11（567.1，522.4，44.7），

M12（567.2，522.4，48.9），

M13（567.2，522.4，53.2），M14（567.2，522.3，55.1）。

如上圖所示：從圖形上不但可以直觀地看出古塔空間中軸線、各個塔層段的中心點，更可以通過塔層中心點的連線判斷出古塔的變形趨勢。

三、古塔的變形

我們利用上述所求得的古塔中心軸線，建立空間三維坐標系，再通過數據擬合出四條空間直線，即可求得斜率。

對于古塔彎曲的分析要從古塔的曲率入手。首先我們用Matlab軟件分別對1986年7月、1996年8月、2009年3月和2011年3月這四次觀測的數據進行擬合，擬合成四個方程組。我們以1986年的多項式為例，進行空間曲線的曲率計算。

1986年：

y■=-0.574605251x■■+650.946992376x■-183835.161349422z■=-34.4235380x■■+39129.9111344x■-11119834.7987881

首先討論其在P點（x1，y1，15）的曲率，將z■=15代入方程，解得x■=566.791589879，y■=522.62741530，再把曲線看做是質點在空間中的運動軌跡，則質點的速度向量及加速度向量可以由曲線方程中的一階及二階倒數給出：速度向量v=（1，-0.4159，107.9675），加速度向量a=（0，-1.1492，-68.8471）。加速度可以分解為兩部分：一個是與速度方向平行的切向加速度，另一個是與速度方向垂直的法向加速度。由圓周運動中速度與加速度的關系就可以確定圓周運動的半徑，也就是曲率半徑，從而求出曲率。圓周運動關系式為：|a|=|v|■/2。其中r是圓周運動的半徑。切向加速度大小就是加速度向量在速度向量上的投影值，由向量的運算規則可以求出切向加速度大小為：（v*a）/|v|，則法向加速度大小為：■，曲率半徑為：r=|v|^3/■。曲率k■=■/|v|^3。再把速度向量、加速度向量代入，解得k■=0.000444。同理可解得1996年、2009年、2011年的曲率分別為k■=0.000481，k■=0.000141，k■=0.000142411。

類似地，可求得各年份在塔高30米、45米的曲率分別為：

K■=0.000423，K■=0.000459，K■=0.000142，K■=0.000144076.（30米處）

K■′=0.000405，K■′=0.000441，K■′=0.000145，K■′=0.000147087.（45米處）

由下圖可分析出從1986年至1996年曲率大幅降低，1996年至2009年、2009年至2011年都呈增長趨勢，可能是隨著觀測年份的增加，古塔的年齡增長，故彎曲的程度會隨之提高。這充分證明隨著年份的增加，塔身離1986年的基準越大，且傾斜的幅度使先升高后減小，因此若塔繼續保持這種勢態塔就不會倒塌。

四、模型評價

在建模過程中我們運用算術平均和加權平均，使多項數據簡化，方便于模型求解。在數據處理上，運用Matlab的線性擬合、二次多項式擬合、三次多項式擬合等方法求得空間直線方程。并給出了斜率公式、曲率公式計算古塔的傾斜、彎曲和扭曲等變形情況，還通過畫圖、列表等方法對傾斜、彎曲、扭曲這三種情況進行了數據分析。不足的地方是，由于時間關系未能對其他模型進行優化處理。

參考文獻：

[1]史峰.MATLAB函數速查手冊.北京：中國鐵道出版社，2011.

[2]薛毅.數學建?；A.北京：北京工業大學出版社，2005.

[3]張德豐.MATLAB數值分析與應用.北京：國防工業出版社，2009.

[4]葉其孝.大學生數學建模競賽輔導材料[M].長沙：湖南教育出版社，2002.

[5]張序.測量學.南京：東南大學出版社，2007.

基金項目：中山職業技術學院教研教改課題《與電子類專業教學相結合的〈電路數學〉探究型項目式教學的研究》（編號：JYC1110），主持人：歐冰。