1、2013高教社杯全國大學生數學建模競賽承 諾 書我們仔細閱讀了全國大學生數學建模競賽章程和全國大學生數學建模競賽參賽規則（以下簡稱為“競賽章程和參賽規則”，可從全國大學生數學建模競賽網站下載）。我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網上咨詢等）與隊外的任何人（包括指導教師）研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽章程和參賽規則的，如果引用別人的成果或其他公開的資料（包括網上查到的資料），必須按照規定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽章程和參賽規則，以保證競賽的公正、公平性。如有違反競賽章程和參賽

2、規則的行為，我們將受到嚴肅處理。我們授權全國大學生數學建模競賽組委會，可將我們的論文以任何形式進行公開展示（包括進行網上公示，在書籍、期刊和其他媒體進行正式或非正式發表等）。我們參賽選擇的題號是（從A/B/C/D中選擇一項填寫）： C 我們的參賽報名號為（如果賽區設置報名號的話）： 5339 所屬學校（請填寫完整的全名）： 參賽隊員 (打印并簽名) ：1. 2. 3. 指導教師或指導教師組負責人 (打印并簽名)： （論文紙質版與電子版中的以上信息必須一致，只是電子版中無需簽名。以上內容請仔細核對，提交后將不再允許做任何修改。如填寫錯誤，論文可能被取消評獎資格。） 日期： 2013 年 09 月

3、 16 日賽區評閱編號（由賽區組委會評閱前進行編號）：2013高教社杯全國大學生數學建模競賽編 號 專 用 頁賽區評閱編號（由賽區組委會評閱前進行編號）：賽區評閱記錄（可供賽區評閱時使用）：評閱人評分備注全國統一編號（由賽區組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：古塔的變形數學模型摘要:本文是研究關于古塔變形類型以及變形分析的模型，用Matlab畫出古塔的三維結構可以看出它是近似于正八邊形的形狀。因此，問題一我們用每層各個測量點坐標的平均值作為塔每層的中心坐標，再用中心坐標的三個坐標值分別對時間t做回歸來得到確定古塔各層中心位置的通用方法。對于問題二，我們分別研究該

4、塔傾斜、彎曲、扭曲等變形情況，通過建立數學模型來確定變形的程度。首先，用各層中心點坐標的z坐標值與其相應點的x,y坐標值做多元線性回歸。然后得到的回歸方程所表示的回歸平面與z軸正方向的夾角就可以表示古塔的傾斜程度大小。最后根據各層中心的分布和變化趨勢方向，確定古塔的傾斜方向。用古塔各層中心點進行平面擬合，從效果上觀察，較為精確地反映了實例中的問題，由此也說明了我們所建模型的合理性。古塔的傾斜變形必然會導致在同一層中，測點存在高程的絕對差h，如果古塔只存在傾斜變形的話，每層的h值會相等；如果古塔存在傾斜變形的同時也存在彎曲變形的話，則每層的h值會發生改變。所以相鄰兩層的高程絕對差的變化量，表示古

5、塔每層彎曲程度大小。再根據每層出現高程絕對差h的兩個測量點的連線，確定每層彎曲方向。古塔的扭曲變形，首先每層選取兩對相同的對測量點，并做連線。然后通過每層對測量點的連線，分別與第一層相同對測量點的連線所成的角度的平均值來衡量古塔的扭曲情況。對于該塔的變形趨勢的研究，將傾斜指標、彎曲指標、扭曲指標對時間的回歸。再用得到的回歸方程預測未來幾年的數據，結合用excel畫出的圖來預測古塔在未來時間里的變形趨勢。關鍵字： 線性回歸 變化趨勢 擬合 預測1、 問題重述由于長時間承受自重、氣溫、風力等各種作用，偶然還要受地震、颶風的影響，古塔會產生各種變形，諸如傾斜、彎曲、扭曲等。為保護古塔，文物部門需適時

