(江蘇建筑職業(yè)技術(shù)學(xué)院 建筑設(shè)備與市政工程學(xué)院，江蘇 徐州　221116)

1　研究背景

雙曲線(xiàn)回歸函數(shù)(以下簡(jiǎn)稱(chēng)雙曲線(xiàn)函數(shù))常用于擬合工程沉降量隨時(shí)間的變化規(guī)律[1-4]。對(duì)于該函數(shù)的回歸計(jì)算，數(shù)理統(tǒng)計(jì)教科書(shū)、以往科技文獻(xiàn)以及生產(chǎn)實(shí)際中常用的方法是：首先，通過(guò)變量代換轉(zhuǎn)化為線(xiàn)性模型，并利用線(xiàn)性回歸方法推求線(xiàn)性回歸系數(shù)；然后，根據(jù)線(xiàn)性回歸系數(shù)反求雙曲線(xiàn)函數(shù)的回歸系數(shù)。這種方法稱(chēng)為線(xiàn)性化回歸方法。該方法看似合理，其實(shí)不然。對(duì)于雙曲線(xiàn)函數(shù)，本文列舉了線(xiàn)性化回歸方法出現(xiàn)回歸失真實(shí)例，并基于雙曲線(xiàn)函數(shù)的因變量與將該函數(shù)線(xiàn)性化后的因變量二者的殘差平方和之間的關(guān)系式，分析了線(xiàn)性化回歸方法可能導(dǎo)致回歸失真的原因；提出應(yīng)采用高斯-牛頓法等非線(xiàn)性回歸方法進(jìn)行雙曲線(xiàn)函數(shù)回歸計(jì)算，并給出了借助MATLAB軟件進(jìn)行求解的方法。

2　采用線(xiàn)性化回歸方法的失真現(xiàn)象與分析

丹土一級(jí)公路03A標(biāo)伍家?guī)X隧道右線(xiàn)k2+546斷面測(cè)點(diǎn)E累計(jì)沉降值隨時(shí)間變化實(shí)測(cè)值見(jiàn)表1，通過(guò)回歸計(jì)算，用于預(yù)測(cè)隧道在該測(cè)點(diǎn)的沉降趨勢(shì)。

采用雙曲線(xiàn)函數(shù)擬合表1數(shù)據(jù)，其回歸模型為

(1)

式中：A，B為待估參數(shù)。

表1　伍家?guī)X隧道測(cè)點(diǎn)E累計(jì)沉降值隨時(shí)間變化的實(shí)測(cè)數(shù)據(jù)

2.1　線(xiàn)性化回歸方法的失真現(xiàn)象

設(shè)v=1/y,x=1/t,將式(1)線(xiàn)性化：

(2)

然后，將表1的樣本觀(guān)察值(ti,yi)轉(zhuǎn)化為(vi,xi),i=1～12，對(duì)式(2)線(xiàn)性回歸得：A=1.074 94，B=0.123 04，擬合直線(xiàn)見(jiàn)圖1，線(xiàn)性回歸的相關(guān)系數(shù)r=0.975 8，樣本點(diǎn)n=12，查相關(guān)系數(shù)檢驗(yàn)臨界值r0.01=0.707 9，可見(jiàn)線(xiàn)性回歸的相關(guān)系數(shù)r大于相關(guān)系數(shù)檢驗(yàn)臨界值r0.01=0.707 9，線(xiàn)性相關(guān)顯著。

根據(jù)上述線(xiàn)性回歸結(jié)果，求得非線(xiàn)性回歸系數(shù)A=1.074 94，B=0.123 04，進(jìn)而用式(3)計(jì)算曲線(xiàn)回歸的相關(guān)指數(shù)(也稱(chēng)確定性系數(shù))R2：

(3)

圖1　實(shí)例1線(xiàn)性化的擬合直線(xiàn)

時(shí)間t/d累計(jì)沉降值yi/mmvi=1/yiy4i^vi=Ax+B(vi-^vi)2(vi-^vi)2y4i^yi=tiA+Bti(yi-^yi)210.811.23460.43051.19800.00130.00060.83470.000621.370.72993.52280.66050.00480.01701.51400.020732.740.365056.36410.48140.01350.76352.07750.438943.260.3067112.94590.39180.00720.81652.55250.500654.030.2481263.76680.33800.00812.13122.95831.148574.310.2320345.07150.27660.00200.68593.61530.482694.150.2410296.61450.24250.00000.00074.12410.0007124.290.2331338.71090.21260.00040.14214.70330.1708144.170.2398302.37380.19980.00160.48355.00450.6963174.210.2375314.14370.18630.00260.82545.36851.3421204.260.2347329.33540.17680.00341.10625.65651.9503224.280.2336335.56380.17190.00381.27935.81732.3633合計(jì)0.04888.25199.1154

