唐貞云，王志宇，杜修力

(1. 北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)試驗(yàn)室，北京 100124；2. 中國地震局工程力學(xué)研究所，中國地震局地震工程與工程振動重點(diǎn)實(shí)驗(yàn)室，黑龍江，哈爾濱 150080)

在土結(jié)動力相互作用問題的研究中，基礎(chǔ)部分通常采用頻響函數(shù)進(jìn)行描述，其在連續(xù)時間情況下可以被近似為有理近似函數(shù)，該函數(shù)的穩(wěn)定性與精度決定了基礎(chǔ)時域模型的穩(wěn)定性與精度，進(jìn)而影響土結(jié)相互作用系統(tǒng)動力時程分析。目前已有大量關(guān)于基礎(chǔ)頻響函數(shù)的研究，如：楊冬英等[1 ? 6]對多種不同土體條件中的樁基礎(chǔ)計(jì)算了其動力剛度；Valeti等[7 ? 11]分別針對筏板基礎(chǔ)、嵌入式基礎(chǔ)、條形基礎(chǔ)和樁筏基礎(chǔ)提出了基礎(chǔ)動力剛度數(shù)值計(jì)算方法；陳少林等[12 ? 15]研究了不同土體條件和基礎(chǔ)物理參數(shù)對基礎(chǔ)動力剛度的影響。為便于時域動力分析，在上述頻響函數(shù)基礎(chǔ)上，基于連續(xù)時間有理近似函數(shù)建立了大量時域集中參數(shù)模型，如：Wolf模型[16]、Wu-Lee模型[17]、Du-Zhao模型[18]、Wang-Liu模型[19]等。如何基于連續(xù)時間有理近似函數(shù)建立精確、穩(wěn)定的基礎(chǔ)時域動力模型是實(shí)現(xiàn)土結(jié)動力相互作用時域分析的關(guān)鍵。大量學(xué)者基于有理近似函數(shù)開發(fā)了頻響函數(shù)時域動力模型參數(shù)辨識方法，如：Wolf等[16 ? 17,19 ? 23]分別采用有理近似函數(shù)擬合了基礎(chǔ)頻響函數(shù)，并將其展開為部分分式的形式，建立了不同的一階與二階彈簧—阻尼模型，給出了確保辨識模型穩(wěn)定的后驗(yàn)方法，但難以保證辨識效率；Wu等[24 ? 26]分別采用多項(xiàng)式消除法將有理近似函數(shù)展開為連分式，建立了嵌套的集中參數(shù)模型，但未對模型穩(wěn)定性進(jìn)行評價導(dǎo)致了部分模型在時程計(jì)算時發(fā)散；Okada等[27]采用最小二乘法擬合了基礎(chǔ)頻響函數(shù)，并建立了時域微分方程求解基礎(chǔ)響應(yīng)，但其參數(shù)取值范圍為整個實(shí)數(shù)域且未對模型穩(wěn)定性進(jìn)行判斷或約束；王丕光等[28 ? 29]采用有理近似函數(shù)對水體剛度頻響函數(shù)進(jìn)行了有約束的擬合，保證了函數(shù)的穩(wěn)定性；Han等[30 ? 31]實(shí)現(xiàn)了剛性基礎(chǔ)及地下管道動力剛度矩陣的有理近似，但未給出有效的穩(wěn)定參數(shù)識別方法。以上研究表明，滿足精度要求的頻響函數(shù)的有理近似模型不一定具有穩(wěn)定性，大多數(shù)識別方法均未關(guān)注穩(wěn)定性問題從而可能導(dǎo)致時程計(jì)算發(fā)散。為了保證辨識模型的穩(wěn)定性，杜修力等[32 ? 33]利用罰函數(shù)給失穩(wěn)參數(shù)組合增加懲罰值，保證了識別結(jié)果的穩(wěn)定性，成為最常用的識別方法。

關(guān)于連續(xù)時間有理近似函數(shù)的參數(shù)識別研究，多數(shù)的研究中未考慮模型穩(wěn)定性，少數(shù)學(xué)者采用后驗(yàn)方法保證穩(wěn)定性，但該方法得到有效結(jié)果的效率過低，在研究中難以應(yīng)用。而罰函數(shù)法雖然可以得到穩(wěn)定的結(jié)果，但在擬合過程中并不能保證每一組參數(shù)組合均可穩(wěn)定，因此造成了計(jì)算效率不高，對多自由度系統(tǒng)適用性不強(qiáng)。本文基于線性控制理論將連續(xù)有理近似函數(shù)分解為一系列一階和二階子系統(tǒng)的組合，并從系統(tǒng)控制理論角度給出了有理近似函數(shù)穩(wěn)定的理論約束條件，保證了辨識函數(shù)的絕對穩(wěn)定性，基于此建立了一種新的連續(xù)時間有理函數(shù)參數(shù)識別方法。