6、對古塔進行觀測，了解各種變形量，以制定必要的保護措施。某古塔已有上千年歷史，是我國重點保護文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進行了4次觀測。請你們根據附件1提供的4次觀測數據，討論以下問題：1. 給出確定古塔各層中心位置的通用方法，并列表給出各次測量的古塔各層中心坐標。2. 分析該塔傾斜、彎曲、扭曲等變形情況。3. 分析該塔的變形趨勢。2、 問題分析（1） 、對問題一的分析 問題一中確定古塔各層中心位置的通用方法。因為古塔各層為近似正八邊形，根據正八邊形圖形特征，可以用每次測量時，古塔各層測量點坐標的平均值作為各層中心點坐標。然

7、后將各層中心點坐標對時間回歸，可得到各層中心點坐標對時間的回歸方程。根據方程就可以確任意時間各層中心點坐標。（2） 、對問題二的分析問題二要求我們確定塔的傾斜、彎曲、扭曲等變形情況。為了簡化模型，我們分別對古塔的傾斜、彎曲、扭曲進行討論。對于傾斜，首先根據不同年份，各層中心點坐標的z坐標值與其相應點的x,y坐標值做多元線性回歸。然后得到的回歸方程所表示的回歸平面與z軸正方向的夾角就可以表示古塔的傾斜程度大小。最后根據各層中心的分布和變化趨勢方向，確定古塔的傾斜方向。對于古塔的彎曲，首先求出每層高程絕對差，然后相鄰兩層的高程絕對差的變化量，表示古塔每層彎曲程度大小。再根據每層出現高程絕對差的兩個

8、測量點的連線，確定每層彎曲方向。對于古塔的扭曲變形，首先每層選取兩對相同的對測量點，并做連線。然后通過每層對測量點的連線，分別與第一層相同對測量點的連線所成的角度的平均值來衡量古塔的扭曲情況。（3） 、對問題三的分析問題三要求我們分析該塔的變形趨勢，這個問題屬于預測的數學問題。對于這個問題我們一般用回歸的方法來求解，得出傾斜指標、彎曲指標、扭曲指標對時間的回歸方程，并作出各自的圖像，觀察趨勢。3、 模型假設1. 假設古塔只存在傾斜，彎曲，扭曲的三種變形情況；2. 假設在1986年到2011年沒有對古塔進行人為的保護，如加固或修補；3. 忽略1986年與1996年觀測的第13層第5個測量點所少數

9、據；4. 假設古塔的變形是連續的；四、符號說明觀測點古塔層數時間，并以1986年為第一年，即t=1古塔同一層測點的最大高程差古塔相鄰兩層h值的差兩對角線的夾角古塔的傾斜角五、模型的建立與求解5.1問題一的求解：為觀察同層各觀測點的大概位置，做出1986年古塔內同層觀測點連線的俯視圖進行分析，做出下圖：圖1-1 圖1-1是通過1986年每層各測量點的坐標點連起來的（用CAD制）圖，每層所測的點相交構成一個多邊形，得到每層的近似平面圖，可以近似地把每層當作正八邊形。根據正八邊形圖形特征，古塔各層測量點坐標（，）的平均值作為各層中心點坐標。即：， ， 算出的各層中心坐標如下：表1-1.所測年數各層中

10、心坐標表1986年1996年樓層iXYZ樓層iXYZ1566.6648522.71051.78741566.665522.71021.71022566.7196522.66847.32032566.7205522.66747.31463566.7735522.627312.75533566.7751522.625612.75084566.8161522.594417.07834566.8183522.592217.07515566.8621522.559121.72055566.8649522.556321.7166566.9084522.524426.23516566.9118522.521

11、26.22957566.9468522.508129.83697566.9506522.504229.83238566.9843522.492433.35098566.9884522.488133.34549567.0218522.476436.85499567.0265522.471436.848310567.0569522.462440.172110567.062522.457240.167611567.1045522.42344.408811567.1102522.417344.435412567.1518522.383648.711912567.1578522.377548.70741

