2012高教社杯全國大學生數學建模競賽承諾書我們仔細閱讀了中國大學生數學建模競賽的競賽規則.我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網上咨詢等）與隊外的任何人（包括指導教師）研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽規則的,如果引用別人的成果或其他公開的資料（包括網上查到的資料），必須按照規定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽規則，以保證競賽的公正、公平性。如有違反競賽規則的行為，我們將受到嚴肅處理。我們授權全國大學生數學建模競賽組委會，可將我們的論文以任何形式進行公開展示（包括進行網上公示，在書籍、期刊和其他媒體進行正式或非正式發表等）。我們參賽選擇的題號是（從A/B/C/D中選擇一項填寫）：C 我們的參賽報名號為（如果賽區設置報名號的話）：所屬學校（請填寫完整的全名）：臨滄師范高等專科學校參賽隊員(打印并簽名)：1.饒兵2.田興平3.熊生指導教師或指導教師組負責人(打印并簽名)：穆志勇日期：2013年9賽區評閱編號（由賽區組委會評閱前進行編號）：2012高教社杯全國大學生數學建模競賽編號專用頁賽區評閱編號（由賽區組委會評閱前進行編號）：賽區評閱記錄（可供賽區評閱時使用）：評閱人評分備注全國統一編號（由賽區組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：2013高教社杯全國大學生數學建模競賽C題論文摘要：文物是人類的重要歷史文化遺產，它不僅反映過去，還可以借鑒未來，文物保護是一項重要的工作，不僅要有高度的文物保護意識，還應該有強有力的措施和先進的科學技術手段。本題給出了某古塔在1986年7月、1996年8月、2009年3月和2011年3月進行的4次觀測數據，數據包含了古塔每一層以及塔尖上的采樣點的三維坐標值（x,y,z）。要求確定古塔各層的中心坐標，分析古塔的傾斜和變形情況，最后再分析古塔的變形趨勢，以便于加強防護措施，保護古塔。本文根據題目附件給出的三維坐標數據，做出了各次測量數據采樣點的三維散點圖，畫出各條相連邊的連線，求出該空間多邊形的中心點坐標值。算出了各樓層的中心坐標值之后，可以畫出4次測量的各樓層中心連線圖，再連接上塔的頂點，就刻畫出古塔的外形。根據該連線圖，我們可以分析古塔的傾斜、扭曲、彎曲等變形情況。再將4次測量畫出的中心連線圖作比較，就可以知道古塔的變形趨勢。關鍵詞：文物保護、古塔、垂徑定理、中心、線性回歸一、問題重述古塔的變形由于長時間承受自重、氣溫、風力等各種作用，偶然還要受地震、颶風的影響，古塔會產生各種變形，諸如傾斜、彎曲、扭曲等。為保護古塔，文物部門需適時對古塔進行觀測，了解各種變形量，以制定必要的保護措施。某古塔已有上千年歷史，是我國重點保護文物。管理部門委托測繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對該塔進行了4次觀測。請你們根據附件1提供的4次觀測數據，討論以下問題：1.給出確定古塔各層中心位置的通用方法，并列表給出各次測量的古塔各層中心坐標。2.分析該塔傾斜、彎曲、扭曲等變形情況。3.分析該塔的變形趨勢。二、問題分析根據數據的散點圖初步判斷古塔是至少具有八個側面的多邊形棱錐體（或者棱臺體），也有可能是圓錐體（或者圓臺體）。每層樓具有八個采樣點，將八個點投影在水平面（x,y平面）上，采樣點的在水平面的坐標值為（x,y），將相連的兩個采樣點用直線連接，連線就是外接圓的玄，做出八條邊（亦及八條外接圓的玄）的中垂線，中垂線的交點（或者多個交點的平均值）就是圓心的水平坐標（xo,yo）,而（外接圓的圓心）是該樓層的中心在水平面上的投影，也就說樓層中心的水平坐標值也是（xo,yo）。再根據該樓層八個采樣點的垂直坐標（z軸）的平均值，來近似的估算中心點的垂直坐標zo的值。最后就得出了該樓層的中心三維坐標值（xo,yo,zo）。三、模型的假設3.1、假設古塔每一層的水平橫截都是中心對稱圖形。（保證中心點坐標的平均值居中）3.2、古塔每層的重心就是它的幾何中心。四、定義與符號說明直線AB的斜率直線AB的中垂線的斜率AB與BC的中垂線的交點xy平面內直線和y軸的夾角θ空間直線和xy水平面的夾角五、模型建立與求解5.1畫散點圖，對塔的形體做初步判斷在第u次測量中，根據數據附件中給出的每次測量的數據，畫出了各次采樣點的散點圖，并將古塔同一樓層的點用線連起來，例如2011年的測量數據采樣點連線圖如下：根據圖像初步判斷塔體是至少具有八個側面的多邊形棱錐體（或者棱臺體），也有可能是圓錐體（或者圓臺體）。并且假設其每個水平截面都是中心對稱圖形。5.2求出古塔每層樓的中心坐標值并畫出連線圖在第u次測量過程中，第j樓的八個觀測點的坐標值分別為：（xj,yj,zj）,樓層序數j=1,2,...8。提取八個采樣點在水平面（x,y）上的投影坐標值：（xj,yj）,j=1,2,...8，連接相連的兩個點，畫得一個八邊形，如圖所示：做出第一條變AB的中垂線，中垂線與AB交于PAB點，坐標值為（x12,y12）,根據中點坐標公式可得：，。直線AB的斜率為：，那么直線AB的中垂線的斜率為：根據直線的點斜式方程，已知直線上的一個點，又已知直線的斜率，則AB的中垂線的方程可以確定。方程為：，帶入的值得：。同理可算得BC的中點PCB的坐標值和BC邊中垂線的斜率。PCB（x23,y23）中，,。也可以確定了BC的中垂線的方程。AB和BC直線的中垂線交于點，已知兩直線的方程，聯立兩條中垂線的方程，可以解得它們的交點坐標值：。以此類推，可以算得所有邊（八條邊）的中垂線的方程，然后出它們兩兩相交的交點坐標。如果八邊形的外接圓存在，那么所有中垂線的交點都應該重合于圓心，但是由于測量誤差和形變，這些中垂線的交點并不重合，因此，取所有交點坐標的平均值，作為外接圓的圓心值，也就是八邊形的中心坐標值。實現的方法是用組合，先算出第1條中垂線（即AB邊的中垂線）與其它7條邊的中垂線的交點，再算第2條中垂線（即BC邊的中垂線）與第3條到第8條中垂線的交點，依次類推，算出所有中垂線之間的交點坐標。如果是八條邊則產生個交點，如果是L條邊則產生個交點。總的交點的坐標值為：,其中i,j為邊的序號，并且要求另外是等價的，只能出現一次。對于中心點的垂直坐標，用,其中L是邊的條數，也是頂點的個數。根據上述方法，算得每層樓的中心點坐標，然后用直線連接，如圖所示：5.3數據和分析從圖中可以發現，第1、2次（1986年和1996年）測量的數據非常接近，圖像幾乎重合；第3、4次（2009年和2011年）測量的數據非常接近，圖像幾乎重合。而在第2次（1996年）測量值和第3次（2009年）測量值之間，存在著較大差異。第1次1986年測量各樓層中心坐標值列表樓層坐標軸12345678910111213樓頂x566.68566.73566.77566.81566.85566.89566.93566.96567567.04567.11567.17567.21567.25y522.7522.66522.63522.59522.56522.53522.5522.47522.44522.42522.39522.36522.32522.24z1.78747.320212.75517.07821.72126.23529.83733.35136.85540.17244.44148.71252.83455.123第2次1996年測量各樓層中心坐標值列表樓層坐標軸12345678910111213樓頂x566.68566.73566.77566.81566.85566.9566.93566.97567.01567.04567.11567.18567.22567.25y522.7522.66522.62522.59522.56522.52522.5522.47522.44522.42522.39522.35522.31522.24z1.7837.314612.75117.07521.71626.2329.83233.34536.84840.16844.43548.70752.8355.12第3次2009年測量各樓層中心坐標值列表樓層坐標軸12345678910111213樓頂x566.76566.81566.84566.88566.91567.02567.06567.12567.19567.24567.28567.32567.37567.34y522.69522.65522.61522.58522.54522.52522.51522.5522.49522.43522.38522.33522.28522.21z1.76457.30912.73217.0721.70926.21129.82533.3436.84440.16144.43348.752.81855.091第4次2011年測量各樓層中心坐標值列表樓層坐標軸12345678910111213樓頂x566.77566.81566.85566.88566.91567.02567.06567.12567.19567.24567.28567.32567.37567.34y522.69522.65522.61522.58522.54522.52522.51522.5522.49522.43522.38522.33522.28522.21z1.76327.290512.72717.05221.70426.20529.81733.33736.82240.14444.42548.68452.81355.0875.3.1靜態分析由于第1、2次（1986年和1996年）測量的數據非常接近，所以選取第2次（1996年）測量的數據為例進行分析。根據1996年數據畫得古塔各樓層中心連線圖，如圖紅色所示，為了便于觀察，相對前一個圖而言，將圖像旋轉了大約180度。在xy水平面（z=0)、xz側面（y=552.8時的平面）和yz側面（x=566.5時的平面）的投影都接近直線，因此可以用線性擬合的方法，將數據擬合為空間直線。因此對連線做了空間直線擬合，擬合后的圖像如圖所示。另外因為樓頂坐標（上圖中的五角星）不是某一樓層的中心，因此我們不放在擬合直線的參考范圍之內。算出了該擬合直線的空間直線方程為x=Az+C1,y=Bz+C2,據此我們可以算出直線在xy水平面內的投影方程為，根據此方程可以算得古塔中心連線的傾斜方向與x軸和y軸的夾角。與y軸的夾角，=-54.439°；那么與x軸的夾角等于144.439°。為了算得空間直線與xy水平面的夾角，我們在空間直線上任取兩點（x1,y1,z1）,(x2,y2,z2),課算得空間直線和它在xy平面內的夾角，可證明這個夾角就是空間直線和xy平面的夾角。設夾角為：θ，則，就得，任取兩點帶入，可算得89.276°，也就說該直線偏離了垂直（z軸）方向0.72353°。第3、4次（2009年和2011年）測量的數據非常接近，因此取第3次（2009年）測量的數據為例進行分析。根據數據畫得如下的圖像：從圖中可以發現，各個樓層的中心連線圖是扭曲的，在xy水平面（z=0)、xz側面（y=552.8時的平面）和yz側面（x=566.5時的平面）的投影都是彎曲的，說明各樓層的中心分布不是線性的，不能簡單的用直線來代替。但我們還是先用線性擬合直線的算法來算一個平均值，后面再具體分析扭曲的情況。用前面的算法，可以算得xy水平投影與y軸夾角為-60.120°，與x軸夾角為：150.12°根據擬合直線的方程，可算得擬合直線和水平面的夾角為89.180°，也就說該直線偏離了垂直方向（z軸方向）0.81989°。根據中心連線及其投影圖，我們將它分為三段來進行計算：1-5樓、5-9樓、9-13樓。分別用直線代替了這三條線段，如圖所示：根據擬合直線的方程算得：夾角名稱樓層段水平投影和y軸的夾角（度）水平投影和x軸的夾角（度）所在直線和xy水平面的夾角（度）偏離垂直線的角度（與z軸的夾角）1-5層-44.0669134.066989.39860.60135-9層-79.6511169.651188.97191.02809-13層-40.2861130.286189.02820.9717根據圖表數據分析，中間段5-9層的角度與上下兩層的角度相差較大，偏離垂直線的角度也是最大，傾斜嚴重，最有可能在此段發生坍塌。5.3.2動態分析。為了比較1996年和2009年的數據，從總體擬合直線上來比較兩組數據，列表如下：數據名測量序數水平投影和y軸的夾角（度）水平投影和x軸的夾角（度）所在直線和xy水平面的夾角（度）偏離垂直線的角度（與z軸的夾角）2009年（第3次）-60.120150.1289.1800.819891996年（第2次）-54.439144.4389.2760.72353差值-5.681-5.690-0.0960.09636可以發現總體上它們的傾斜的角度的差值并不明顯。將1996年和2009年的數據擬合圖的起點重合后，兩次擬合連線的對圖比如下：將2009年（第3次）分為三段擬合直線的夾角，分別和1996年（第2次）的擬合直線的夾角值進行比較。古塔2009年測量所得角度和1996年測量所得角度分層比較表：角度量分段數水平投影和y軸的夾角之差水平投影和x軸的夾角之差所在直線和xy水平面的夾角之差偏離垂直線（z軸）的角度差1-5層10.372-10.3630.1226-0.12225-9層-25.21225.221-0.30410.30449-13層14.153-14.143-0.24780.2481由表格可以看出，5-9層在1996年至于2009年期間，發生了較大的形變，角度變化最大，致使整個塔身發生了扭曲。六、對模型的評價6.1模型優點簡化了塔的模型，確定了每層塔的中心（重心），以連線代替塔體。可以充分展示塔身的傾斜、變形、彎曲和扭曲情況。6.2模型缺點實際的塔體情況可能很復雜，其切面不一定是中心對稱圖形。只能通過平均值來盡量接近中心值。6.3改進如果能增加測量點的數據個數是最好的，在限制現有數據個數的情況下，可以進一步用數據篩選方法，淘汰偏差較大的數據。七、參考文獻[1]李柏年.MATLAB數據分析方法[M].北京:機械工業出版社,2012:12-19.[2]王能超.數值分析簡明教程[M].北京：高等教育出版社，2003：36-37.[3]于曉秋，李欣，張宏禮.MATLAB數值實驗[M].沈陽：遼寧科學技術出版社，2004：34-40.八、附錄附錄一，主程序（保存為文件C2013.m）：gcsj1=xlsread('D:\gcsj.xls','sheet1','C4:E111');%從電子表格中讀取對應測量序數的數據gcsj2=xlsread('D:\gcsj.xls','sheet1','I4:K111');gcsj3=xlsread('D:\gcsj.xls','sheet1','O4:Q111');gcsj4=xlsread('D:\gcsj.xls','sheet1','U4:W111');symsYoXoZoyoxoXYyo2xo2%畫散點連線圖foru=1:4figure;fork=0:12x=gcsj((k*8+1):(8*k+8),1);%讀取第k+1樓層的坐標數據值y=gcsj((k*8+1):(8*k+8),2);z=gcsj((k*8+1):(8*k+8),3);x=x(~isnan(x));%去除空白數據NaNy=y(~isnan(y));z=z(~isnan(z));r1=size(x);r=r1(1);plot3(x,y,z,'r-*');%畫采樣坐標點連線圖gridonholdonplot3([x(1),x(r)],[y(1),y(r)],[z(1),z(r)],'r-*');%畫第1點和最后1個的連線圖xlabel('x軸');ylabel('y軸');zlabel('z軸');endend%計算古塔每個樓層的中心坐標值figure;foru=1:4switchu%讀取第u次測量數據case1gcsj=gcsj1;case2gcsj=gcsj2;case3gcsj=gcsj3;case4gcsj=gcsj4;endYo=0;Xo=0;fork=0:12xq=gcsj((k*8+1):(8*k+8),1);%讀取第k+1樓層坐標數據yq=gcsj((k*8+1):(8*k+8),2);zq=gcsj((k*8+1):(8*k+8),3);Yo=0;Xo=0;ll=0;x=xq(~isnan(xq));%去除空白數據NaNy=yq(~isnan(yq));z=zq(~isnan(zq));r1=size(x);r=r1(1);forj=1:r-2forn=j+1:r-1K1=(x(j)-x(j+1))/(y(j+1)-y(j));%算中垂線斜率K2=(x(j+2)-x(j+1))/(y(j+1)-y(j+2));equation1=sym('Y-(y(j)+y(j+1))/2-K1*(X-(x(j+1)+x(j))/2)=0');%解中垂線交點方程equation2=sym('Y-(y(j+2)+y(j+1))/2-K2*(X-(x(j+1)+x(j+2))/2)=0');[xo,yo]=solve(equation1,equation2,'X','Y');yo2=eval(yo);xo2=eval(xo);Yo=Yo+yo2;%累加中心點坐標值Xo=Xo+xo2;ll=ll+1;%累加中心點個數endK1=(x(r-1)-x(r))/(y(r)-y(r-1));K2=(x(1)-x(r))/(y(r)-y(1));equation1=sym('Y-(y(r-1)+y(r))/2-K1*(X-(x(r)+x(r-1))/2)=0');equation2=sym('Y-(y(1)+y(r))/2-K2*(X-(x(1)+x(r))/2)=0');[xo,yo]=solve(equation1,equation2,'X','Y');yo2=eval(yo);xo2=eval(xo);Yo=Yo+yo2;Xo=Xo+xo2;ll=ll+1;%累加中心點個數endYok(u,k+1)=Yo/ll;%計算中心點坐標平均值Xok(u,k+1)=Xo/ll;Zok(u,k+1)=sum(z)/r;ends=size(gcsj);%開始計算頂點坐標xd=gcsj(105:s(1),1);yd=gcsj(105:s(1),2);zd=gcsj(105:s(1),3);xd2=xd(~isnan(xd));%去除空白數據NaNyd2=yd(~isnan(yd));zd2=zd(~isnan(zd));p1=size(xd2);p=p1(1);xt(u)=sum(xd2)/p;%算得頂點的坐標值（xt,yt,zt）yt(u)=sum(yd2)/p;zt(u)=sum(zd2)/p;huatu(u,Xok,Yok,Zok,xt,yt,zt);%畫各樓層中心連線圖end%1996年數據分析figure;title('1996年測量數據樓層中心連線分析圖');plot3(Xok(2,:),Yok(2,:),Zok(2,:),'r-o')%畫中心坐標連線圖holdon;gridon;xlabel('x軸');ylabel('y軸');zlabel('z軸');plot3([Xok(2,13),xt(2)],[Yok(2,13),yt(2)],[Zok(2,13),zt(2)],'r-pentagram');%最后和頂點連線plot(Xok(2,:),Yok(2,:),'k-o')%中心坐標連線在xy平面上的投影圖plot([Xok(2,13),xt(2)],[Yok(2,13),yt(2)],'k-pentagram');%和頂點的連線在xy平面上的投影plot3(566.5*ones(1,13),Yok(2,:),Zok(2,:),'g-o')%畫連線在zy平面上的投影plot3([566.5,566.5],[Yok(2,13),yt(2)],[Zok(2,13),zt(2)],'g-pentagram');%最后和頂點連線plot3(Xok(2,:),522.8*ones(1,13),Zok(2,:),'b-o')%畫連線在zx平面上的投影plot3([Xok(2,13),xt(2)],[522.8,522.8],[Zok(2,13),zt(2)],'b-pentagram');%最后和頂點連線%1996年數據分析擬合figure;nihe(Xok(2,:),Yok(2,:),Zok(2,:));%調用擬合函數title('1996年測量數據樓層中心連線分析圖');%2009年數據分析打印figure;plot3(Xok(3,:),Yok(3,:),Zok(3,:),'r-o')%畫中心坐標連線圖holdon;gridon;xlabel('x軸');ylabel('y軸');zlabel('z軸');plot3([Xok(3,13),xt(3)],[Yok(3,13),yt(3)],[Zok(3,13),zt(3)],'r-pentagram');%最后和頂點連線plot(Xok(3,:),Yok(3,:),'k-o')%中心坐標連線在xy平面上的投影圖plot([Xok(3,13),xt(3)],[Yok(3,13),yt(3)],'k-pentagram');%和頂點的連線在xy平面上的投影plot3(566.5*ones(1,13),Yok(3,:),Zok(3,:),'g-o')%畫連線在zy平面上的投影plot3([566.5,566.5],[Yok(3,13),yt(3)],[Zok(3,13),zt(3)],'g-pentagram');%最后和頂點連線plot3(Xok(3,:),522.8*ones(1,13),Zok(2,:),'b-o')%畫連線在zx平面上的投影plot3([Xok(3,13),xt(3)],[522.8,522.8],[Zok(3,13),zt(3)],'b-pentagram');%最后和頂點連線title('2009年測量數據樓層中心連線分析圖');figure;nihe(Xok(3,:),Yok(3,:),Zok(3,:));%總體擬合直線，并算出夾角title('2009年測量數據樓層中心連線擬合圖');legend('中心連線圖','擬合后的直線圖');%2009年分段討論：figure;holdon;gridon;xlabel('x軸');ylabel('y軸');zlabel('z軸');%1-5樓數據擬合nihe(Xok(3,1:5),Yok(3,1:5),Zok(3,1:5));%5-9樓數據擬合nihe(Xok(3,5:9),Yok(3,5:9),Zok(3,5:9));%9-13樓數據擬合nihe(Xok(3,9:13),Yok(3,9:13),Zok(3,9:13));title('2009年測量數據樓層中心連線分段擬合圖');Xokc=Xok(3,:)-Xok(2,:);%算兩次測量坐標的差值Zokc=Zok(3,:)-Zok(2,:);Yokc=Yok(3,:)-Yok(2,:);%1996年與2009年擬合直線比較畫圖figure;XOKK=Xok(2,:)+Xok(3,1)-Xok(2,1);%將坐標起點重合以便于比較YOKK=Yok(2,:)+Yok(3,1)-Yok(2,1);nihe2(XOKK,YOKK,Zok(2,:));%1-5樓數據擬合nihe2(Xok(3,1:5),Y