1 基礎(chǔ)頻響函數(shù)的連續(xù)時間有理近似

土結(jié)動力相互作用研究中，常采用頻響函數(shù)對基礎(chǔ)進(jìn)行描述：

式中：pj和qj分別為分母多項(xiàng)式和分子多項(xiàng)式的待定系數(shù)，且均為實(shí)數(shù)；a0=ωd/Vs為無量綱頻率，其大小與地基長度d和土體剪切波速Vs相關(guān)；N為有理函數(shù)的階數(shù)。采用式(2)的系數(shù)可建立對應(yīng)的集中參數(shù)模型，包括常見的Wu-Lee模型[17]、Du-Zhao模型[18]，如圖1、圖2所示。雖然采用連續(xù)時間有理近似函數(shù)可以較好地?cái)M合基礎(chǔ)頻響函數(shù)S(ω)，但由于部分滿足精度的有理函數(shù)不具有時域穩(wěn)定性，將其用于土結(jié)動力相互作用系統(tǒng)時域分析時會使動力時程計(jì)算發(fā)散。根據(jù)線性系統(tǒng)的穩(wěn)定性理論[34]，連續(xù)時間有理近似函數(shù)的穩(wěn)定性根據(jù)輸入與輸出的不同，可由其極點(diǎn)或是零點(diǎn)確定，如圖3所示：當(dāng)輸入為位移，輸出為力，即擬合動力剛度函數(shù)時，當(dāng)且僅當(dāng)其極點(diǎn)全部位于S平面的左半平面，也即全部極點(diǎn)的實(shí)部均為負(fù)數(shù)時，函數(shù)穩(wěn)定；而當(dāng)輸入為力，輸出為位移，即擬合動力柔度函數(shù)時，當(dāng)且僅當(dāng)其零點(diǎn)全部位于S平面的左半平面可使函數(shù)穩(wěn)定。假設(shè)式(2)的極點(diǎn)或零點(diǎn)采用sj表示，則其穩(wěn)定條件可表示為：

圖1 Wu-Lee模型Fig.1 Wu-Lee model

圖2 Du-Zhao模型Fig.2 Du-Zhao model

圖3 S平面穩(wěn)定區(qū)域Fig.3 Stable area in S plane

以文獻(xiàn)[17]和文獻(xiàn)[32]中的圓形基礎(chǔ)搖擺頻響函數(shù)為例，其模型如圖4所示，基礎(chǔ)受到力矩(M)并產(chǎn)生轉(zhuǎn)角(θ)，圖中物理量分別表示：基礎(chǔ)的直徑d，土體的質(zhì)量密度ρ，剪切模量G，泊松比ν=0.5。兩文獻(xiàn)的辨識結(jié)果如圖5所示，圖5(a)、圖5(b)分別為復(fù)剛度的實(shí)部與虛部，圖5(c)、圖5(d)分別為以兩者擬合參數(shù)計(jì)算的無量綱位移，由于文獻(xiàn)[17]未對函數(shù)的穩(wěn)定性進(jìn)行判斷，而文獻(xiàn)[32]采用罰函數(shù)對輸出結(jié)果進(jìn)行了約束，因此對于同一頻響函數(shù)，前者進(jìn)行時程計(jì)算是發(fā)散的，而后者是收斂的。

圖4 圓形基礎(chǔ)的搖擺運(yùn)動Fig.4 Rocking motion of circular foundation

圖5 圓形基礎(chǔ)動力響應(yīng)Fig.5 Dynamic response of circular foundation

由上述案例可知，對于連續(xù)時間有理近似函數(shù)，保證其在頻域內(nèi)的擬合精度并不能達(dá)到實(shí)際應(yīng)用的要求，還需要保證有理近似函數(shù)穩(wěn)定性以用于時程分析。

2 穩(wěn)定參數(shù)識別方法

為了確?；谶B續(xù)有理近似的基礎(chǔ)時域動力模型在時程分析中的有效性，需要在時域參數(shù)辨識過程中確保模型對基礎(chǔ)阻抗描述的精確性和模型的時域穩(wěn)定性。前述可知，既有參數(shù)識別方法要么只考慮精度不考慮穩(wěn)定性，要么通過后驗(yàn)法判斷模型穩(wěn)定性，即先滿足頻域響應(yīng)精度，而后進(jìn)行穩(wěn)定性判別，當(dāng)發(fā)現(xiàn)不穩(wěn)定時重新識別參數(shù)，直到獲得同時滿足精度和穩(wěn)定性條件的參數(shù)后停止識別。如果預(yù)先確定滿足穩(wěn)定性的參數(shù)條件，然后在該范圍內(nèi)尋找滿足精度要求的最優(yōu)參數(shù)將會更有效。