12、3567.085522.740352.834313567.0912522.73452.83塔尖567.2473522.243855.1233塔尖567.2544522.236755.119752009年2011年樓層iXYZ樓層iXYZ1566.7268522.70151.76451566.727522.70141.763252566.764522.66937.3092566.764522.6697.29053566.8001522.638412.73233566.8004522.638712.72694566.8293522.613217.06984566.8297522.612717.05

13、25566.8604522.586621.70945566.861522.58621.70396566.947522.534226.2116566.9478522.533526.20457566.9792522.512329.82467566.98522.511529.8178567.0305522.479733.33998567.0313522.478833.33669567.0816522.446636.84389567.0825522.445736.822310567.137522.393740.161110567.1381522.392640.144111567.1799522.354

14、744.432611567.181522.353544.424912567.2225522.31648.699812567.2238522.314748.683913567.2712522.271552.818413567.2725522.270152.8131塔尖567.336522.214855.091塔尖567.3375522.213555.087用上表得到的數據，把每層中心點的X，Y，Z坐標分別對時間t（設1986年為第一年，即t=1）做回歸，得到下表的一系列回歸方程，用以下的方程就能算出古塔任意一年任意一層的中心坐標,即為確定古塔各層中心位置的通用方法。表1-2.各層坐標與時間的回歸

15、表中心點坐標層數iXiYiZi1X=566.65164+0.00286tY=522.71231-0.00041tZ=1.76325-0.00024t2X=566.71052+0.00203tY=522.66783+0.00044tZ=7.322956-0.00093t3X=566.7685+0.001214tY=522.62418+0.00054tZ=12.75904-0.00115t4X=566.81415+0.00059tY=522.58954+0.00088tZ=12.75904-0.00115t5X=566.86355-0.00009tY=522.55215+0.00128tZ=21.

16、72178-0.0006t6X=566.90148+0.00176tY=522.52099+0.00047tZ=26.23885-0.00121t7X=566.94137+0.00147tY=522.50569+0.00022tZ=29.83865-0.00071t8X=566.97606+0.00210tY=522.49343-0.00056tZ=33.35139-0.00053t9X=567.01090+0.00272tY=522.48078-0.00134tZ=36.85757-0.00099t10X=567.04181+0.00366tY=522.47502-0.00313tZ=40.

17、17504-0.00089t11X=567.09066+0.00322tY=522.43534-0.00311tZ=44.41597-0.00061t12X=567.13906+0.00922tY=522.39564-0.00308tZ=48.71487-0.00091t13X=567.04774+0.00853tY=522.83776-0.02153tZ=52.83664-0.00082t14X=567.23116+0.00404tY=522.24712-0.00129tZ=55.12923-0.00155t5.2問題二的求解：問題二要求我們確定塔的傾斜、彎曲、扭曲等變形情況。為了簡化模型，

18、我們分別對古塔的傾斜、彎曲、扭曲進行討論。5.2.1傾斜變形求解：1) 傾斜度大小的求解 首先用第t年中,各層中心點坐標的對，做多元線性回歸。 = （1） 該回歸方程在空間直角坐標系中是一個平面，表示各層中心近似所處的平面。運用Excel軟件根據（1）式求解各年回歸方程。統計各年回歸方程系數可得下表：表2-1 回歸方程系數年份t19861-82920.6116.501332.34341199611-82194.2115.423532.12215200924-57953.592.5080710.58585201126-5715091.765449.853791運用Matlab做出2009年回歸方

19、程對應的回歸平面（過程見附錄1）。 圖2-1由圖2-1可以直觀地看出各層中心點貼近回歸平面，證明上面所建立模型的準確性。 上述方法所得的回歸平面與z軸正方向的夾角可以表示塔的傾斜角。原理解釋如下： 圖2-2 如圖，空間直角坐標系（由CAD制作）中有下列關系： 其中，平面ABC是同一年塔內各層中心點的回歸平面。AB垂直于OD，AB垂直于OC，即AB垂直于CD，即ODC為平面OAB與平面ABC夾角。所以平面ABC與z軸夾角為： 即角為回歸平面與z軸正方向的夾角可以表示塔的傾斜角。令=0，即方程為=0為AB所在直線方程，所以OD為點O到AB的距離，根據點到直線的距離公式可得：令x和y等于0，可以得出

