時(shí)域穩(wěn)定的基礎(chǔ)頻響連續(xù)有理近似參數(shù)識(shí)別方法

唐貞云，王志宇，杜修力

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

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

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

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

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

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

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

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

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

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

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

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]，每個(gè)多階系統(tǒng)均由一些列一階和二階子系統(tǒng)組合而成?；谠摳拍?，將式(2)的分母多項(xiàng)式通過部分分式展開分解為l個(gè)二階項(xiàng)與N?2l個(gè)一階項(xiàng)的乘積，而分子多項(xiàng)式不變，如式(4)所示：

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

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

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

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

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

由于共軛根為復(fù)數(shù)，根據(jù)穩(wěn)定性條件：極點(diǎn)實(shí)部小于0。可以得到：

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

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

圖7 實(shí)數(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í)別方法

得到了新的有理近似函數(shù)形式及其穩(wěn)定參數(shù)取值范圍后，如何取得最優(yōu)的擬合結(jié)果是一個(gè)函數(shù)尋優(yōu)問題。優(yōu)化算法可以分為全局優(yōu)化算法和局部?jī)?yōu)化算法，常見的全局優(yōu)化算法有進(jìn)化算法、遺傳算法、粒子群優(yōu)化算法等等，而常見的局部算法有單純形法、序列二次規(guī)劃算法等等。本文以Wang等[35]將遺傳算法與單純形法結(jié)合形成的混合算法為基礎(chǔ)，采用序列二次規(guī)劃算法(SQP算法)代替單純形法，形成遺傳-SQP算法用于求解連續(xù)時(shí)間有理近似函數(shù)的待定系數(shù)。遺傳算法作為全局優(yōu)化算法，不需要給出合理的初始值，可以從多個(gè)任意初始值開始尋優(yōu)，搜索范圍大，利于全局擇優(yōu)。而SQP算法作為局部?jī)?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é)分別采用簡(jiǎn)單與復(fù)雜的基礎(chǔ)頻響函數(shù)進(jìn)行對(duì)比研究。具體參數(shù)辨識(shí)流程如圖9所示，首先根據(jù)頻響函數(shù)曲線獲得函數(shù)的實(shí)部與虛部，然后根據(jù)第二節(jié)穩(wěn)定范圍編寫程序，通過遺傳算法與SQP算法迭代識(shí)別，在滿足精度要求后將結(jié)果輸出。本文進(jìn)行參數(shù)識(shí)別時(shí)采用式(4)給出的新形式有理函數(shù)進(jìn)行辨識(shí)，辨識(shí)的目標(biāo)函數(shù)為：

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

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

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

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

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

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

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

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

圖11 圓形基礎(chǔ)的水平運(yùn)動(dòng)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ǔ)動(dòng)力剛度隨頻率變化劇烈，因此利用連續(xù)時(shí)間有理近似函數(shù)擬合具有挑戰(zhàn)性[32]。本文分別采用算例1中的兩種方法，取階數(shù)N=3/5/7/9/11階進(jìn)行擬合。函數(shù)的擬合精度采用式(12)進(jìn)行評(píng)價(jià)，辨識(shí)結(jié)果如圖12、圖13所示，辨識(shí)效率與誤差和曲線參數(shù)如表3、表4所示。

表3 圓形基礎(chǔ)水平頻響函數(shù)辨識(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ù)辨識(shí)結(jié)果Fig.12 Identification results of proposed method for rational functions of different orders

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

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

4 結(jié)論

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

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

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