2.1 穩(wěn)定參數(shù)約束條件

根據(jù)系統(tǒng)穩(wěn)定性理論[34]，連續(xù)有理近似函數(shù)穩(wěn)定的充分必要條件為全部極點(diǎn)的實(shí)部小于0，因此若是可以通過約束其參數(shù)的方法規(guī)避失穩(wěn)值，即可保證擬合過程的穩(wěn)定。由式(2)可知，需要先得到有理近似函數(shù)的參數(shù)才能獲得其極點(diǎn)，這就是現(xiàn)有方法無法提前確定穩(wěn)定參數(shù)范圍的原因所在。根據(jù)線性控制系統(tǒng)[34]，每個多階系統(tǒng)均由一些列一階和二階子系統(tǒng)組合而成?；谠摳拍?，將式(2)的分母多項(xiàng)式通過部分分式展開分解為l個二階項(xiàng)與N?2l個一階項(xiàng)的乘積，而分子多項(xiàng)式不變，如式(4)所示：

溫度升高會影響雞的休息行為、生理行為，最終導(dǎo)致生產(chǎn)性能下降。但目前關(guān)于偏熱刺激對肉雞休息行為、生理行為和生產(chǎn)性能的影響研究相對較少，因此本次研究對上述幾項(xiàng)指標(biāo)進(jìn)行了研究現(xiàn)狀，具體研究內(nèi)容介紹如下。

式中：s=ia0；xkj為部分分式展開后的分母多項(xiàng)式系數(shù)，k=0,1,2；N為原函數(shù)的階數(shù)；l為二階系統(tǒng)的個數(shù)。由于二階系統(tǒng)可以看作兩個一階系統(tǒng)，因此可以認(rèn)為：當(dāng)N為偶數(shù)時，N=2l，即式(4)中沒有一階項(xiàng)；當(dāng)N為奇數(shù)時，令N?1=2l，即式(4)中僅存在一個一階項(xiàng)。因此，根據(jù)求根公式可以得到式(4)的極點(diǎn)為：

式中：第一個根為二階項(xiàng)的共軛根；第二個根為一階項(xiàng)的根。一階項(xiàng)根的取值范圍可以簡單地得到：

而由于共軛根為實(shí)數(shù)或是復(fù)數(shù)時，所需進(jìn)行穩(wěn)定性判斷的參數(shù)不同，因此將共軛根分為實(shí)數(shù)根與復(fù)數(shù)根兩部分討論：

1)共軛根為復(fù)數(shù)

由于共軛根為復(fù)數(shù)，根據(jù)穩(wěn)定性條件：極點(diǎn)實(shí)部小于0?？梢缘玫剑?/p>

圖6 復(fù)數(shù)根時穩(wěn)定參數(shù)界限Fig.6 Stability parameter boundary for complex roots

2)共軛根為實(shí)數(shù)

圖7 實(shí)數(shù)根時穩(wěn)定參數(shù)界限Fig.7 Stability parameter boundary for real roots

3)穩(wěn)定參數(shù)界限

圖8 穩(wěn)定參數(shù)界限Fig.8 Boundary of stable parameters

2.2 參數(shù)識別方法

得到了新的有理近似函數(shù)形式及其穩(wěn)定參數(shù)取值范圍后，如何取得最優(yōu)的擬合結(jié)果是一個函數(shù)尋優(yōu)問題。優(yōu)化算法可以分為全局優(yōu)化算法和局部優(yōu)化算法，常見的全局優(yōu)化算法有進(jìn)化算法、遺傳算法、粒子群優(yōu)化算法等等，而常見的局部算法有單純形法、序列二次規(guī)劃算法等等。本文以Wang等[35]將遺傳算法與單純形法結(jié)合形成的混合算法為基礎(chǔ)，采用序列二次規(guī)劃算法(SQP算法)代替單純形法，形成遺傳-SQP算法用于求解連續(xù)時間有理近似函數(shù)的待定系數(shù)。遺傳算法作為全局優(yōu)化算法，不需要給出合理的初始值，可以從多個任意初始值開始尋優(yōu)，搜索范圍大，利于全局擇優(yōu)。而SQP算法作為局部優(yōu)化算法可以從遺傳算法得到的全局最優(yōu)解中選擇初始值，避免了因初值選擇不當(dāng)導(dǎo)致的擬合效果差。此外，采用SQP算法代替單純形法可以將無約束優(yōu)化問題變?yōu)橛屑s束優(yōu)化問題，即將式(10)所示參數(shù)約束條件加入尋優(yōu)算法中，計(jì)算流程圖如圖9所示。

