確定性周期與隨機激勵聯合作用下非線性系統非平穩響應的統計線性化方法

孔 凡，韓仁杰，張遠進

（1.武漢理工大學土木工程與建筑學院，湖北武漢 430070；2.合肥工業大學土木與水利工程學院，安徽合肥 230009；3.武漢理工大學安全科學與應急管理學院，湖北武漢 430070）

引言

工程動力作用具有強烈的隨機性，通常利用非平穩隨機過程描述［1］。以完全非平穩地震動為例：它不僅在強度上具有明顯的上升和衰減過程，而且在不同時間段上能量的頻率分布也不同［2］。因此，地震動非平穩性體現為幅值非平穩和頻率非平穩。學者們提出了各種模型以全面反映地震動的非平穩特性［3-6］。近年來，小波變換以其時-頻聯合分辨特性在地震動建模方面得到了廣泛發展，詳見文獻［7-10］。

復雜非線性系統在隨機激勵下的隨機動力響應是隨機振動研究的重要課題。具體而言，利用Fokker-Planck-Kolmogorov（FPK）方程［11］、隨機平均法［12］、統計線性化（或非線性化）［13］、矩截斷［14］、信息理論（熵）方法［15］以及最近發展的概率密度演化方法［16］等，人們已解析或數值地得到了很多非線性隨機動力系統的精確解或近似解。然而，高效地計算非平穩激勵下復雜非線性系統的非平穩響應仍然是工程隨機振動理論發展面臨的巨大挑戰［17］，具有重要的實際意義。統計線性化方法是隨機振動理論中具有廣泛適用性的非線性方法之一。它可以用于不同的無記憶非線性［18］，甚至滯回非線性系統［19］、平穩或非平穩隨機激勵［17］、多自由度［20］或連續系統［21］，最近的文獻綜述詳見文獻［22］。

上述研究僅考慮激勵為隨機過程的情況。然而，工程實際中，結構往往會同時受到確定性周期和隨機激勵的聯合作用。例如，旋轉式飛機［23］經常受到色噪聲和諧波激勵聯合作用；俘能系統的非線性響應［24］；近斷層脈沖地震作用下的結構等。為此，人們發展了聯合激勵下非線性系統隨機動力響應的系列方法：Chang［25］使用高斯與非高斯線性化求得了Duffing 系統在諧波與高斯白噪聲激勵聯合作用下的平穩響應；Chen 等［26］給出了分數階Duffing 系統在確定性激勵與平穩白噪聲聯合作用下平穩響應的隨機平均法；Zhu 等［27］利用諧波平衡和矩截斷方法研究了Duffing 系統在諧波與隨機激勵聯合作用下的平穩響應；Zhang 等［28］給出了周期與隨機激勵聯合作用下多自由度齒輪結構平穩響應的統計線性化方法。然而，考察聯合激勵下非平穩隨機作用帶來的響應非平穩性，尚未引起人們足夠的重視。這種隨機響應的非平穩性主要來自于三個方面：突加激勵非平穩、與諧波響應耦合的非平穩、隨機激勵自身的變慢非平穩。

注意到統計線性化方法的廣泛適用性，本文提出了一種用于求解確定性周期與非平穩隨機激勵聯合作用下、單自由度非線性系統非平穩響應的統計線性化方法，可視為作者發展的平穩方法的非平穩推廣［29-30］。該方法將系統響應分解為確定性諧波和零均值隨機分量之和，將原非線性運動方程等效地化為一組耦合的、分別以確定性和隨機動力響應為未知量的非線性微分方程。然后，利用統計線性化方法將非平穩隨機激勵作用下的非線性隨機動力方程化為等效線性方程，得到關于線性隨機響應二階矩的李雅普諾夫微分方程。最后，聯立李雅普諾夫微分方程與諧波激勵作用下的確定性微分方程，通過數值算法（如龍格-庫塔法）對這組耦合非線性常微分方程進行求解。數值算例中，利用蒙特卡洛模擬驗證此方法的適用性和精度。