結(jié)果為R2=0.448 5。繪出雙曲線(xiàn)回歸函數(shù)的擬合曲線(xiàn)如圖2所示。

圖2　實(shí)例1采用線(xiàn)性化回歸方法所得的擬合曲線(xiàn)

上述計(jì)算結(jié)果表明，采用線(xiàn)性化回歸方法，盡管變量代換后線(xiàn)性回歸的擬合效果較好，但由其求得的雙曲線(xiàn)回歸方程擬合實(shí)測(cè)數(shù)據(jù)的效果并不好，擬合曲線(xiàn)未能反映實(shí)測(cè)點(diǎn)的分布情況。盡管對(duì)雙曲線(xiàn)線(xiàn)性化后的相關(guān)系數(shù)平方r2=0.952 2,接近1，而曲線(xiàn)回歸的相關(guān)指數(shù)R2卻較小，出現(xiàn)了擬合失真現(xiàn)象。

2.2　失真現(xiàn)象分析

(4)

(5)

3　雙曲線(xiàn)函數(shù)的非線(xiàn)性回歸方法及實(shí)現(xiàn)

為提高雙曲線(xiàn)函數(shù)的擬合精度，本文提出應(yīng)采用高斯-牛頓法、麥夸爾特法等非線(xiàn)性回歸方法。

一般地，對(duì)于一元非線(xiàn)性模型

y=f(x,θ)+ε

式中，f為一般函數(shù)；θ為P維參數(shù)向量，即θ=(θ1,θ2,…,θp)′；ε為隨機(jī)誤差項(xiàng)，且ε服從N(0，σ2)。設(shè)y和x具有n組觀(guān)察值(xi，yi)，i=1～n。高斯-牛頓法求解非線(xiàn)性模型參數(shù)的“最小二乘”估計(jì)的參數(shù)遞推公式[6]，寫(xiě)成矩陣形式為

θ(k+1)=θ(k)+[J′(θ(k))J(θ(k))]-1×

