楊晶晶　歐冰

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

關(guān)鍵詞： 數(shù)學(xué)模型 數(shù)據(jù)擬合 斜率 曲率

一、問題的提出與分析

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

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

二、古塔中心坐標(biāo)、中心軸線的求解

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

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

在上述模型的求解過程中，我們使用的是算術(shù)平均，經(jīng)過實(shí)踐，我們發(fā)現(xiàn)運(yùn)用加權(quán)平均的處理方法可以減少塔中心坐標(biāo)的誤差，使模型更加貼近真實(shí)情況。我們對(duì)年份進(jìn)行一個(gè)加權(quán)平均處理，隨著年份的增加，所占權(quán)重越?。晃覀儼褭?quán)重設(shè)為1986年占40%，1996占30%，2009年占20%，2011年占10%。經(jīng)過類似計(jì)算處理可得到各塔層中心坐標(biāo)為：

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）。

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

三、古塔的變形

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

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

1986年：

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

首先討論其在P點(diǎn)（x1，y1，15）的曲率，將z■=15代入方程，解得x■=566.791589879，y■=522.62741530，再把曲線看做是質(zhì)點(diǎn)在空間中的運(yùn)動(dòng)軌跡，則質(zhì)點(diǎn)的速度向量及加速度向量可以由曲線方程中的一階及二階倒數(shù)給出：速度向量v=（1，-0.4159，107.9675），加速度向量a=（0，-1.1492，-68.8471）。加速度可以分解為兩部分：一個(gè)是與速度方向平行的切向加速度，另一個(gè)是與速度方向垂直的法向加速度。由圓周運(yùn)動(dòng)中速度與加速度的關(guān)系就可以確定圓周運(yùn)動(dòng)的半徑，也就是曲率半徑，從而求出曲率。圓周運(yùn)動(dòng)關(guān)系式為：|a|=|v|■/2。其中r是圓周運(yùn)動(dòng)的半徑。切向加速度大小就是加速度向量在速度向量上的投影值，由向量的運(yùn)算規(guī)則可以求出切向加速度大小為：（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年都呈增長(zhǎng)趨勢(shì)，可能是隨著觀測(cè)年份的增加，古塔的年齡增長(zhǎng)，故彎曲的程度會(huì)隨之提高。這充分證明隨著年份的增加，塔身離1986年的基準(zhǔn)越大，且傾斜的幅度使先升高后減小，因此若塔繼續(xù)保持這種勢(shì)態(tài)塔就不會(huì)倒塌。

四、模型評(píng)價(jià)

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

參考文獻(xiàn)：

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

[2]薛毅.數(shù)學(xué)建?；A(chǔ).北京：北京工業(yè)大學(xué)出版社，2005.

[3]張德豐.MATLAB數(shù)值分析與應(yīng)用.北京：國防工業(yè)出版社，2009.

[4]葉其孝.大學(xué)生數(shù)學(xué)建模競(jìng)賽輔導(dǎo)材料[M].長(zhǎng)沙：湖南教育出版社，2002.

[5]張序.測(cè)量學(xué).南京：東南大學(xué)出版社，2007.

基金項(xiàng)目：中山職業(yè)技術(shù)學(xué)院教研教改課題《與電子類專業(yè)教學(xué)相結(jié)合的〈電路數(shù)學(xué)〉探究型項(xiàng)目式教學(xué)的研究》（編號(hào)：JYC1110），主持人：歐冰。