劉之林，李興莉

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

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

1.1 問題提出

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

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

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

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

1.2 模型假設(shè)

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

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

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

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

2 模型的建立與求解

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

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

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

2.2 計算各層中心點的坐標(biāo)

2.2.1 模型建立

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

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

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

2.2.2 模型求解

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

2.2.3 誤差分析

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

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

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

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

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

2.3 計算各層傾斜變形量

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

圖1 古塔各層測量點在地平面的投影

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

設(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ù)解析幾何知識，在ΔO'1OkO''k中有：

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

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

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

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

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

2.4 計算各層的扭曲變形量

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

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

圖2 第K層與第一層測量點的扭曲角

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

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

2.5 計算古塔的彎曲變形量

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

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

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

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

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

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

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

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

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

2.6.1 模型建立

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

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

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

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

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

計算原始數(shù)據(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）

計算均值數(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ù)。利用最小二乘法可以求得：

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

其中

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

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

第三步，檢驗預(yù)測值。

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

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

（2）級比偏差檢驗：令ρ(k)為級比偏差

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

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

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

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

計算級比，得到

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

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

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

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

還原模擬式為：

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

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

表7 傾斜角分析

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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

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

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