劉之林，李興莉

（重慶房地產(chǎn)職業(yè)學(xué)院，重慶401331）

1 問(wèn)題的提出與模型假設(shè)

1.1 問(wèn)題提出

以國(guó)內(nèi)某古塔為例，管理部門委托測(cè)繪公司先后于1986年7月、1996年8月、2009年3月和2011年3月對(duì)該塔進(jìn)行了4次觀測(cè)。根據(jù)觀測(cè)數(shù)據(jù)提出了以下問(wèn)題：

（1）給出確定古塔各層中心位置的通用方法，并列表給出各次測(cè)量的古塔各層中心坐標(biāo)。

（2）分析該塔傾斜、彎曲、扭曲等變形情況。

（3）分析該塔變形的趨勢(shì)（有關(guān)數(shù)據(jù)參見2013年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽C題）。

1.2 模型假設(shè)

（1）假設(shè)古塔各層質(zhì)量分布是均勻的。

（2）假設(shè)每年自然環(huán)境因素對(duì)古塔變形的影響基本相同。

（3）假設(shè)題中所給測(cè)量數(shù)據(jù)準(zhǔn)確無(wú)誤。

（4）假設(shè)每次測(cè)量時(shí)各測(cè)量點(diǎn)的位置不變。

2 模型的建立與求解

2.1 缺失數(shù)據(jù)的處理

由于在1986年和1996年的測(cè)量數(shù)據(jù)中缺少第十三層第五號(hào)測(cè)量點(diǎn)的數(shù)據(jù)，為了能夠用一個(gè)通用算法求出各年第一層至第十三層各層的中心點(diǎn)坐標(biāo)，本研究利用SPSS軟件分別對(duì)第一層至第十二層的第五號(hào)測(cè)量點(diǎn)的橫坐標(biāo)、縱坐標(biāo)、豎坐標(biāo)分別進(jìn)行擬合，并預(yù)測(cè)出第十三層第五號(hào)測(cè)量點(diǎn)的測(cè)量數(shù)據(jù)。

首先，利用SPSS軟件所提供的十一種模型對(duì)所給數(shù)據(jù)進(jìn)行初步分析，發(fā)現(xiàn)對(duì)于橫坐標(biāo)和縱坐標(biāo)，二次模型效果最好；對(duì)于豎坐標(biāo)，立方模型效果最好。以1986年數(shù)據(jù)為例：擬合優(yōu)度（調(diào)整R方） 分別為 0.979，0.995，1.000， 顯著性檢驗(yàn)值（Sig．）均為 0。經(jīng)過(guò)計(jì)算，1986 年和 1996 年第十三層第五號(hào)測(cè)量點(diǎn)的坐標(biāo)預(yù)測(cè)值分別為（567.986,519.737,53.012）,(567.992,519.73,53.014)

2.2 計(jì)算各層中心點(diǎn)的坐標(biāo)

2.2.1 模型建立

由于存在變形因素的影響，古塔各層的中心點(diǎn)產(chǎn)生了微小的變化，每層的中心點(diǎn)到該層各測(cè)量點(diǎn)的距離只能近似相等，即中心點(diǎn)到各測(cè)量點(diǎn)的距離兩兩差值的絕對(duì)值之和應(yīng)該接近于零。顯然，此和式的值越小，中心點(diǎn)的位置就越精確。所以，古塔各層樓面上得到所有測(cè)量點(diǎn)的距離兩兩差值的絕對(duì)值之和最小的點(diǎn)即為該層樓面的中心點(diǎn)。

設(shè)第 k 層（k=1，2，…，13）的中心點(diǎn)的坐標(biāo)為（x(k)，y(k),z（k）），第 k 層第 i個(gè)（i=1,2,…，8）測(cè)量點(diǎn)的坐標(biāo)為(x(k，i),y(k,i),z(k,i))，d(k,i)為第 k 層中心點(diǎn)到該層第 i個(gè)測(cè)量點(diǎn)距離；mx(k)、my(k)、mz(k)分別為第k層各測(cè)量點(diǎn)橫坐標(biāo)、縱坐標(biāo)、豎坐標(biāo)的最小值，Mx(k)、My(k)、Mz(k)分別為第 k 層各測(cè)量點(diǎn)橫坐標(biāo)、縱坐標(biāo)、豎坐標(biāo)的最大值。

我們以第k層中心點(diǎn)的三個(gè)坐標(biāo)值x(k)、y(k)、z(k)為決策變量，以第k層中心點(diǎn)到各測(cè)量點(diǎn)的距離兩兩差值之和最小為目標(biāo)函數(shù)，以x(k)、y(k)、z(k)分別界于第k層各測(cè)量點(diǎn)的橫坐標(biāo)、縱坐標(biāo)、豎坐標(biāo)的最小值與最大值之間為約束條件，建立非線性規(guī)劃模型：