20、 根據正切性質得： 根據反三角函數，可知： (2) 根據（2）式可以得算出古塔每年的傾斜角，列表如下：表2-2 塔的傾斜角度，單位：（°）年份1986199620092011塔的傾斜角0.47390.47820.61530.6208夾角的值可以表示古塔的傾斜程度大小。2) 傾斜方向的求解根據古塔各層中心點在水平面xoy中的投影的分布和變化趨勢，來確定古塔的傾斜方向。下面以2009年數據為例。用Matlab作2009年各測點與中心點的平面圖（過程見附錄2）： 87654321 圖2-3 由圖2-3可以看出，各層中心點都大致分布在第2，6個測量點的對角線上。再根據中心點投影位置隨樓層的增

21、加而自測量點2向測量點6移動。可以知道古塔的傾斜方向大致是沿測量點2向測量點6方向傾斜。5.2.2彎曲變形求解1) 彎曲程度大小的求解對于古塔的彎曲情況，我們通過每層平面傾斜的變化程度初步分析，然后再結合整棟古塔，得出古塔的大概外形，從局部到整體分析古塔的變形。首先，計算第t年，i層測量點高程的絕對差，為：它能直觀地反映在各層內最大傾斜程度，但絕對差不能全面的表現出彎曲的情況。而第t年i層到i+1層高程的絕對差的變化量，可以反映相鄰兩層的彎曲的大小程度。公式如下：運用excel計算與排列每年各層的測量點絕對差，得下表：表2-3每年各層的測量點絕對差，單位：m絕對差層數/年份1986199620

22、09201110.0490.0550.0790.07820.050.0490.0480.05730.050.0520.0760.08340.050.0460.0440.07350.050.0510.050.05160.1730.1660.1980.19770.1740.1670.160.15380.1740.180.1860.18590.1740.1670.1730.198100.2060.2050.2040.237110.1960.1950.1920.207120.1960.1950.1880.168130.1940.1930.1920.194通過上表中的絕對差計算第t年i層到i+1層高程的

23、絕對差的變化量。得出1986年，1996年，2009年，2011年各個值，得出下表： 表2-4相鄰兩層高程的絕對差的變化量表，單位：m層數1986年1996年2009年 2011年120.001-0.006-0.031-0.0212300.0030.0280.026340-0.006-0.032-0.014500.0050.006-0.022560.1230.1150.1480.146670.0010.001-0.038-0.0447800.0130.0260.032890-0.013-0.0130.0139100.0320.0380.0310.0391011-0.01-0.01-0.012-

24、0.03111200-0.004-0.0391213-0.002-0.0020.0040.026從所得各年值的統計可得出結論：在各年中都是56層間的值為最大，從1986年到2011年，值呈某部分增大，而且總體的總值都在變大。對統計數據進行分析：在56層間的傾斜程度變化很大，表現為56層間產生大的彎曲，在1986年到2011年塔的相鄰層之間的傾斜變化程度越來越大。表現為隨著時間的推移塔的彎曲程度越來越嚴重。的值主要在第五與第六層之間有較大值，其他值都較小。所以古塔主要彎曲的地方是在第5層與第6層。所以用第5層到第6層間的高程絕對差的增量的值來表示古塔的彎曲大小程度。由上表可以知在每年中，每層的高

25、程的絕對差都是較小的數，即表示該塔同層內高程起伏較小，基本處于同一水平面；在第5層至第6層高程的絕對差異常增大，且使在第6層以后每層的絕對差都穩定在較大數值；在2011年中，在相鄰兩層的高程的絕對差的數值變化較大。 圖2-4圖2-4是用Matlab作出古塔的三維結構圖（過程在附件2），可以看出古塔中間的地方比較彎曲，由此也說明了我們所建模型的合理性。2) 彎曲方向的求解列出在第5、6層各點高程從高到低的測量點序號： 表2-5 第5、6層各年各測量點的排序1986年1996年2009年2011年第5層222288811117777533666534565344481986年1996年2009年2