1 動力學方程

單自由度非線性系統在確定性周期和非平穩隨機激勵聯合作用下的運動方程為：

式中m，c和k分別為質量、阻尼和剛度系數；和分別為響應的位移、速度和加速度；為和的非線性函數；Q(t)為零均值均勻調制非平穩隨機激勵，可寫為：

式中a(t)為調制函數；Qs(t)為零均值平穩隨機過程；F(t)為確定性周期激勵，寫為Fourier 展開形式為：

式中Fl為第l次諧波幅值。

激勵的非平穩性導致響應的非平穩性。這里，響應的非平穩性主要來自于：突加激勵、諧波激勵諧和變化以及隨機激勵幅值慢變非平穩。假設響應x(t)可分解為均值μx(t)和方差為的零均值非平穩隨機過程之和，即：

同時，

式中 D[·]表示方差。將式（4）代入式（1）中得：

為方便計，省略了響應量的時間參數。對式（6）兩邊求數學期望得：

式中 E[·]表示數學期望，用式（6）減式（7）得：

式中

2 統計線性化方法

利用統計線性化方法將式（8）化為線性方程：

式中ce(t)，ke(t)為等效線性參數。當隨機激勵Q(t)為高斯過程時，響應x^ (t)也為高斯過程；即使Q(t)為非高斯過程，當系統阻尼較小時，響應也可近似為高斯過程。因此，ce(t)，ke(t)均為時間t的函數，且有：

以下考慮Qs(t)為白噪聲和色噪聲兩種情況。

2.1 Qs(t)為白噪聲

此時，Q(t)為均勻調制白噪聲，即：

式中w(t)為零均值高斯白噪聲，功率譜密度為S0。將式（10）化為狀態方程的形式：

與式（13）對應的李雅普諾夫方程為：

式中v為響應的協方差矩陣，即：

vij(t)為響應(qi，qj)，i，j=1，2 的方差/協方差，且vij(t)=vji(t)；Θ(t)為激勵Q(t)的協方差矩陣，即：

結合響應的高斯性可知，期望：

可得

所有未知量的初始值均為零，可利用龍格-庫塔法求解式（16）～（20），得響應的確定性分量時程與隨機分量協方差矩陣。

2.2 Qs(t)為色噪聲

此時，Q(t)為均勻調制色噪聲，平穩隨機過程Qs(t)的功率譜密度為S(ω)。將式（2）中的Qs(t)表示為白噪聲經過成型濾波器的形式，即：

式中Φ為濾波器的前置輸出，n為成型濾波器維度；Φ(i)表示Φ的第i階導數；λ0，λ1，…，λn-1和ν0，ν1，…，νn-1為濾波器參數；w(t)為零均值高斯白噪聲，功率譜密度為S0。令：

則式（22）可以寫為：

式中

令：

則式（21）可寫為：

聯立式（8）、式（24）和式（26），得：

與式（27）對應的李雅普諾夫方程為：

式中

且Vqq=E{qqT}，VqΦ=VΦq=E{qΦT}VΦΦ=E{ΦΦT}。

同樣地，將式（28）與式（7）聯立，可利用數值方法求解得到所有未知量需要注意的是，VΦΦ的值僅與色噪聲本身有關，可直接確定。具體的計算步驟將在算例中進一步說明。

3 數值算例

作為算例，不失一般性地考慮具有立方非線性剛度的Duffing 振子，即：

式中ε為表示非線性強度的小量。式（7）中非線性項的數學期望為：

從而，式（7）可寫為：

同樣地，式（8）可化為：

依據式（11），等效線性參數為：

式（13）中，系數矩陣：

如果確定性激勵為單諧波的形式，即：

式中a0為確定性激勵幅值，ω0為確定性激勵頻率。

假定調制函數：

式中A為調制函數峰值，μ標識調制函數下降速度。

3.1 Qs(t)為白噪聲

此時，Q(t)為均勻調制白噪聲，即：

將式（33）～（37）代入式（16）～（20）中，利用數值方法求得y，即所有確定性與隨機響應分量。