2.2.2 模型求解

根據(jù)所給數(shù)據(jù)利用MATLAB編程計(jì)算，可以獲得1986年、1996年、2009年、2011年古塔各層的中心位置的坐標(biāo)。其中1986年各層中心點(diǎn)的坐標(biāo)如表1所示。

2.2.3 誤差分析

由于古塔變形，中心點(diǎn)到各測(cè)量點(diǎn)的距離不盡相等，模型的計(jì)算也有可能存在偏差。本研究以每層中心點(diǎn)到各測(cè)量點(diǎn)的距離兩兩差值的平均值作為計(jì)算該層中心點(diǎn)的系統(tǒng)誤差，再計(jì)算出各層系統(tǒng)誤差的平均值，并以此平均值作為計(jì)算古塔中心點(diǎn)坐標(biāo)的平均誤差。

表1 1986年古塔各層中心坐標(biāo)

設(shè)εk為計(jì)算第k層中心點(diǎn)坐標(biāo)的系統(tǒng)誤差，ε為各層系統(tǒng)誤差的平均值，則

利用Matlab計(jì)算結(jié)果如表2所示：

表2 計(jì)算中心點(diǎn)坐標(biāo)的平均系統(tǒng)誤差(單位：m)

2.3 計(jì)算各層傾斜變形量

古塔各層的測(cè)量點(diǎn)在地平面上的投影為一個(gè)八邊形（如圖1），設(shè)第一層中心點(diǎn)為O1，過(guò)O1作第一層所在平面的垂線

圖1 古塔各層測(cè)量點(diǎn)在地平面的投影

O1P（初始中軸線），第k層的中心點(diǎn)為Ok，Ok在第一層所在平面上的投影點(diǎn)為O'k，Ok在O1P上的投影為O''k，Ok到O''k的距離dk即為第k層中心點(diǎn)的傾斜位移，O1Ok與O1P的夾角θk即為第k層相對(duì)于初始中軸線的傾斜角度。

設(shè) Ok的坐標(biāo)為（x(k),y(k),z(k)），O1的坐標(biāo)為（x(1),y(1),z(1)），于是 O''k的坐標(biāo)為（x(1),y(1),z(k)）， lk為Ok到O1的距離，hk為O''k到O1的距離，則

根據(jù)解析幾何知識(shí)，在ΔO'1OkO''k中有：

利用Matlab計(jì)算，可求得各層中心點(diǎn)相對(duì)于初始中軸線的傾斜角度θk與傾斜位移dk，結(jié)果如表3所示。

表3 各層傾斜角度與傾斜位移

出于對(duì)建筑物安全系數(shù)的考慮，θk選取最大值，在現(xiàn)實(shí)生活中，對(duì)一棟建筑物選取最大的傾斜角度，有助于提醒工作人員及時(shí)維修或采取其他措施以達(dá)到保障建筑物安全的目的。在這里選取 θk（k=1,2,…，13）的最大值作為此古塔的傾斜角。古塔各層最大傾斜角度與傾斜位移如表4。

表4 古塔各層最大傾斜角度與傾斜位移

以上結(jié)果表明，從第二層到十三層，各層的傾斜角度、傾斜位移都有呈逐年增加的趨勢(shì)。最大傾斜角度的年均增長(zhǎng)量為0.01187度，最大傾斜位移的年均增長(zhǎng)量為0.01615m。

2.4 計(jì)算各層的扭曲變形量

古塔的扭曲度是指古塔各層的測(cè)量點(diǎn)相對(duì)于底層各對(duì)應(yīng)測(cè)量點(diǎn)存在一個(gè)水平旋轉(zhuǎn)角度。假設(shè)第k層中心點(diǎn)在第一層所在平面上的投影為O'k（k=1,2,…13），第一層各測(cè)量點(diǎn)記為 Nki（k=1,2,…,13,i=1,2,…,nk）。如果古塔沒(méi)有發(fā)生扭曲變換，那么OK、Nki兩點(diǎn)的連線在第一層所在平面上的投影O'kN'ki應(yīng)該與O1N1i重合，但是由于存在扭曲變形，它們之間是不重合的，其夾角即為第k層第i個(gè)測(cè)量點(diǎn)的扭曲度。

設(shè)αki為第k層第i個(gè)測(cè)量點(diǎn)相對(duì)于第一層第i個(gè)測(cè)量點(diǎn)的扭曲角度（如圖2），根據(jù)平面知識(shí)，可知

