2013數(shù)學(xué)建模C題古塔問題探究

1、古塔的變形監(jiān)測(cè)與預(yù)測(cè)摘要：本文針對(duì)古塔長(zhǎng)時(shí)間承受自重和外力作用引起的組合變形問題，采用數(shù)據(jù)處理、幾何分析及曲線擬合的方法，利用matlab、lingo等數(shù)學(xué)軟件編程、計(jì)算，給出了確定古塔各層中心位置的通用方法，根據(jù)所給的各次測(cè)量數(shù)據(jù)計(jì)算出古塔各層的中心坐標(biāo)（見表9），分析了該塔的傾斜、彎曲、扭曲等變形情況并對(duì)變形趨勢(shì)做了預(yù)測(cè)。得出如下結(jié)論：1、隨著時(shí)間的增長(zhǎng)，古塔的傾斜弧度與時(shí)間呈二次拋物線分布，總的變形是先減小后增大趨勢(shì)。2、隨著時(shí)間的增長(zhǎng)，古塔的彎曲曲率與時(shí)間呈三次函數(shù)分布，曲率隨時(shí)間先增大，后減小，再增大的趨勢(shì)。3、隨著時(shí)間的增長(zhǎng)，古塔的扭曲角弧度與時(shí)間呈二次拋物線分布，總的扭曲變形是先

2、減小后增大趨勢(shì)。關(guān)鍵詞：中心坐標(biāo) 投影 組合變形 傾斜 彎曲 扭曲 1、 問題重述由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生各種變形，諸如傾斜、彎曲、扭曲等。為保護(hù)古塔，防止其變形而導(dǎo)致破壞，文物部門需適時(shí)對(duì)古塔進(jìn)行觀測(cè)，通過測(cè)量所得數(shù)據(jù)可以了解各種變形量，以制定必要的保護(hù)措施。某古塔已有上千年歷史，是我國(guó)重點(diǎn)保護(hù)文物。管理部門委托測(cè)繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)。請(qǐng)根據(jù)提供的4次觀測(cè)數(shù)據(jù)，討論以下問題：1. 列表給出各次測(cè)量的古塔各層中心坐標(biāo)，給出確定古塔各層中心位置的通用方法。2. 分析該

3、塔傾斜、彎曲、扭曲等變形情況。3. 分析該塔的變形趨勢(shì)。 二、模型假設(shè)1、 假設(shè)塔的每層是一個(gè)平面，不考慮樓層厚度；2、 假設(shè)古塔在2011年3月前沒有進(jìn)行過大的維修工程，如地基處理；古塔加固；3、假設(shè)觀測(cè)過程中測(cè)的數(shù)據(jù)精確度高，誤差很小，可以忽略不計(jì)；4、假設(shè)把各種外界因素對(duì)塔的影響看做一個(gè)整體來研究；5、假設(shè)古塔在受到外界因素干擾時(shí)，未受到嚴(yán)重破壞；三、符號(hào)說明：第i層觀測(cè)點(diǎn)的橫坐標(biāo)；：第i層觀測(cè)點(diǎn)的縱坐標(biāo)；：第i層觀測(cè)點(diǎn)的豎坐標(biāo)；G：表示的各層中心坐標(biāo)；：第1層塔中心塔尖中心連線的傾斜角；：塔尖到第1層豎坐標(biāo)差；：曲線彎曲曲率；：扭轉(zhuǎn)角；四、模型準(zhǔn)備 觀察所給的數(shù)據(jù)，發(fā)現(xiàn) 1986年和1

4、996年的測(cè)量數(shù)據(jù)中缺失了第13層觀測(cè)點(diǎn)5的數(shù)據(jù)，對(duì)此我們觀察當(dāng)年已有的13層觀測(cè)點(diǎn)5數(shù)據(jù)，計(jì)算各觀測(cè)點(diǎn)的平均坐標(biāo)，發(fā)現(xiàn)觀測(cè)值與平均值之差在允許范圍內(nèi)，因此我們選取平均值作為缺失數(shù)據(jù)的補(bǔ)充，1986年觀測(cè)點(diǎn)5補(bǔ)充數(shù)據(jù)（566.195,520.132,52.795），1996年觀測(cè)點(diǎn)5補(bǔ)充數(shù)據(jù)（566.2314,520.1482.52.731）。2009年和2011年觀測(cè)點(diǎn)只有1個(gè)，我們只能將此觀測(cè)數(shù)據(jù)作為塔尖的中心坐標(biāo)。五、模型建立與求解（一）問題一 古塔有13個(gè)塔層和1個(gè)塔尖，通過對(duì)數(shù)據(jù)的分析以及matlab作圖（如圖1）可以知道古塔1到13層的結(jié)構(gòu)相似，而塔尖和塔層觀測(cè)點(diǎn)數(shù)目不同，而且結(jié)構(gòu)