3.1.1 典型參數設置的情況

選取系統參數m=1，c=0.4，k=1，ε=0.5；確定性激勵參數a0=1，ω0=1 rad/s；隨機激勵參數本文利用蒙特卡洛模擬（Monte Carlo simulation，MCS）驗證所建議方法的適用性。其中，利用譜表現方法生成白噪聲樣本10000 條，乘以調制函數得均勻調制白噪聲。利用本文所建議方法與MCS 得到了位移均值及其標準差，其對比如圖1（a），（b）所示。

圖1 Duffing 系統在諧波與調制白噪聲聯合作用下的位移Fig.1 Displacement of the considered Duffing system subjected to combined harmonic and modulated white noise

從圖1（a），（b）可以看出，在考慮的參數設置情況下，本文建議方法（Proposed Method，PM）與蒙特卡洛模擬的均值幾乎完全吻合。在隨機響應分量標準差方面，所建議方法能很好地捕捉由于突加激勵（由0 上升的過程）、確定性諧和激勵分量（下降過程中呈諧和變化）和隨機激勵分量調制函數帶來的非平穩性（慢變的上升-下降過程）。

3.1.2 非線性強度的影響

下面分析此方法在不同非線性強度下的適用性。其他參數與3.1.1 節選取相同，非線性強度系數ε取0～1。定義確定性諧和分量平均功率為：

式中T為響應持時。定義隨機響應分量標準差的時間平均為：采用所建議方法與蒙特卡洛模擬計算得到的確定性響應分量平均功率對比如圖2所示，標準差時間平均對比如圖3所示。

圖2 非線性強度對Duffing 系統確定性響應平均功率的影響Fig.2 Averaged power of the deterministic response component of the Duffing system with different degrees of nonlinearity

圖3 非線性強度對Duffing 系統隨機響應標準差時間平均的影響Fig.3 Time-averaged standard deviation of the stochastic response component of the Duffing system with different degrees of nonlinearity

從圖2，3 可見，在其他參數不變的情況下，確定性響應平均功率隨著非線性強度增加而降低；位移標準差的時間平均隨非線性強度增加而降低，速度響應標準差的時間平均隨非線性強度增加而略有減小。所有非線性強度情況下，兩種方法所得結果均吻合較好。

3.1.3 激勵幅值的影響

為研究確定性激勵幅值對此方法適用性的影響，保持其他參數不變，確定性激勵幅值a0取0～2 rad/s。利用本文所建議方法與蒙特卡洛模擬計算得到的確定性響應分量平均功率對比如圖4所示；時變標準差的時間平均對比如圖5所示。

圖4 諧波激勵幅值對Duffing 系統確定性響應平均功率的影響Fig.4 Averaged power of the deterministic response component versus the harmonic excitation amplitude

圖5 諧波激勵幅值對Duffing 系統隨機響應標準差時間平均的影響Fig.5 Influence of the harmonic excitation amplitude on the time-averaged standard deviation of the response

由圖4，5 可見，在其他參數不變的情況下，確定性響應平均功率隨諧波激勵幅值增大而增大；位移標準差的時間平均隨諧波激勵幅值增大而減小，速度標準差的時間平均隨諧波激勵幅值增大而增大。不同激勵幅值下，兩種方法所得的結果均吻合較好。

3.1.4 確定性諧和激勵頻率的影響

為研究確定性激勵頻率對此方法適用性的影響，保持其他參數不變，確定性激勵頻率ω0取0.1～2 rad/s。采用本文所建議方法與蒙特卡洛模擬得到的確定性響應分量的平均功率對比如圖6所示；隨機響應非平穩標準差的時間平均如圖7所示。

圖6 諧波激勵頻率對Duffing 系統確定性響應平均功率的影響Fig.6 Influence of the harmonic excitation frequency on the averaged power of the deterministic response component

圖7 諧波激勵頻率對Duffing 系統隨機響應標準差時間平均的影響Fig.7 Influence of the harmonic excitation frequency on the time-averaged standard deviation of the response