圖2 第K層與第一層測(cè)量點(diǎn)的扭曲角

其中 kli、kki分別為 O1N1i與 O'kN'ki的斜率。

我們定義每層8個(gè)測(cè)量點(diǎn)相應(yīng)扭曲度的最大值為該層的扭曲度，第一層至第十三層各層扭曲度的平均值為古塔的扭曲度，利用Matlab計(jì)算，各年的平均扭曲度分別為：3.8023771，3.8158038，3.0679394，3.0683079。

2.5 計(jì)算古塔的彎曲變形量

由于存在扭曲變形，古塔各層中心點(diǎn)不在同一平面。為了便于研究古塔的彎曲程度，我們對(duì)各層中心點(diǎn)進(jìn)行坐標(biāo)變換：先將坐標(biāo)原點(diǎn)平移至第一層中心點(diǎn)處，然后將各層中心點(diǎn)繞軸旋轉(zhuǎn)至坐標(biāo)平面x'o'z內(nèi)（如圖3），變換公式為：

圖3 坐標(biāo)變換圖

其中（x01,y01,z01）為第一層中心點(diǎn)坐標(biāo)。經(jīng)過(guò)平移變換和旋轉(zhuǎn)變換以后的各層中心點(diǎn)坐標(biāo)如表5所示（以1986年和1996年數(shù)據(jù)為例）。

表5 變換以后的各層中心點(diǎn)坐標(biāo)

經(jīng)過(guò)以上坐標(biāo)平移和旋轉(zhuǎn)變換以后，各層中心點(diǎn)均旋轉(zhuǎn)至xoz平面上，為了計(jì)算古塔的彎曲程度，我們利用Matlab對(duì)其進(jìn)行曲線擬合。根據(jù)各層中心點(diǎn)散點(diǎn)圖的特點(diǎn)，可選用二次函數(shù)來(lái)進(jìn)行擬合，擬合結(jié)果如表6。

表6 各層中心點(diǎn)所在曲線的擬合參數(shù)值

根據(jù)解析幾何的知識(shí)，曲線的彎曲度可用曲率來(lái)表示，計(jì)算公式為

根據(jù)擬合得到的函數(shù)，利用Matlab計(jì)算求得1986年、1996年、2009年、2011年的平均曲率分別為：0.00003474995，0.00002862823，0.00037136045，0.00036879992。

2.6 古塔的變形趨勢(shì)分析與預(yù)測(cè)

2.6.1 模型建立

在分析古塔的變形趨勢(shì)時(shí)，由于原始數(shù)據(jù)序列中只有四個(gè)數(shù)據(jù)，且時(shí)間間隔不等，用任何常規(guī)的模型都很難得到較高的擬合優(yōu)度。為了更好地分析、預(yù)測(cè)古塔未來(lái)的變化趨勢(shì)，我們采用非等距灰色GM（1,1）模型進(jìn)行分析和預(yù)測(cè)。

第一步，數(shù)據(jù)檢驗(yàn)與處理。

對(duì)原始數(shù)據(jù)序列進(jìn)行檢驗(yàn)和預(yù)處理以保證所建模型的可行性以及預(yù)測(cè)的準(zhǔn)確性。設(shè)原始數(shù)據(jù)序列為:

X(0)={x(0)(k1),x(0)(k2),…,x(0)(kn)}，時(shí)間間距為

Δki=ki-ki-1(i=2,3,…,n)。

計(jì)算原始數(shù)據(jù)序列的級(jí)比

第二步，建立非等距的GM(1,1)模型 。

將數(shù)據(jù)序列X(0)一次帶權(quán)累加生成數(shù)列X(1)={x(1)(k1),x(1)(k2),…,x(1)(kn)}，其中1,2,3…n（其中 Δk1=1）

計(jì)算均值數(shù)列：

z(1)(ki)=0.5x(1)(ki)+0.5x(1)(ki-1)(i=2,3,…,n)。

然后，由一次累加生成序列X(1)建立灰模型GM（1,1），灰微分方程為：

白化微分方程為：

其中a（發(fā)展系數(shù)）與b（灰作用量）為待辨識(shí)參數(shù)。利用最小二乘法可以求得：

[a,b]=(BTB)-1BTY

其中

于是（2）式的解即時(shí)間響應(yīng)方程為：

再做累減生成，可得到模型的還原值，公式如下：

第三步，檢驗(yàn)預(yù)測(cè)值。

（1）殘差檢驗(yàn)：令 ε(k)為殘差

一般要求 ε(k)＜0.2,如果 ε(k)＜0.1 就比較好。