26、011年第6層22333322111144448888555577776666從上表可得出結論：每年的第5層樓到第6層樓，測量點6相對其他點降低的程度最大，即可知在第5層到第6層樓中古塔向測量點6方向彎曲。5.2.3扭曲變形求解對于古塔的扭曲變形，首先每層選取兩對相同的對測量點，并做連線。然后通過每層對測量點的連線，分別與第一層相同對測量點的連線所成的角度的平均值來衡量古塔的扭曲情況。由于古塔首層是與穩定的地基連接，所以首層的扭曲變形幾乎可以忽略不計，所以取每層相同對角線在平面xoy投影的夾角作為扭曲變形的量度：夾角越大，則扭曲變形越嚴重。具體原理如下：如圖2-5所示，以第1與第13層為例，首

27、先選取兩對測量點，分別為2-6、4-8。取第一層的2-6對角線為L1，第13層的2-6對角線為L2，L1與L2的夾角為。取第一層的4-8對角線為L3，第13層的4-8對角線為L4，L3與L4的夾角為。用與的平均值作為古塔第13層相對于第一層扭曲變形的量度。12345678L3L4 圖2-5 同理，其他層2-6測量點對角線與第一層2-6測量點對角線L1的夾角為，其他層4-8測量點對角線與第一層4-8測量點對角線L3的夾角為。它們的平均值反映其他層相對第一層的扭曲度。即： 計算結果如下表表2-7 每層相對于第一層扭曲的度，單位：（°）1986年1996年2009年2011年0.50485

28、60.5047130.3108580.3107641.0344411.0344410.6422510.6419261.4917171.4917170.9242250.9242252.0181542.0178871.2512371.2509742.5700912.5699332.0318552.0325022.8833382.8830062.3507992.3507993.2271833.2271832.9175442.9175443.6229583.6219163.5455453.5466064.0536174.0534273.5283223.5279533.9985913.9985913.5

29、016613.5016613.9416813.9412963.4748983.4748983.8800773.8802833.4496713.448519 得出結論：第2層到第10層隨著層數的增高的數值越來越大，第11層到第13層隨著層數的增高的數值越來越小。表現為在第2層到第10層的各層塔向逆時針扭曲，第11層到第13層的各層塔向順時針扭曲。5.3問題三的求解：5.3.1塔傾斜變形的趨勢：要確定塔的變形趨勢，對于傾斜變形，我們利用問題二中所求出的塔的傾斜角來對時間做回歸，得到，可以知道回歸方程為： 利用上述公式可以得到2011年之后每年的傾斜角，列表如下： 圖3-1 塔的傾斜角與時間的變化曲

30、線 由上圖所得出來的擬合曲線圖我們可以得出結論：隨著時間的增長，塔傾斜所成的傾斜角度值會越來越大。預測在2025年塔的傾斜角度值達到0.7º，2039年塔的傾斜角度值達到0.8°。5.3.2塔彎曲變形的趨勢： 討論該塔彎曲程度的變化趨勢，即可變為討論不同時間第5層到第6層間的高程絕對差的增量隨時間變化的變化。因為每年都是在第5層到第6層為最大彎曲，也為整個塔的主要彎曲的表現，所以可用對t進行回歸來說明塔隨時間的變化導致其彎曲程度的變化，對回歸方程求解；= 運用Excel軟件的LINEST函數求解該回歸方程系數為： =0.11461 =0.001186可以得出： （4） 用（