3 數(shù)值算例驗(yàn)證

為了驗(yàn)證本文方法的有效性，本節(jié)分別采用簡單與復(fù)雜的基礎(chǔ)頻響函數(shù)進(jìn)行對比研究。具體參數(shù)辨識流程如圖9所示，首先根據(jù)頻響函數(shù)曲線獲得函數(shù)的實(shí)部與虛部，然后根據(jù)第二節(jié)穩(wěn)定范圍編寫程序，通過遺傳算法與SQP算法迭代識別，在滿足精度要求后將結(jié)果輸出。本文進(jìn)行參數(shù)識別時采用式(4)給出的新形式有理函數(shù)進(jìn)行辨識，辨識的目標(biāo)函數(shù)為：

圖9 參數(shù)識別流程Fig.9 Parameter identification process

式中：xj為待識別參數(shù)；ω為外荷載頻率；R(ω)為連續(xù)時間有理近似函數(shù)；S(ω)為動力剛度頻響函數(shù)。

而為了對辨識結(jié)果的精度進(jìn)行定量評估，取精確解與識別函數(shù)之差絕對值積分面積占精確解絕對值積分面積的比值作為誤差評價指標(biāo)，其中實(shí)部與虛部各占50%，如下式所示：

式中：E表示識別模型誤差；R(ω)Re、S(ω)Re、R(ω)Im和S(ω)Im分別表示有理函數(shù)與動力剛度精確值的實(shí)部與虛部。

算例1. 圓形基礎(chǔ)搖擺頻響函數(shù)

對于簡單的頻響函數(shù)，本文以圖4、圖5所示的圓形基礎(chǔ)搖擺運(yùn)動及無約束擬合失穩(wěn)案例為例，其基礎(chǔ)頻響函數(shù)具體數(shù)值取自文獻(xiàn)[32]。參數(shù)辨識分別采用罰函數(shù)法[32]和本文方法進(jìn)行對比，有理函數(shù)階數(shù)取為N=3/4/5/6/7階，辨識結(jié)果如表1、表2和圖10所示，其中，表1計(jì)算效率單位為秒。由圖10(a)、圖10(b)可以看到，兩種方法辨識結(jié)果均可較精確地描述動力剛度，采用式(12)計(jì)算所得誤差均小于1%。事實(shí)上，當(dāng)采用N=3階有理函數(shù)時，辨識誤差約為6.5%，但當(dāng)階數(shù)增加，擬合精度迅速增加，在N=4階及以上時，誤差均低于1%，證明兩種方法擁有很高的精度。由于靜剛度不影響穩(wěn)定性，因此采用表2參數(shù)，將EL-Centro波作為地震動輸入，利用Du-Zhao模型[18]和Newmark-β法進(jìn)行無量綱時程計(jì)算。由圖10(c)、圖10(d)可以看到，兩種方法均可保證時程計(jì)算的穩(wěn)定性，且時程曲線完全一致。由表1可知，在計(jì)算效率方面上，本文方法略優(yōu)于罰函數(shù)法。兩種方法分別采用3階～7階函數(shù)對頻響函數(shù)進(jìn)行識別，擬合用時均隨階數(shù)增加而增加。但罰函數(shù)法擬合用時隨階數(shù)增加是成倍增長的，而本文方法效率損失隨階數(shù)增加較低，用時從罰函數(shù)法的約100%逐漸降低到約65%。原因在于本文方法取值區(qū)間內(nèi)可保證函數(shù)穩(wěn)定，而罰函數(shù)法需要求解高階多項(xiàng)式判斷穩(wěn)定性，且階數(shù)越高效率損失越大。

圖10 圓形基礎(chǔ)搖擺頻響 5階有理函數(shù)辨識結(jié)果Fig.10 Identification results of a 5th-order rational function for rocking frequency response of a circular foundation

表1 圓形基礎(chǔ)搖擺頻響函數(shù)辨識效率Table 1 Identification efficiency of rocking frequency response of a circular foundation

表2 5階有理函數(shù)參數(shù)Table 2 Parameters of 5th-order rational function

算例2. 圓形基礎(chǔ)水平頻響函數(shù)