（2）級(jí)比偏差檢驗(yàn)：令ρ(k)為級(jí)比偏差

一般要求 ρ(k)＜0.2,如果 ρ(k)＜0.1 就比較好。

2.6.2 用灰預(yù)測(cè)對(duì)傾斜角度的分析、預(yù)測(cè)

第一步：對(duì)原始數(shù)據(jù)列。

X(0)=0.79868455,0.79876666,0.82121955,0.82326289

計(jì)算級(jí)比，得到

λ=(0.9998972015,0.9726590910,0.99755179957),經(jīng)驗(yàn)證級(jí)比級(jí)數(shù)據(jù)均落在可容覆蓋

(0.6703200460,1.3956124251)內(nèi)，說(shuō)明原始數(shù)據(jù)序列適合灰模型。

第二步：將原始數(shù)據(jù)序列進(jìn)行加權(quán)累加得到。X(1)=（0.79868455，8.78635117，19.46220533，21.10873112）

通過(guò)Matlab編程求解，a=-0.0016479784，b=0.7928815402，傾斜角的時(shí)間響應(yīng)方程為：

還原模擬式為：

經(jīng)計(jì)算得到預(yù)測(cè)還原值：

第三步：殘差檢驗(yàn)。將原始數(shù)據(jù)序列與預(yù)測(cè)還原值分別代入上式中計(jì)算出殘差與級(jí)比偏差，結(jié)果如下表7。

表7 傾斜角分析

從上表可以看出預(yù)測(cè)的結(jié)果和原始數(shù)據(jù)非常吻合，所以采取灰預(yù)測(cè)是一個(gè)有效合理的預(yù)測(cè)方法。

第四步：預(yù)測(cè)最大傾斜角度的變化趨勢(shì)。古塔變形的年均變化量是很小，因此如果預(yù)測(cè)的時(shí)間太短就缺乏實(shí)際意義。我們以原始數(shù)據(jù)序列的時(shí)間區(qū)間向外推1/3作為預(yù)測(cè)的時(shí)間點(diǎn)，即對(duì)2019年的最大傾斜角度進(jìn)行預(yù)測(cè)。

經(jīng)Matlab計(jì)算，2019年古塔的最大傾斜角度約為0.833081。

2.6.3 對(duì)傾斜位移、扭曲度、彎曲度的分析與預(yù)測(cè)

將1986年、1996年、2009年、2011年各年的最大傾斜位移、扭曲度、彎曲度分別代入上面建立的灰模型中，即可計(jì)算出相應(yīng)的待辨識(shí)參數(shù)的值以及時(shí)間響應(yīng)方程，并可預(yù)測(cè)出2019年的變形量，預(yù)測(cè)結(jié)果如下表8。

表8 傾斜位移、扭曲度、彎曲度的預(yù)測(cè)結(jié)果

3 結(jié)論

本研究建立的計(jì)算中心點(diǎn)坐標(biāo)以及確定三種變形量的模型具有簡(jiǎn)單明了、處理方法巧妙，且速度快、效率高、實(shí)用性強(qiáng)的特點(diǎn)，所得結(jié)果對(duì)于古塔管理的相關(guān)部門制定必要的保護(hù)措施具有一定的指導(dǎo)意義。但在對(duì)古塔的變形趨勢(shì)進(jìn)行預(yù)測(cè)時(shí)，由于原始數(shù)據(jù)序列的數(shù)據(jù)量太少，而且2009年的數(shù)據(jù)存在異常變化，雖然我們選用了比較科學(xué)的預(yù)測(cè)方法——灰模型，但仍不可避免地存在一定的偏差。

[1]陳東佐,康玉慶.中國(guó)古塔的維修與保護(hù)[J].太原大學(xué)學(xué)報(bào),2006(4).

[2]李文軍,付世祿.淺析傾斜建筑物糾偏[J].中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào),2000(1).

[3]張志勇,楊祖櫻.MATLAB 教程：R2011a[M].北京：北京航空航天大學(xué)出版社,2011.

[4]劉圣保,張公讓,李巧巧,湯義強(qiáng).非等間距GM（1,1）模背景值的改進(jìn)及其最優(yōu)化[J].合肥工業(yè)大學(xué)學(xué)報(bào),2010（11）：1749-1752.

[5]同濟(jì)大學(xué)數(shù)學(xué)教研室.高等數(shù)學(xué)：上冊(cè),第四版[M].北京：高等教育出版社,2002.

[6]羅應(yīng)婷,楊鈺娟.SPSS統(tǒng)計(jì)分析從基礎(chǔ)到實(shí)踐：第2版[M].北京：電子工業(yè)出版社,2010.