J′(θ(k))[(y-f(θ(k))]

(6)

(7)

式中，k為遞推次數(shù)；向量y=(y1,y2,…,yi，…,yn)′，yi為因變量第i個(gè)觀(guān)察值；向量f(θ(k))=[f1(θ(k)),f2(θ(k)), …,fi(θ(k)) , …,fn(θ(k))]′，fi(θ(k))為由非線(xiàn)性方程及第k次迭代參數(shù)計(jì)算的因變量y的第i個(gè)估計(jì)值，i=1～n。

實(shí)例1擬合雙曲線(xiàn)函數(shù)，采用高斯-牛頓法回歸計(jì)算的步驟如下。

(1)分別對(duì)式(1)中參數(shù)A、B求偏導(dǎo)數(shù)：

(9)

(10)

(2)可利用實(shí)測(cè)值中任兩組關(guān)系值求待估參數(shù)A，B的初值θ(0)。例如采用(4，3.26)、(14，4.17)，得θ(0)=(0.39，0.21)。

(3)利用θ(0)、式(8)～(9)及n組實(shí)測(cè)值(ti，yi)，i=1～n，計(jì)算偏導(dǎo)數(shù)矩陣J(θ(0))及f(θ(0))，進(jìn)而根據(jù)式(6)進(jìn)行迭代計(jì)算，直到θ(k)收斂穩(wěn)定，即|θ(k+1)-θ(k)|小于或等于預(yù)先指定的小正數(shù)(例如δ=0.000 5)，從而得到非線(xiàn)性回歸系數(shù)A，B的估計(jì)值。

對(duì)實(shí)例1采用高斯-牛頓法參數(shù)的遞推結(jié)果見(jiàn)表3，誤差平方和2.156 6遠(yuǎn)小于線(xiàn)性化方法的誤差平方和 9.115 4。進(jìn)一步計(jì)算相關(guān)指數(shù)R2=0.869 5, 遠(yuǎn)大于采用線(xiàn)性化方法所得雙曲線(xiàn)回歸方程求得的相關(guān)指數(shù)的R2=0.448 5，擬合曲線(xiàn)如圖3中實(shí)線(xiàn)所示。可見(jiàn)，擬合效果顯著優(yōu)于線(xiàn)性化回歸方法。

表3　實(shí)例1高斯-牛頓法參數(shù)的遞推結(jié)果

圖3　實(shí)例1由不同方法所得擬合曲線(xiàn)比較

上述計(jì)算也可調(diào)用MATLAB軟件中nlinfit函數(shù)進(jìn)行非線(xiàn)性回歸計(jì)算，具體方法如下[7-8]。

function yhat=volumsq(beta,t)

yhat=t./(beta(1)+ beta(2)t)

(2)在命令窗口輸入

t=[1,2,3,4,5,7,9,12,14,17,20,22]

y=[0.81,1.37,2.74,3.26,4.03,4.31,4.15,4.29,4.17,4.21,4.26,4.28]

beta0=[0.39,0.21]′

[beta]= nlinfit(t′,y′,′volumsq′,beta0);

beta

得結(jié)果beta=[0.537 2,0.193 9]′。

實(shí)例2 文獻(xiàn)[3]中溫嶺東海塘軟基公路工程，由觀(guān)測(cè)的時(shí)間t與實(shí)測(cè)累計(jì)沉降量y的關(guān)系值點(diǎn)繪關(guān)系點(diǎn)見(jiàn)圖4[3]，其規(guī)律可用雙曲線(xiàn)函數(shù)式(1)進(jìn)行擬合。采用線(xiàn)性化回歸方法、高斯-牛頓法的計(jì)算結(jié)果見(jiàn)表4；1/y與1/t線(xiàn)性回歸直線(xiàn)如圖5所示，線(xiàn)性化方法和高斯-牛頓法所得雙曲線(xiàn)分別如圖4虛線(xiàn)和實(shí)線(xiàn)所示。由計(jì)算結(jié)果可見(jiàn)，盡管線(xiàn)性相關(guān)系數(shù)的平方r2=0.988 8(r=0.994 4，線(xiàn)性相關(guān)高度顯著)，但由線(xiàn)性化方法所得雙曲線(xiàn)擬合效果并不理想，而采用高斯-牛頓法進(jìn)行非線(xiàn)性回歸顯著優(yōu)于線(xiàn)性化回歸方法。

表4　實(shí)例2不同回歸方法求得雙曲線(xiàn)回歸的計(jì)算結(jié)果

圖4　實(shí)例2不同方法所得的擬合曲線(xiàn)比較

圖5　實(shí)例2線(xiàn)性化的擬合直線(xiàn)

4　結(jié)　語(yǔ)

從理論上分析了雙曲線(xiàn)函數(shù)因變量與其線(xiàn)性化后的因變量的殘差平方和之間的關(guān)系式，該式表明，線(xiàn)性化回歸方法無(wú)法滿(mǎn)足雙曲線(xiàn)函數(shù)的因變量的殘差平方和為最小，并且結(jié)合雙曲線(xiàn)函數(shù)采用線(xiàn)性化回歸方法出現(xiàn)回歸失真的實(shí)例，分析了線(xiàn)性化回歸方法導(dǎo)致其回歸失真的原因。

雙曲線(xiàn)函數(shù)回歸計(jì)算的實(shí)例表明，采用高斯-牛頓法進(jìn)行非線(xiàn)性回歸計(jì)算，其擬合效果顯著優(yōu)于線(xiàn)性化回歸方法。因此，應(yīng)采用高斯-牛頓法等非線(xiàn)性回歸方法進(jìn)行雙曲線(xiàn)函數(shù)的回歸計(jì)算。本文給出了借助MATLAB軟件中的nlinfit函數(shù)進(jìn)行非線(xiàn)性回歸計(jì)算的方法，易于實(shí)現(xiàn)。

參考文獻(xiàn):

[1]倪錦初,李曉偉,張治軍，等.高填方渠道施工期沉降預(yù)測(cè)分析[J]. 人民長(zhǎng)江,2014,45(6):85-88.

[2]王元柱,梁城.隧道監(jiān)控量測(cè)數(shù)據(jù)的回歸分析[J].土工基礎(chǔ),2013,27(5):70-72.

[3]金錫斐,江維文.雙曲線(xiàn)法在軟土地基沉降預(yù)測(cè)中的適用性研究[J].浙江水利科技,2008(5): 47-48.

[4]李彥杰，黃斌，譚凡，等.堆載預(yù)壓法在在深圳機(jī)場(chǎng)軟基處理中的應(yīng)用[J]，人民長(zhǎng)江，2010，41(S):21-23.

[5]張子賢.可線(xiàn)性化的非線(xiàn)性回歸的有關(guān)問(wèn)題與幾種回歸方法的比較[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí), 2015,45(18): 167-173.

[6]袁志發(fā),周靜芋.多元統(tǒng)計(jì)分析[M].北京:科學(xué)出版社,2002.

[7]周品.MATLAB概率與數(shù)理統(tǒng)計(jì)[M].北京:清華大學(xué)出版社, 2012.

[8]謝中華,李國(guó)棟,劉煥進(jìn),等.MATLAB從零到進(jìn)階[M].北京:北京航空航天大學(xué)出版社, 2012.