5、各異，由此我們分別討論塔層和塔尖，找到它們的中心。圖11 第1層到13層1.1 問題分析每層塔的8個(gè)觀測(cè)點(diǎn)的z不全相同，它們構(gòu)成一個(gè)不規(guī)則的曲面。首先，將8個(gè)觀測(cè)點(diǎn)投影到xoy面，投影點(diǎn)分別為。利用matlab擬合一個(gè)斜平面，滿足：8個(gè)觀測(cè)點(diǎn)與斜平面的交點(diǎn)與對(duì)應(yīng)觀測(cè)點(diǎn)的距離之和最小。依次連接這8個(gè)交點(diǎn)，所得圖形的中心就是要找的塔層中心。再將xoy面的投影點(diǎn)依次連接，利用lingo軟件擬合出一個(gè)到8個(gè)頂點(diǎn)的距離之和最小的圓，圓心G可以近似看做該八邊形的中心。最后，過圓心作平行于z軸的直線，它與斜平面的交點(diǎn)就是我們所要找的塔層的中心點(diǎn)。1.2模型建立設(shè)擬合的斜平面方程為： （1）建立模型： min

6、 （）設(shè)xoy面內(nèi)擬合的圓的方程為： （2） 建立模型：min 1.3模型求解根據(jù)已知數(shù)據(jù)，利用matlab擬合可以算出A，B，C，得到斜平面的方程，做出圖像（如圖2）；利用lingo軟件擬合，可以得到圓的方程，如圖3。將圓心坐標(biāo)，代入斜平面的方程可解出，即得到中心坐標(biāo)的通式 圖2圖3將四年測(cè)量的數(shù)據(jù)代入之后得到1到13層的中心坐標(biāo)（如表1）表1年份19861996坐標(biāo)層數(shù)XYZXYZ1566.665522.7091.7874566.6652522.70871.7832566.7225522.67147.3203566.7234522.67047.31463566.7787522.634812

7、.7553566.7803522.633112.75084566.8227522.605817.7553566.8248522.603617.07525566.87522.574621.7206566.8727522.571821.71616566.9188522.544326.2351566.9222522.540826.22947566.9522522.527129.8369566.9561522.523229.83238566.9847522.510333.351566.9888522.50633.34559567.0176522.493436.8551567.0223522.48843

8、6.848410567.0476522.479240.1723567.0527522.473940.167811567.1014522.438844.441567.1071522.43344.435512567.1554522.398448.7118567.1614522.392248.707313567.202522.3852.8288567.2082522.373752.8243塔尖567.3785522.113455.1818567.3375522.266355.244年份20092011 坐標(biāo)層數(shù)XYZXYZ1566.7447522.70031.7645566.7448522.7002

9、1.76322566.7790522.67197.3089566.7792522.67177.29053566.8117522.645412.7322566.8121522.645712.72684566.8386522.623217.0697566.839522.622717.05195566.8666522.600221.7093566.8672522.599621.70386567.0103522.639626.2085566.8967522.7426.20037567.2446 522.791929.8176566.9164522.746629.81328566.9645522.712

10、833.3353566.9653522.71233.3329567.0128522.677236.8391567.0138522.676336.816710567.1112522.521340.1583567.0686522.622540.138111567.1909522.3743.4323567.192522.368944.424512567.2331522.330248.6994567.2344522.328948.683613567.2817522.284652.818567.283522.283252.8127塔尖567.336522.214855.091567.3375522.21