可見，諧波激勵頻率對確定性響應平均功率和隨機響應標準差平均有較大影響，存在使他們達到極值的諧波激勵頻率。其中，使確定性響應平均功率達到最大的諧波激勵頻率約為1.4 rad/s，而使隨機響應標準差時間平均值達到最大的諧波激勵頻率約為1.5 rad/s。此外，不同激勵頻率下本文所建議的方法與蒙特卡洛模擬所得結果均吻合較好。

3.1.5 隨機激勵強度的影響

為研究隨機激勵強度對此方法適用性的影響，保持其他參數不變，功率譜強度S0取0～1。采用本文所建議方法與蒙特卡洛模擬得到的確定性響應分量平均功率對比如圖8所示；隨機響應時變標準差的時間平均對比如圖9所示。

圖8 隨機激勵功率譜強度對Duffing 系統確定性響應功率的影響Fig.8 Influence of the stochastic excitation power spectral density on the averaged power of the deterministic response component

圖9 隨機激勵功率譜強度對Duffing 系統隨機響應標準差時間平均的影響Fig.9 Influence of the stochastic excitation power spectral density on the time-averaged standard deviation of the response

可見，隨機激勵功率譜強度對確定性響應平均功率和隨機響應標準差平均有較大影響。由圖8，9 可知隨著隨機激勵功率譜強度的增大，確定性響應平均功率逐漸減小，隨機響應標準差的時間平均逐漸增大。所有隨機激勵強度下兩種方法所得結果均吻合較好。

綜上，對于Duffing 系統在確定性諧波與調制白噪聲聯合作用下的響應，本文所建議的方法在不同參數設置情況下均有良好的適用性。

3.2 Qs(t)為色噪聲

令Q(t)為均勻調制色噪聲，即：

式中Qs(t)為平穩色噪聲。設功率譜密度為金井清譜［31］的形式：

式中ζg為場地阻尼比，ωg為特征周期。則成型濾波器的維度n=2，式（21）和（22）退化為：

式中

則式（24）中

式（26）中

且式（28）中

將式（35），（49）和式（50）代入式（28）中，并與式（31）聯立，可得一組微分方程，具體步驟可見文獻［18］。

用龍格-庫塔法解此非線性常微分方程組可得響應確定性與隨機分量。

3.2.1 典型參數設置的情況

選取系統參數m=1，c=0.4，k=1，ε=0.5；確定性激勵參數a0=0.5，ω0=1.2 rad/s；隨機激勵參數利用譜表現方法生成色噪聲樣本10000條。

利用本文所建議方法與蒙特卡洛模擬得到了位移的非平穩均值及標準差，其對比如圖10（a），（b）所示。

從圖10（a），（b）可以看出，在考慮的參數設置情況下，本文建議方法與蒙特卡洛模擬的均值較吻合。在隨機響應分量標準差方面，所建議方法能很好地捕捉突加激勵、確定性激勵分量和隨機激勵分量調制函數帶來的非平穩性。

圖10 Duffing系統在諧波與調制色噪聲激勵聯合作用下的位移Fig.10 Displacement of the considered Duffing system subjected to combined harmonic and modulated colored noise

3.2.2 非線性強度的影響

同樣地，保持其他參數不變，非線性強度系數ε取0～1。采用所建議方法與蒙特卡洛模擬計算得到的確定性響應分量的平均功率對比如圖11所示；標準差時間平均對比如圖12所示。

圖12 非線性強度對Duffing 系統隨機響應標準差時間平均的影響Fig.12 Time-averaged standard deviation of the stochastic response component of the Duffing system with different degrees of nonlinearity

從圖11，12可以看出，在其他參數不變的情況下，確定性響應平均功率隨著非線性強度增加先增大后降低；位移標準差的時間平均隨非線性強度增加而降低；速度標準差的時間平均在非線性程度較小時略有升高，之后隨非線性程度增加而減小。在所有非線性強度情況下，兩種方法所得結果均吻合較好。