31、4）式在excel軟件中畫出彎曲程度的擬合曲線，如下： 圖3-2 塔的彎曲程度與時間的變化曲線預測在2011年以后，每年塔的彎曲程度都會增大。預測在2020年塔彎曲程度為0.15612，即為塔在第5層到第6層的測量點最大高程差的變化增加了0.15612米，在2050年他彎曲程度為0.17391，即為塔在第5層到第6層的測量點最大高程的變化增加了0.17391米。5.3.3塔扭曲的變形趨勢：由于每層相對于第一層都存在一定的扭曲，因此取各層對第一層的扭曲度數的平均值作為古塔的扭曲程度。用對時間t的回歸求回歸方程，這樣就得到古塔扭曲變形的趨勢?；貧w方程為： 運算得到與的值(過程見附錄.3)，則回歸方

32、程解析式為; 圖3-3由圖3-3可以知道用excel作出來的曲線的斜率為負值，因此古塔的扭曲程度隨時間的增加越來越小。六 、模型的評價與推廣6.1 模型的評價 本模型從實際出發，分析了多種情況，將古塔的變形分成傾斜、彎曲、扭曲三種類型來分別討論。先定性地分析古塔的變形，再定量解出古塔變形了多少。在通過結合回歸思想，將變形的數據轉化為回歸方程，建立了古塔變形趨勢的模型，使得對古塔的變形了解得更全面，能為古塔有效地制出保護措施提供依據。同時模型也存在著不足的地方，由于題目給出的數據不是很全面，所以很難算出古塔準確的變形，只能大概地算出一些近似值。但這模型也有很大的意義，不但對該古塔適用，也能把它用

33、于其他古塔的變形分析，從而保證了對古塔進行維護的快速性和高效性。6.2 模型的推廣 由圖1-1中同層視為平面時可看出古塔每層的相鄰連線圖像相似于是正八邊形，但相對于正八邊形確有一定差距，但從古塔的設計理念應該使古塔完全成為正八邊形從而增強整個塔的穩固性，但是實際塔的形狀確實產生了誤差從而使整個塔的形狀產生的變化。因為原本的形狀為8條等邊的正八邊形，而古塔各層形狀的變化必引起8條邊的長度變化，使8條邊的長度變為不等，即可用八邊形邊長的長度的標準差來判斷該層形狀的變化程度，用各邊長長度大小的比較來判斷該層變化大概形狀。七、參考文獻【1】數學建模論文格式要求 360問答 【2】李大潛主編 中國大學生

34、數學建模競賽 高等教育出版社 (1998)附件附錄一：close allclcclearx1=561.4782 563.9092 567.9981 571.3812 572.0142 569.5162 565.4062 562.1112; 561.7463 564.1008 568.0171 571.243 571.8122 569.4066 565.4625 562.3235; 562.0086 564.289 568.0358 571.107 571.6131 569.3 565.5175 562.5296; 562.2183 564.4385 568.0507 570.9988 571.

35、4555 569.2151 565.562 562.6956; 562.4427 564.5997 568.0667 570.8829 571.2857 569.1237 565.6087 562.8727; 562.8807 564.9208 568.1219 570.6945 571.0345 568.984 565.7192 563.2215; 563.0424 565.0389 568.1419 570.6257 570.941 568.93 565.7612 563.3528; 563.3012 565.2289 568.1741 570.5146 570.7907 568.8459

36、 565.8278 563.5608; 563.5593 565.4173 568.2063 570.4043 570.6414 568.7614 565.8943 563.7683; 563.7822 565.5507 568.1927 570.2892 570.5166 568.7563 566.0125 563.9956; 564.0405 565.6943 568.1455 570.1131 570.3292 568.7197 566.135 564.2616; 564.2984 565.8377 568.0986 569.9368 570.1419 568.6831 566.2569

37、 564.5268; 564.5544 565.9834 568.0609 569.7744 569.9704 568.6554 566.3814 564.7894y1=521.4177 518.0897 517.4028 519.8707 523.9807 527.3387 527.9887 525.5227; 521.4074 518.2185 517.5595 519.9587 523.9056 527.1539 527.7745 525.3761; 521.3974 518.345 517.7137 520.0452 523.8316 526.975 527.5655 525.2338