11、3555.0871.4模型分析與評(píng)價(jià)該模型的建立，采用投影、擬合的方法將不規(guī)則轉(zhuǎn)化為規(guī)則，將空間轉(zhuǎn)化到平面，近似找到塔層的中心坐標(biāo)，由于第1到13層的結(jié)構(gòu)相似，因此方法具有通用性。即使1986年和1996年的測(cè)量數(shù)據(jù)中缺失了第13層觀測(cè)點(diǎn)5的數(shù)據(jù)，但是我們根據(jù)其他各層的數(shù)據(jù)猜測(cè)得到，不影響方法的使用。但是也因?yàn)椴聹y(cè)和兩次擬合，塔層中心坐標(biāo)的計(jì)算產(chǎn)生誤差。2 塔尖2.1 問題分析考慮到塔尖的橫截面面積相對(duì)塔層橫截面較小，且已知的測(cè)點(diǎn)數(shù)據(jù)表明，塔尖各測(cè)點(diǎn)變形相對(duì)較小，因此我們采用求平均值的方法來計(jì)算塔尖的中心坐標(biāo)。2.2 模型建立將塔尖的4個(gè)觀測(cè)點(diǎn)坐標(biāo)設(shè)為，中心坐標(biāo)的計(jì)算式為： 2.3 模型求解將

12、四次觀測(cè)數(shù)據(jù)分別代入平均值計(jì)算的中心坐標(biāo)公式： 算出塔尖的中心坐標(biāo)（如表2）表2年份19861996坐標(biāo)層數(shù)XYZXYZ塔尖567.3785522.113455.1818567.3375522.266355.244年份20092011 坐標(biāo)層數(shù)XYZXYZ塔尖567.336522.214855.091567.3375522.213555.0872.4 模型分析與評(píng)價(jià)該模型可以簡(jiǎn)單、方便地計(jì)算出塔尖的中心坐標(biāo)，但是2009年和2011年觀測(cè)點(diǎn)只有1個(gè)，我們只能將此觀測(cè)數(shù)據(jù)作為塔尖的中心坐標(biāo)。而且平均值計(jì)算的方法有它的缺點(diǎn)，觀測(cè)點(diǎn)的數(shù)目越少，結(jié)果的誤差就會(huì)越大。這是該模型有待改進(jìn)的地方。（二）問題

13、二1 問題分析在古塔剛修建成功時(shí)，不考慮修建過程中的誤差，各層中心的連線（中性軸）應(yīng)該是垂直于底平面的一條垂線。但是由于長(zhǎng)時(shí)間承受自重、氣溫、風(fēng)力等各種作用，偶然還要受地震、颶風(fēng)的影響，古塔會(huì)產(chǎn)生組合變形，包括傾斜、彎曲、扭曲等，此時(shí)其中性軸也會(huì)發(fā)生相應(yīng)的變形（見圖4）。因此，我們可以通過對(duì)古塔中性軸的變形來討論古塔的變形，分別單獨(dú)討論傾斜、彎曲和扭曲。圖4（中心點(diǎn)連線的圖）1.1 傾斜變形1.1.1問題分析塔的傾斜變形是一種塑性變形，是塔基固定，塔身和塔尖產(chǎn)生傾斜的一種變形。這種變形反映了整座塔的變形，因此我們考慮塔尖中心相對(duì)于第1層中心的傾斜。連接第1層的中心和塔尖的中心，將其投影到過且平

14、行于z軸的直線上，的投影點(diǎn)到的距離記作，到的距離記作，傾斜角。（如圖5） 圖51.1.2模型建立1.1.3 模型求解將第一題中計(jì)算出來的數(shù)據(jù)帶到以上公式，分別計(jì)算出四年古塔的傾斜角度（表3）表3年份1986199620092011傾斜位移值0.92940.80480.76510.7669傾斜度0.01680.01460.01390.0139傾斜角59'49.99''51'44.87''49'19.17''49'26.28''1.1.4模型分析與評(píng)價(jià)根據(jù)上述表格所得的數(shù)據(jù)可知，古塔的傾斜角度是很微小的