圖11 非線性強度對Duffing 系統確定性響應平均功率的影響Fig.11 Averaged power of the deterministic response component of the Duffing system with different degrees of nonlinearity

3.2.3 激勵幅值的影響

為研究確定性激勵幅值對此方法適用性的影響，保持其他參數不變，確定性激勵幅值a0取0～2。本文所建議方法與蒙特卡洛模擬計算得到的確定性響應分量的平均功率對比如圖13所示，時變標準差的時間平均值對比如圖14所示。

圖13 諧波激勵幅值對Duffing 系統確定性響應平均功率的影響Fig.13 Averaged power of the deterministic response component versus the harmonic excitation amplitude of the Duffing system

圖14 諧波激勵幅值對Duffing 系統隨機響應標準差時間平均的影響Fig.14 Influence of the harmonic excitation amplitude on the time-averaged standard deviation of the Duffing system

由圖13，14 可見，在其他參數不變的情況下，確定性響應平均功率隨諧波激勵幅值增大而增大；位移標準差的時間平均隨諧波激勵幅值增大而減小，速度標準差的時間平均隨諧波激勵幅值增大先增大后減小，這一點與圖5所示的隨機激勵為白噪聲的情況有所不同。在所有確定性激勵幅值的情況下，兩種方法所得的結果均吻合較好。

3.2.4 確定性諧和激勵頻率的影響

為研究確定性激勵頻率對此方法適用性的影響，保持其他參數不變，確定性激勵頻率ω0取0.1～2 rad/s。本文所建議方法與蒙特卡洛模擬得到的確定性響應分量的平均功率對比如圖15所示；隨機響應非平穩標準差的時間平均如圖16所示。

圖15 諧波激勵頻率對Duffing 系統確定性響應平均功率的影響Fig.15 Influence of the harmonic excitation frequency on the averaged power of the deterministic response component

圖16 諧波激勵頻率對Duffing 系統隨機響應標準差時間平均的影響Fig.16 Influence of the harmonic excitation frequency on the time-averaged standard deviation of the response

可見，諧波激勵頻率對確定性響應平均功率和隨機響應標準差的時間平均有較大影響，存在使它們達到極值的諧波激勵頻率。對于前者，達到極值的頻率為1.2 rad/s；對于后者，達到極值的頻率為1.3 rad/s。此外，不同激勵頻率下本文所建議方法與蒙特卡洛模擬所得結果均吻合較好。

綜上所述，對于Duffing 系統在確定性諧波與隨機激勵聯合作用下的響應，本文所建議方法在不同參數設置的情況下，均有良好的適用性。所建議方法的計算效率較蒙特卡洛模擬有顯著優勢。以調制白噪聲與確定性激勵聯合作用下的響應計算為例，所建議方法僅需0.01 s，而10000 個樣本的蒙特卡洛模擬需要1.62 s，效率提高162 倍。

4 結論與展望

提出了一種用于求解確定性諧波與非平穩隨機激勵聯合作用下，單自由度非線性隨機動力系統非平穩響應的統計線性化方法。首先，將系統響應分解為確定性周期和零均值隨機分量之和，得到了與原非線性運動方程等效的、兩個耦合的且分別以確定性和隨機動力響應為未知量的非線性運動微分方程。隨后，利用非平穩統計線性化方法，將得到的非線性隨機運動方程化為了等效線性方程。最后，同時考慮與等效線性隨機微分方程對應的李雅普諾夫方程和諧波激勵下的確定性運動方程，利用數值方法聯立求解了所有確定性未知量與隨機未知量。蒙特卡洛模擬驗證了此方法的適用性和精度。

數值算例表明，無論隨機噪聲為調制白噪聲或調制色噪聲，該方法均能準確地求解系統的確定性周期響應；求解的隨機動力響應分量的非平穩二階矩也具有較好的精度。此外，該方法適用于不同的參數設置情況。可進一步將該方法拓展應用于（多自由度）滯回非線性系統和分數階非線性系統、（不可分）完全非平穩激勵等情況。