為驗(yàn)證函數(shù)階次對于辨識精度的影響，分別采用不同階次的有理函數(shù)對圖11所示的圓形基礎(chǔ)水平運(yùn)動復(fù)雜頻響函數(shù)進(jìn)行參數(shù)識別，其值由文獻(xiàn)[36]的圖7得到。圖中PeiωΔt和ueiωΔt分別為基礎(chǔ)受到的

圖11 圓形基礎(chǔ)的水平運(yùn)動Fig.11 Horizontal motion of circular foundation

水平外荷載與產(chǎn)生的位移?；A(chǔ)參數(shù)為：基礎(chǔ)半徑r=10 m，剪切波速Vs=100 m/s，土體泊松比ν=1/3，土體阻尼比ξ=0.05，土層厚度H=2r。由于其基礎(chǔ)動力剛度隨頻率變化劇烈，因此利用連續(xù)時間有理近似函數(shù)擬合具有挑戰(zhàn)性[32]。本文分別采用算例1中的兩種方法，取階數(shù)N=3/5/7/9/11階進(jìn)行擬合。函數(shù)的擬合精度采用式(12)進(jìn)行評價，辨識結(jié)果如圖12、圖13所示，辨識效率與誤差和曲線參數(shù)如表3、表4所示。

表3 圓形基礎(chǔ)水平頻響函數(shù)辨識精度及效率Table 3 Identification accuracy and efficiency of horizontal frequency response for a circular foundation

表4 不同階次有理函數(shù)參數(shù)Table 4 Parameters of rational functions of different orders

圖12 本文方法不同階次有理函數(shù)辨識結(jié)果Fig.12 Identification results of proposed method for rational functions of different orders

圖13 罰函數(shù)法不同階次有理函數(shù)辨識結(jié)果Fig.13 Identification results of penalty function method for rational functions of different orders

由表3、圖12和圖13可以看到，采用不同階次的有理函數(shù)辨識時，兩種辨識方法精度基本相同，辨識精度隨階數(shù)增加逐漸提高，當(dāng)采用3階函數(shù)擬合時，效果較差，整體誤差均約有18%，難以擬合復(fù)雜函數(shù)。但階數(shù)增加到7階時，整體誤差降低到約7%，此時采用有理函數(shù)已可以較好的擬合頻響函數(shù)，而再增加階數(shù)，擬合精度進(jìn)一步提高，但此時提升較小，達(dá)到11階時，整體誤差約為5%。但從辨識效率來看，當(dāng)階數(shù)較低時，兩種方法擬合用時均較短，而隨階數(shù)增加，罰函數(shù)法擬合用時迅速增加。當(dāng)階數(shù)達(dá)到11階時，本文擬合用時僅為98.8 s，而罰函數(shù)法則用時390.44 s，遠(yuǎn)大于本文方法。證明了本文方法在同時保證擬合精度與效率上的優(yōu)勢。

4 結(jié)論

針對基礎(chǔ)頻響連續(xù)時間有理近似函數(shù)的參數(shù)識別難以同時保證穩(wěn)定性、精度及計(jì)算效率的問題，本文提出了時域絕對穩(wěn)定的連續(xù)有理近似函數(shù)參數(shù)識別方法。采用不同基礎(chǔ)頻響函數(shù)，對本文方法和罰函數(shù)擬合方法就穩(wěn)定性、精度及計(jì)算效率等方面進(jìn)行了對比分析，驗(yàn)證了本文方法的有效性與優(yōu)勢。得到主要結(jié)論如下：

(1) 從線性系統(tǒng)理論角度將基礎(chǔ)頻響函數(shù)連續(xù)有理近似分解為一階與二階系統(tǒng)的組合，并根據(jù)其根的穩(wěn)定性條件建立了被辨識參數(shù)的穩(wěn)定界限。據(jù)此，采用遺傳-序列二次規(guī)劃算法建立了時域穩(wěn)定的參數(shù)識別方法，保證了識別函數(shù)時域模型的絕對穩(wěn)定性。

(2) 單自由度基礎(chǔ)頻響函數(shù)對比仿真表明：對于簡單函數(shù)模型本文方法與既有方法擁有同等精度，但在計(jì)算效率上明顯優(yōu)于既有方法，本文方法用時在不同階次均少于既有方法，最大可節(jié)約耗時35%左右；對于復(fù)雜函數(shù)，通過增加有理近似函數(shù)階次可保證模型精度，采用本文方法時也可保證一定的效率，在達(dá)到11階時用時小于100 s，而既有方法需要390 s。這不但同時確保了頻響函數(shù)識別的穩(wěn)定性、精度和效率，更提高了多自由度頻響函數(shù)識別的適用性。