15、。在古塔開始發(fā)生傾斜時(shí)，由于地基基礎(chǔ)下的地質(zhì)分布不同及其他外界不利因素，導(dǎo)致古塔剛開始是向一邊傾斜。但是隨著時(shí)間的變化，以及各種作用下，古塔開始向其另一邊發(fā)生傾斜，從而抑制塔向一邊傾斜的趨勢(shì)。由二力平衡可知，兩邊對(duì)塔產(chǎn)生力的作用，這樣使塔處于相對(duì)平衡的狀態(tài)。由此我們可以從上述計(jì)算出的結(jié)果看出，在沒有大的事故作用下，塔在一定的范圍時(shí)間內(nèi)不會(huì)被外界作用所毀壞。2.1 彎曲變形2.1.1問題分析曲線的彎曲程度可以用曲率來度量，依次連接塔層中心點(diǎn)，得到的是空間折線段，將它擬合成一條曲線，參數(shù)方程為：。代入曲率公式可以計(jì)算出各點(diǎn)的曲率，由于曲線各點(diǎn)的彎曲程度不同，我們用曲線的最大曲率來度量古塔的彎曲變形

16、程度。2.1.2模型建立曲線方程：max 2.1.3模型求解 利用第一題計(jì)算的中心坐標(biāo)，運(yùn)用matlab軟件擬合出曲線（見圖6）的方程，計(jì)算出參數(shù)方程的一階導(dǎo)數(shù)、二階導(dǎo)數(shù)代入曲率公式，得到再利用matlab對(duì)k求導(dǎo)數(shù)，計(jì)算出當(dāng)時(shí)，k取得最大值。圖6將第一題求得的中心坐標(biāo)帶入，按照上述方法計(jì)算，可以分別計(jì)算不同年古塔的最大曲率（見表4）。表4年份1986199620092011曲率1.0303×10-97.8736×10-103.6160×10-78.003×10-73.1扭曲變形3.1.1問題分析 扭曲變形是指各塔層橫截面繞軸線的相對(duì)轉(zhuǎn)動(dòng)，轉(zhuǎn)動(dòng)后與其原來

17、的位置形成一定的角度，稱這個(gè)角為相對(duì)扭轉(zhuǎn)角?？紤]古塔的扭曲是繞著中性軸的轉(zhuǎn)動(dòng)，而每層扭轉(zhuǎn)的程度不一樣，因此分層討論。假設(shè)第1層為基準(zhǔn)面，不發(fā)生扭曲，考慮其他各層相對(duì)于它繞著中性軸的扭曲程度。而塔尖的數(shù)據(jù)不全，且近似看作在中性軸上，因此只考慮第2層到第13層的扭曲。每一個(gè)塔層有8個(gè)觀測(cè)點(diǎn)，任選一個(gè)作為扭曲變形的研究對(duì)象，不妨假設(shè)觀測(cè)點(diǎn)1。將第1層的觀測(cè)點(diǎn)1投影到第i層中心點(diǎn)所在的水平面，投影點(diǎn)記作，分別連接與以及第i層觀測(cè)點(diǎn)1（）與，它們之間的夾角就是要找的相對(duì)扭轉(zhuǎn)角（如圖7）。圖73.1.2模型建立已知坐標(biāo)，, ，利用余弦定理列出計(jì)算式(i =2、3.13) 。 ， 3.1.3模型求解將題中所

18、給的數(shù)據(jù)以及第一題計(jì)算的中心坐標(biāo)代入上式，利用matlab可以計(jì)算出每一層相對(duì)于第一層的扭轉(zhuǎn)角（見表5）。表5 1986年的值1996年的值2009年的值2011年的值 0 0.0469 0.0945 0.1334 0.1766 0.2205 0.2810 0.3444 0.4131 0.4832 0.5666 0.6581 0.7597 0 0.0470 0.0962 0.1334 0.1768 0.2223 0.2811 0.3451 0.4133 0.4852 0.5673 0.6616 0.7613 0 0.0529 0.1054 0.1498 0.2012 0.2950 0.3227

19、 0.4093 0.4807 0.5412 0.6106 0.6957 0.7886 0 0.0493 0.1053 0.1497 0.2011 0.3021 0.3460 0.4134 0.4806 0.5461 0.6103 0.6954 0.7883 （三）問題三3.1問題分析 變形趨勢(shì)就是一種慣性，沒有外力不會(huì)輕易改變狀態(tài)。由于古塔的變形是多種影響的組合變形，如傾斜、彎曲、扭曲等等。這些影響來自各個(gè)方面，不便于綜合確定各種影響的程度，這樣就造成了無法整體去研究它的變形趨勢(shì)。因此我們還是分開考慮，分別考慮傾斜、彎曲、扭曲對(duì)塔的變形趨勢(shì)。這些變形趨勢(shì)都是建立在問題二的基礎(chǔ)之上，模型相同，因