38、; 521.3894 518.4456 517.8365 520.114 523.773 526.831 527.397 525.119; 521.3808 518.5539 517.9679 520.1879 523.7099 526.6769 527.2189 524.9969; 521.363 518.8401 518.318 520.3849 523.6306 526.288 526.7638 524.6852; 521.3565 518.9453 518.4445 520.456 523.6 526.1379 526.5906 524.5678; 521.346 519.1146 5

39、18.6487 520.573 523.5537 525.9031 526.3165 524.382; 521.3356 519.2826 518.8526 520.6886 523.5066 525.6676 526.0426 524.1966; 521.3422 519.3767 518.9714 520.7323 523.4077 525.462 525.8122 524.045; 521.3802 519.5287 519.1543 520.8171 523.3109 525.2229 525.5456 523.8776; 521.4182 519.6805 519.3369 520.

40、9021 523.2143 524.9848 525.2806 523.7107; 521.4476 519.8195 519.5048 520.9766 523.1144 524.7485 525.0184 523.5424z1=1.777 1.802 1.773 1.757 1.765 1.723 1.76 1.759; 7.322 7.335 7.298 7.287 7.292 7.302 7.31 7.326; 12.745 12.769 12.734 12.725 12.739 12.693 12.724 12.729; 17.081 17.088 17.065 17.052 17.

41、049 17.048 17.092 17.083; 21.721 21.741 21.697 21.691 21.695 21.699 21.706 21.725; 26.251 26.285 26.293 26.248 26.182 26.095 26.132 26.202; 29.864 29.894 29.893 29.851 29.776 29.734 29.75 29.835; 33.377 33.421 33.406 33.374 33.299 33.235 33.259 33.348; 36.882 36.919 36.911 36.874 36.799 36.746 36.76

42、6 36.853; 40.195 40.21 40.226 40.217 40.165 40.022 40.109 40.145; 44.467 44.477 44.492 44.477 44.438 44.3 44.389 44.421; 48.737 48.745 48.758 48.746 48.702 48.57 48.653 48.687; 52.859 52.873 52.878 52.873 52.796 52.686 52.773 52.809x=x1,x1(:,1)y=y1,y1(:,1)z=z1,z1(:,1)d=567.336 522.2148 55.091hold on

43、for i=1:13plot3(x(i,:),y(i,:),z(i,:),'k.') plot3(x(i,:),y(i,:),z(i,:),'k-') plot3(mean(x1(i,:),mean(y1(i,:),mean(z1(i,:),'k.') cen(i,:)=mean(x1(i,:),mean(y1(i,:),mean(z1(i,:)endplot3(d(1,1),d(1,2),d(1,3),'k.') X,Y=meshgrid(cen(:,1),cen(:,2) Z=-58296.37472+92.78414147.

44、*X+10.94248841.*Y mesh(X,Y,Z) axis(560 575 515 530 0 60)xlabel('x', 'FontWeight', 'bold');ylabel('y', 'FontWeight', 'bold');zlabel('z', 'FontWeight', 'bold'); zoom onaxis squareview(45,45) 附錄二：close allclcclearx1=561.4782 563.90

45、92 567.9981 571.3812 572.0142 569.5162 565.4062 562.1112; 561.7463 564.1008 568.0171 571.243 571.8122 569.4066 565.4625 562.3235; 562.0086 564.289 568.0358 571.107 571.6131 569.3 565.5175 562.5296; 562.2183 564.4385 568.0507 570.9988 571.4555 569.2151 565.562 562.6956; 562.4427 564.5997 568.0667 570

46、.8829 571.2857 569.1237 565.6087 562.8727; 562.8807 564.9208 568.1219 570.6945 571.0345 568.984 565.7192 563.2215; 563.0424 565.0389 568.1419 570.6257 570.941 568.93 565.7612 563.3528; 563.3012 565.2289 568.1741 570.5146 570.7907 568.8459 565.8278 563.5608; 563.5593 565.4173 568.2063 570.4043 570.6414 568.7614 565.8943