20、此在以下變形趨勢(shì)中同樣適用。3.1.1傾斜變形趨勢(shì)分析在問題二傾斜變形模型求解中，我們計(jì)算出了塔的整體傾斜角（見表3），我們?cè)儆脮r(shí)間（單位是年）與傾斜角（單位是弧度）建立平面直角坐標(biāo)系，以O(shè)為坐標(biāo)原點(diǎn)，時(shí)間（單位是年）為X軸,傾斜角（單位是弧度）為Y軸，坐標(biāo)系表示傾斜角與時(shí)間變化的關(guān)系。所以我們可以得出四個(gè)坐標(biāo)點(diǎn)，再通過Excel處理擬合曲線作圖（見表6、圖8）。表6年份1986199620092011傾斜角(弧度)0.99720.86250.8220.824圖8根據(jù)上述圖表預(yù)測(cè)，隨著時(shí)間的增長(zhǎng)，古塔的傾斜弧度與時(shí)間呈二次拋物線分布，當(dāng)（年）時(shí)塔的傾斜角最小，即傾斜程度最小，塔的變形最小?？偟?/p>

21、變形是先減小后增大趨勢(shì)。 彎曲變形趨勢(shì)分析在問題二彎曲變形模型中，我們求解出了曲線的最大曲率，曲線的最大曲率能夠幫助我們分析曲線的彎曲程度?；诖饲闆r，我們?cè)趩栴}二的彎曲變形模型基礎(chǔ)上去分析彎曲變形趨勢(shì)。首先，建立平面直角坐標(biāo)系，O為坐標(biāo)原點(diǎn)，OX軸代表時(shí)間（年），OY軸代表曲率，此坐標(biāo)系是表示曲率與時(shí)間的變化關(guān)系。 所以我們可以得出四個(gè)坐標(biāo)點(diǎn)，再通過Excel處理擬合曲線作圖（見表7、圖9）。 表7 時(shí)間1986199620092011曲率（10-10）1.03037.873636168003圖9根據(jù)上述圖表預(yù)測(cè)，隨著時(shí)間的增長(zhǎng)，古塔的彎曲曲率與時(shí)間呈三次函數(shù)分布，曲率隨時(shí)間先增大，后減小，

22、再增大的趨勢(shì)。扭曲變形趨勢(shì)分析在問題二扭曲變形模型中，已經(jīng)求解出了各層相對(duì)于第1層的扭曲角（見表5），觀察表的數(shù)據(jù)，可以找出一規(guī)律：隨著塔層的增高，扭曲角增大，即扭曲程度增大。最大扭曲角來最能夠代表整座塔的變形，因此我們還是選取它來研究扭曲變形趨勢(shì)。同樣建立平面直角坐標(biāo)系，O為坐標(biāo)原點(diǎn)，OX軸代表時(shí)間（年），OY軸代表最大扭曲角（弧度）。此坐標(biāo)系是表示扭曲角隨時(shí)間的變化關(guān)系。所以我們可以得出四個(gè)坐標(biāo)點(diǎn)，再通過Excel處理擬合曲線作圖（見表8、）。 表8年份1986199620092011扭曲角（弧度）0.75970.76130.78860.7883圖10根據(jù)上述圖表預(yù)測(cè)，隨著時(shí)間的增長(zhǎng)，古塔

23、的扭曲角弧度與時(shí)間呈二次拋物線分布，當(dāng)（年）時(shí)塔的扭曲角最小，即是扭曲程度最小。總的扭曲變形是先減小后增大趨勢(shì)。六、模型推廣與評(píng)價(jià)一、優(yōu)點(diǎn)：本文通過對(duì)數(shù)據(jù)的處理、分析，利用MATLAB和LINGO軟件對(duì)古塔每年每層的測(cè)量點(diǎn)數(shù)據(jù)進(jìn)行擬合，有效地將不規(guī)則立體轉(zhuǎn)化為規(guī)則的平面圖形，利用規(guī)則平面圖形的中心來找到塔層的中心。這樣的處理很方便、明了，便于實(shí)現(xiàn)。此模型可用于大多數(shù)文物建筑，通用性較強(qiáng)。二、缺點(diǎn)：在數(shù)據(jù)的處理中，由于對(duì)數(shù)據(jù)采用了“四舍五入”方法，并且預(yù)測(cè)補(bǔ)充缺失數(shù)據(jù)，多次擬合，從而使數(shù)據(jù)誤差增大，降低了其精確性。 7、 參考文獻(xiàn)1 趙靜 但琦，數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)第三版，北京：高等教育出版社，2

24、008年1月。2 王兵團(tuán)，數(shù)學(xué)建?；A(chǔ)，北京:清華大學(xué)出版社，2004年11月。3 姜啟源,謝金星,葉俊,數(shù)學(xué)模型第三版，北京：高等教育出版社,2003年。4 于潤(rùn)偉 朱曉慧，MTLAB基礎(chǔ)及應(yīng)用第二版，北京：機(jī)械工業(yè)出版社，2011年1月。5 同濟(jì)大學(xué)航空航天與力學(xué)學(xué)院基礎(chǔ)力學(xué)教學(xué)研究部，材料力學(xué)，同濟(jì)大學(xué)出版社，2008年8月。6 同濟(jì)大學(xué)數(shù)學(xué)系，高等數(shù)學(xué)第六版上下冊(cè)，高等教育出版社，2007年4月。7 廖武鵬 鄧俊曄 王丹，管道切片的三維重建，工程數(shù)學(xué)學(xué)報(bào)，第19卷：22-34，2002年2月。8 梁?？?，古塔變形測(cè)量方法探討，城市勘測(cè)，第3期：114-118頁,2011年6月。9 胡志

25、曉，古塔傾斜觀測(cè)和數(shù)據(jù)分析，江蘇建筑，第6期：34-44頁，2011年6月。10 盧俊龍 張萌 姚謙峰，磚石古塔糾偏工程的有限元分析八、附錄1、用matlab實(shí)現(xiàn)空間斜平面擬合： x=565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5'y=528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356'z=1.792 1.818 1.783 1.769 1.772 1.77 1.794 1.801'scatter3(x,y,z,'

26、filled')hold on X = ones(8,1) x y; b = regress(z,X) ； xfit = min(x):0.1:max(x);yfit = min(y):0.1:max(y);XFIT,YFIT= meshgrid (xfit,yfit); ZFIT = b(1) + b(2) * XFIT + b(3) * YFIT;mesh (XFIT,YFIT,ZFIT) x=566.6650;y=522.7090;Z=b(1)+b(2)*x+b(3)*y2、運(yùn)用lingo軟件對(duì)數(shù)據(jù)處理的程序：model:sets:dianji/1.8/:x,y,d;endset

27、sdata: x=565.454 562.058 561.39 563.782 567.941 571.255 571.938 569.5;y=528.012 525.544 521.447 518.108 517.407 519.857 523.953 527.356;enddatamin=sum(dianji(i):(sqrt(x(i)-x0)2+(y(i)-y0)2)-r)2);end通過此程序可以到：1986、1996、2009、2011年每層的中心坐標(biāo)，下面是1986年所得的數(shù)據(jù)： 1、 Variable Value Reduced Cost X0 566.6650 0.000000

28、 Y0 522.7090 0.000000 R 5.427619 0.000000 2、 Variable Value Reduced Cost X0 566.7225 0.000000 Y0 522.6714 0.000000 R 5.226233 0.000000 3、 Variable Value Reduced Cost X0 566.7787 0.000000 Y0 522.6348 0.000000 R 5.028751 0.000000 4、 Variable Value Reduced Cost X0 566.8227 0.000000 Y0 522.6058 0.000000

29、 R 4.872247 0.000000 5、 Variable Value Reduced Cost X0 566.8700 0.000000 Y0 522.5746 0.000000 R 4.704324 0.000000 6、 Variable Value Reduced Cost X0 566.9188 0.000000 Y0 522.5443 0.000000 R 4.541165 0.000000 7、 Variable Value Reduced Cost X0 566.9522 0.000000 Y0 522.5271 0.000000 R 4.272584 0.000000

30、8、 Variable Value Reduced Cost X0 566.9847 0.000000 Y0 522.5103 0.000000 R 4.011286 0.000000 9、 Variable Value Reduced Cost X0 567.0176 0.000000 Y0 522.4934 0.000000 R 3.750459 0.000000 10、 Variable Value Reduced Cost X0 567.0476 0.000000 Y0 522.4792 0.000000 R 3.503796 0.000000 11、 Variable Value R

31、educed Cost X0 567.1014 0.000000 Y0 522.4388 0.000000 R 3.276693 0.000000 12、 Variable Value Reduced Cost X0 567.1554 0.000000 Y0 522.3984 0.000000 R 3.049347 0.000000 13、 Variable Value Reduced Cost X0 567.2020 0.000000 Y0 522.3800 0.000000 R 2.819595 0.000000按照同樣的方法我們可以的到1996，2009,2011年古塔的每層中心坐標(biāo)（見

32、表9）：年份19861996坐標(biāo)層數(shù)XYZXYZ1566.665522.7091.7874566.6652522.70871.7832566.7225522.67147.3203566.7234522.67047.31463566.7787522.634812.7553566.7803522.633112.75084566.8227522.605817.7553566.8248522.603617.07525566.87522.574621.7206566.8727522.571821.71616566.9188522.544326.2351566.9222522.540826.2294756

33、6.9522522.527129.8369566.9561522.523229.83238566.9847522.510333.351566.9888522.50633.34559567.0176522.493436.8551567.0223522.488436.848410567.0476522.479240.1723567.0527522.473940.167811567.1014522.438844.441567.1071522.43344.435512567.1554522.398448.7118567.1614522.392248.707313567.202522.3852.8288

34、567.2082522.373752.8243塔尖567.3785522.113455.1818567.3375522.266355.244表9年份20092011 坐標(biāo)層數(shù)XYZXYZ1566.7447522.70031.7645566.7448522.70021.76322566.7790522.67197.3089566.7792522.67177.29053566.8117522.645412.7322566.8121522.645712.72684566.8386522.623217.0697566.839522.622717.05195566.8666522.600221.7093

35、566.8672522.599621.70386567.0103522.639626.2085566.8967522.7426.20037567.2446 522.791929.8176566.9164522.746629.81328566.9645522.712833.3353566.9653522.71233.3329567.0128522.677236.8391567.0138522.676336.816710567.1112522.521340.1583567.0686522.622540.138111567.1909522.3743.4323567.192522.368944.424

36、512567.2331522.330248.6994567.2344522.328948.683613567.2817522.284652.818567.283522.283252.8127塔尖567.336522.214855.091567.3375522.213555.0873、根據(jù)空間參數(shù)方程： 通過matlab軟件對(duì)上述數(shù)據(jù)進(jìn)行擬合： x=data(:,1);>> z=data(:,3);>> polyfit(z,x,2);>> y=data(:,2);>> z=data(:,3);>> polyfit(z,y,2);得到具體

37、的空間曲線方程（以1986年的為例）：通過下面這個(gè)程序，畫出1986年的空間曲線：>> t=0:100;>> x=1.556978433478483e-05*t.2+0.0095*t+566.6518;>> y=-3.775091994799030e-06*t.2 -0.0061*t+ 522.7157;>> z=t;>> plot3(x,y,z,'r-');然后用此方法分別對(duì)每年的數(shù)據(jù)進(jìn)行擬合，這樣我們就可以分別得到每年空間曲線方程和空間曲線。用diff命令對(duì)空間曲線方程求導(dǎo)：syms tdiff(1.556978433478483e-05*t2+0.0095*t+566.6518)diff(1.556978433478483e-05*t2+0.0095*t+566.6518,2)代入曲率公式，這樣便可求出此曲線的曲率。4、運(yùn)用matlab軟件求最優(yōu)半徑程序：function R,A,B=circ(x,y,N)x1 = 0; x2 = 0; x3 = 0; y1 = 0; y2 = 0; y3 = 0; x1y1 = 0; x1y2 = 0; x2y1 = 0; for i = 1 : N x1 = x1 + x(i); x2