隨機海浪的濾波器模型及仿真

喬倩， 范菊

(上海交通大學 海洋工程國家重點實驗室， 上海 200240)

安全性在船舶設計工作中至關(guān)重要。由大幅度橫搖引起的船舶傾覆問題，往往會造成不可估量的損失。目前船舶穩(wěn)性還不能完全考慮隨機海浪對船舶大幅橫搖以及傾覆的影響。因此，隨機波浪下的船舶橫搖傾覆問題的研究對指導船舶穩(wěn)性規(guī)范修訂方面具有重要意義[1]。

船舶及海洋結(jié)構(gòu)物在隨機波浪上的運動相應于一個動力系統(tǒng)在隨機干擾下的響應，當研究隨機波浪下船舶大幅非線性橫搖運動所導致的傾覆問題時，不可避免地需要用到概率統(tǒng)計的方法進行求解。在時域內(nèi)研究概率特性最有效的工具是馬爾可夫過程理論[2]，當外干擾為白噪聲或過濾白噪聲時，動態(tài)系統(tǒng)的響應是一個馬爾可夫過程，其轉(zhuǎn)移概率密度滿足一個統(tǒng)稱為Fokker-Planck方程的線性微分方程[3]。對于馬爾可夫過程理論應用于船舶及海洋結(jié)構(gòu)物在隨機波浪上的響應方面，Naess等[4]采用路徑積分法通過解Fokker-Planck方程研究了隨機海浪中結(jié)構(gòu)物慢漂運動的響應。朱新穎等[5]對船舶在白噪聲波浪上的傾覆問題的研究中，利用路徑積分法和有限差分法2種方法在時間域內(nèi)求解相應的Fokker-Planck方程，給出船舶橫搖運動響應隨時間演變的概率密度分布。柴威[2]在時域內(nèi)求解Fokker-Planck方程的基礎(chǔ)上，進一步探討了外部激勵、橫搖阻尼以及非線性復原力矩對船舶橫搖運動響應轉(zhuǎn)移概率密度的影響。沈棟等[6]采用路徑積分法通過解Fokker-Planck方程研究了隨機橫浪下船舶橫搖傾覆問題。何誠亮等[7]對船舶在真實海況下的傾覆問題的研究中，利用路徑積分法和有限差分法在時間域內(nèi)求解相應的Fokker-Planck方程，給出接近真實海況下的船舶大幅橫搖運動響應概率密度分布情況。

在應用馬爾可夫理論到船舶在隨機波浪下的運動問題時，由于實際海浪具有有限的譜寬和給定的譜密度函數(shù)，屬于有色噪聲激勵[8]，因此需要將白噪聲外干擾項通過濾波器轉(zhuǎn)換為有色噪聲激勵。通過引入合適的濾波器模型，可以使該問題的解決更接近實際情況。

Spanos[9]將ARMA濾波算法引入到波浪模擬中，對實際的海浪譜進行擬合。為了更好地將ARMA過程應用到運動微分方程，Brockwell等[10]提出了時間連續(xù)的CARMA過程。隨后，濾波方法被廣泛應用于船舶與海洋工程領(lǐng)域的波浪載荷模擬和非線性系統(tǒng)的響應評估中。沈棟等[6]采用由二階時間連續(xù)的AR模型構(gòu)造的濾波器對白噪聲進行濾波處理，得到濾波后的船舶橫搖的隨機微分方程。何誠亮等[7]利用二階CARMA模型制作濾波器，將輸入的白噪聲轉(zhuǎn)化為有色噪聲以模擬真實海況的海浪譜，改善了二階AR濾波器對海浪譜的擬合程度。Chai等[11]采用高階CARMA模型的濾波器，將白噪聲外載荷轉(zhuǎn)換成有色噪聲載荷。

本文利用CARMA過程，給出構(gòu)造不同階數(shù)的濾波器的統(tǒng)一方法，系統(tǒng)的探討了不同階數(shù)的濾波器模型可行性。對各濾波器參數(shù)進行最小二乘擬合，使濾波器生成譜逼近海浪譜，為模擬隨機海浪的有色噪聲激勵奠定基礎(chǔ)。通過將不同濾波器生成的波浪時歷進行統(tǒng)計值比較，從而為海浪模擬提供有效的濾波器形式。

1 隨機海浪的時間序列

像隨機海浪波面升高的時間序列可用自回歸滑動平均過程ARMA建模。ARMA過程的譜密度由多項式的商組成，能夠很好地逼近給定的海浪譜[9]。此外，因為船舶運動通常是用微分方程來模擬的，所以還需要引入時間連續(xù)的CARMA過程[10]。

ARMA的基本思想是將預測對象看作一串隨時間變化且又相互關(guān)聯(lián)的隨機序列[12]，并采用相應的數(shù)學模型進行描述分析，從而在最小方差的層面達到最佳預測。ARMA模型具有3種表現(xiàn)形式：自回歸模型(auto-regressive，AR)，滑動平均模型(moving-average，MA)以及自回歸滑動平均模型(auto-regressive moving-average，ARMA)。

自回歸滑動平均模型一般為時間離散的差分方程，考慮到在時間連續(xù)的情況下，這一差分方程可以改寫為微分方程的形式[13]。引入時間連續(xù)的自回歸滑動平均模型即CARMA模型(continuous-time ARMA)，該模型對應的線性微分方程表達式為：

x(p)(t)+a1x(p-1)(t)+a2x(p-2)(t)+…+apx(t)=

ε(t)+b1ε(1)(t)+b2ε(2)(t)+…+bqε(q)(t)

(1)

式中上標代表對t的導數(shù)的階數(shù)。

考慮一個CARMA(p,q)過程y(t),用狀態(tài)空間的形式來表達：

y=c′u(t)

(2)

式中狀態(tài)向量u(t)∈Rp，并滿足線性Ito微分方程：

du(t)=Au(t)dt+bdW(t)

(3)

式中：W是維納過程；dW(t)=W(t+dt)-W(t)是單位維納過程的一個增量；dW(t)/dt為白噪聲序列，取其功率譜密度S(ω)≡1。

式中0≤q

上述CARMA(p,q)過程對應的譜密度為：

(4)

式中s=iω。值得注意的是該定義所給出的CARMA模型的狀態(tài)空間表示方法并不是唯一的。

2 濾波器模型

工程上一般將影響船舶運動的波浪當作平穩(wěn)隨機過程來處理。輸入白噪聲作為外載荷干擾，并引入合適的濾波器對其進行濾波處理，獲得符合一定統(tǒng)計特性的有色噪聲，以模擬真實海況的海浪譜。

2.1 ITTC雙參數(shù)譜

本文采用ITTC雙參數(shù)譜[14]作為標準海浪譜。它是第15屆國際船模拖曳水池會議所推薦的，以有義波高和特征周期為參數(shù)，結(jié)構(gòu)簡單，適用于充分成長的海浪。其表達式為：

(5)

式中：H1/3是有義波高；T1是特征周期。

圖1是四級海況對應的ITTC雙參數(shù)譜，其中有義波高H1/3=2 m，特征周期T1=8 s。由圖1可以看出，實際的海浪具有一定的譜密度函數(shù)并且譜寬有限，表現(xiàn)為有色噪聲而非白噪聲，因此在隨機問題的求解中需要引入濾波器模型對白噪聲輸入做濾波處理。

圖1 ITTC雙參數(shù)譜曲線形狀

本文以ITTC雙參數(shù)譜為標準海浪譜，構(gòu)造濾波器模型對海浪譜進行逼近，給出了由幾種不同階數(shù)的CARMA(p,q)過程構(gòu)造的濾波器模型。

2.2 二階濾波器

2.2.1 CARMA(2,0)模型

通常情況下，一階自回歸AR(1)模型在隨機運動及載荷領(lǐng)域的應用非常廣泛，而船舶及海洋結(jié)構(gòu)物所涉及的運動響應至少是二階的。因此，首先采用時間連續(xù)的CARMA(2,0)模型來構(gòu)造濾波器。為了將問題簡化處理，此二階模型只考慮了自回歸部分，沒有考慮滑動平均部分，對應的微分方程為：

(6)

式中：輸出用y表示，代表波面升高；輸入為譜密度恒等于1的白噪聲序列，用dW(t)/dt表示；ai是自回歸參數(shù)，bj是滑動平均參數(shù)。該CARMA(2,0)模型濾波器的傳遞函數(shù)為：

(7)

對應的譜密度為：

(8)

為了近似擬合本研究所采用的四級海況對應的ITTC雙參數(shù)譜，采用最小二乘法來確定各參數(shù)值。該濾波器的參數(shù)取值為a1=-0.4，a2=0.447 4，b0=0.196 4。將該濾波器所對應的功率譜曲線與ITTC雙參數(shù)譜進行對比，如圖2所示。

圖2 CARMA(2,0)濾波器生成譜與ITTC譜對比

該濾波器生成譜可以大致體現(xiàn)隨機海浪的統(tǒng)計特性，然而在低頻部分與實際情況相差較大，因此需要采取一定方法對上述缺陷加以改善。

2.2.2 CARMA(2,1)模型

AR模型適合處理全極點譜，MA模型適于處理全零點譜。而現(xiàn)有的用于描述隨機海浪的波譜同時包含極點和零點，需要用極零點模型表達。因此在前述的二階濾波器基礎(chǔ)之上，增加對滑動平均部分的考慮，采用CARMA(2,1)模型構(gòu)造濾波器，對應微分方程的表達式為：

(9)

該CARMA(2,1)模型濾波器的傳遞函數(shù)為：

(10)

對應的頻譜公式為：

(11)

采用最小二乘法對本研究所采用的ITTC雙參數(shù)譜曲線近似擬合，該CARMA(2,1)濾波器模型的參數(shù)取值為a1=-0.37，a2=0.367 4，b1=0.284 6。

圖3為該濾波器對應的有色噪聲功率譜與ITTC雙參數(shù)譜。同時對比圖2和圖3可以發(fā)現(xiàn)，相較于第1種只考慮自回歸模型的CARMA(2,0)濾波器，第2種同時考慮自回歸與滑動平均的CARMA(2,1)濾波器模型對應的有色噪聲功率譜更接近ITTC雙參數(shù)譜，能夠更好地反映隨機海浪的統(tǒng)計特性。

圖3 CARMA(2,1)濾波器生成譜與ITTC譜對比

從圖3可以看出，雖然滑動平均部分的加入改善了自回歸模型的不足，但是CARMA(2,1)模型對實際海浪譜的擬合仍存在缺陷，在低頻段及高頻段都略高于ITTC雙參數(shù)譜。為了更加精確地模擬描述真實海浪統(tǒng)計特性的海浪譜，需要考慮構(gòu)造更高階數(shù)的濾波器模型。

由于三階濾波器模型對ITTC雙參數(shù)譜擬合程度的改善效果不太明顯，下面給出四階濾波器模型。

2.3 四階濾波器

為了使濾波后的有色噪聲更加符合真實海浪的統(tǒng)計特性，需要用到更高階的濾波器模型。四階濾波器依然是基于自回歸滑動平均模型構(gòu)造得到，這里探討不同四階濾波器模型的擬合效果。

2.3.1 CARMA(4,1)模型

由四階自回歸模型和一階滑動平均模型構(gòu)成的CARMA(4,1)濾波器表達式為：

(12)

對應的頻率響應函數(shù)為：

(13)

對應的功率譜密度函數(shù)為：

(14)

該CARMA(4,1)四階濾波器的各參數(shù)取值為a1=-0.312，a2=1.275，a3=0.073 8，a4=0.372，b1=0.153。濾波器生成譜見圖4。

圖4 CARMA(4,1)濾波器生成譜與ITTC譜對比

從圖4可以看出，四階濾波器的生成譜對ITTC雙參數(shù)譜的擬合度明顯優(yōu)于二階濾波器。但是在細節(jié)方面，低頻區(qū)域略高于ITTC雙參數(shù)譜，高頻區(qū)域略低于ITTC譜。而根據(jù)前述二階濾波器的模擬情況，高階的MA模型恰恰可以解決這個問題。高階的MA模型可以降低低頻區(qū)域的譜密度值，并且提高高頻區(qū)域的譜密度值。下文將采用高階MA模型的四階CARMA濾波器。

2.3.2 CARMA(4,2)模型

CARMA(4,2)濾波器由四階自回歸模型和二階滑動平均模型構(gòu)成，表達式如下：

(15)

傳遞函數(shù)的表達式為：

(16)

對應的功率譜密度函數(shù)為：

(17)

該四階濾波器的各參數(shù)取值為a1=-0.932，a2=1.265，a3=-0.41，a4=0.243，b2=0.201。濾波器生成譜見圖5。

圖5 CARMA(4,2)濾波器生成譜與ITTC譜對比

從圖5可以看出，由四階AR模型和二階MA模型組合而成的CARMA(4,2)濾波器對應的有色噪聲功率譜，除了在低頻區(qū)域僅有一小部分略高于ITTC雙參數(shù)譜，在大部分區(qū)域近乎完全擬合。由此可見，四階濾波器模型已能夠很好地實現(xiàn)對實際海浪譜的模擬。不過，為了改善濾波器對應的功率譜密度在低頻區(qū)域略高于實際海浪譜的這個缺陷，考慮引入三階的MA模型做進一步探討。

2.3.3 CARMA(4,3)模型

由四階自回歸模型和三階滑動平均模型構(gòu)成的CARMA(4,3)濾波器表達式為：

(18)

對應的傳遞函數(shù)的表達式為：

(19)

對應的功率譜密度函數(shù)為：

(20)

CARMA(4,3)濾波器模型的各參數(shù)取值為a1=-0.677 5，a2=0.9，a3=-0.237，a4=0.143，b3=0.183 5。濾波器生成譜見圖6。

圖6 CARMA(4,3)濾波器生成譜與ITTC譜對比

對比圖6與圖5，MA模型階數(shù)的增加使得CARMA(4,3)濾波器生成譜明顯改善了在低頻區(qū)域?qū)γ枋鰧嶋H海浪統(tǒng)計特性的ITTC雙參數(shù)譜的模擬效果。在高頻區(qū)域，該濾波器生成譜不能完全趨近于零，但基本接近零。整體來看，CARMA(4,3)數(shù)學模型濾波器的生成譜與ITTC雙參數(shù)譜高度契合，可以用來實現(xiàn)對真實海浪統(tǒng)計特性的描述。

從上面的分析可以看出，四階濾波器CARMA(4,3)對ITTC雙參數(shù)譜的擬合效果最好。但是濾波器階數(shù)的提高會大大增加實際問題求解過程中的計算成本[15]。例如在利用有限差分法數(shù)值求解船舶橫搖運動響應轉(zhuǎn)移概率密度分布隨時間的變化時，在每一個遞進的時間步長內(nèi)，若只考慮用白噪聲作為船舶橫搖運動微分方程的外干擾輸入項，則對應的二維FPK方程的求解區(qū)域按50×50均勻網(wǎng)格離散化，網(wǎng)格數(shù)為502；若引入二階濾波器模型對白噪聲進行濾波處理，則對應的四維FPK方程的求解區(qū)域按50×50×50×50均勻網(wǎng)格離散化，網(wǎng)格數(shù)為504；若將四階濾波器與船舶橫搖運動微分方程聯(lián)立，對白噪聲進行濾波，則相應的六維FPK方程求求解區(qū)域按50×50×50×50×50×50均勻網(wǎng)格離散化，網(wǎng)格數(shù)為506。區(qū)域網(wǎng)格數(shù)以幾何級數(shù)增長，會大大增加計算耗時。

因此在選擇濾波器數(shù)學模型時，應綜合考慮擬合精度與計算復雜度。低階濾波器建議選擇CARMA(2,1)模型，高階濾波器建議選擇CARMA(4,3)模型，下面將主要針對這2個濾波器進行海浪的模擬及驗證。

3 隨機海浪的仿真與驗證

濾波器生成譜與ITTC雙參數(shù)譜的擬合程度從理論上給出了利用CARMA模型構(gòu)造濾波器來對隨機海浪進行模擬的可行性。利用濾波器對波浪時歷進行數(shù)字模擬，并根據(jù)時歷得出的統(tǒng)計值與原波譜統(tǒng)計值進行對比，可以進一步探討不同階濾波器生成的隨機海浪與原隨機海浪的統(tǒng)計特性差異程度。

3.1 隨機海浪的數(shù)學模型

實際海面上的波浪主要由風產(chǎn)生，是一種非常復雜的隨機過程。由于風向的多變性和風的隨機性，風浪可能向各個方向傳播，起伏的波面就像形狀不等的小丘。在本文中，為使問題簡化，假定海浪只沿一個固定方向傳播，其波峰線和波谷線彼此平行且垂直于波浪的前進方向，這種海浪被稱作“二元不規(guī)則波”，或長峰波海浪。

在隨機海浪研究中，以Longuet-Higgins模型[16]在工程中的應用最為廣泛。Longuet-Higgins模型認為長峰波海浪可以看作由無數(shù)個不同波幅、頻率與初相位的余弦波疊加而成。固定點海浪波面位移可表示為：

(21)

式中ζAi、ωi、εi分別為第i個余弦波的波幅、頻率和初相位。其中εi為0～2π均勻分布的隨機相位。

振幅ζAi與譜之間的關(guān)系：

(22)

3.2 海浪的仿真與驗證

已知有義波高H1/3=2 m，特征周期T1=8 s的四級海況對應的ITTC雙參數(shù)譜，以及前文所推薦的CARMA(2,1)和CARMA(4,3)濾波器對應的有色噪聲功率譜，根據(jù)Longuet-Higgins模型生成波浪時歷，對海浪波面的瞬時值進行模擬。

由ITTC譜和CARMA(2,1)、CARMA(4,3)2種濾波器生成譜對應的波浪時歷如圖7～9所示。

圖7 基于ITTC譜生成的波浪時歷

圖8 基于CARMA(2,1)譜生成的波浪時歷

圖9 基于CARMA(4,3)譜生成的波浪時歷

表1為CARMA(2,1)和CARMA(4,3)2種濾波器所生成波浪時歷的統(tǒng)計值與ITTC譜生成的波浪時歷對應的統(tǒng)計值的對比。

表1 波浪時歷統(tǒng)計值的比較

由表1可以看出，基于自回歸滑動平均模型構(gòu)造的CARMA(2,1)二階濾波器所生成的波浪時歷統(tǒng)計值，比較接近ITTC雙參數(shù)譜對應的統(tǒng)計值。而由CARMA(4,3)四階濾波器生成的波浪時歷統(tǒng)計值更接近ITTC雙參數(shù)譜對應的統(tǒng)計值，對應的方差、有義波高、平均跨零周期的相對誤差分別是4%、1.5%、5.3%。

根據(jù)上面對功率譜和生成的波浪時歷統(tǒng)計值的比較結(jié)果，可以得出在精度要求不高的情況下，可采用CARMA(2,1)二階濾波器得到隨機波浪。而對于CARMA(4,3)濾波器來說，無論是從波譜曲線的形狀趨勢，還是所生成波浪時歷的統(tǒng)計值，該四階濾波器的功率譜都與ITTC雙參數(shù)譜有著更高的契合度。濾波器階數(shù)的提高雖然使其表達式更為復雜，但卻保證了對實際海浪更高精度的模擬效果。值得注意的是，在進行隨機海浪的數(shù)字仿真時，要同時考慮計算復雜度與模擬精度這兩方面的要求來選擇濾波器模型。

4 結(jié)論

1)二階濾波器結(jié)果表明，MA模型的引入可以極大地改善單純用AR模型對海浪譜的擬合程度。具有極零點譜形式的ARMA過程適用于海浪譜的模擬。

2)基于ARMA過程，可以構(gòu)造多種形式的濾波器模型。四階濾波器生成譜與ITTC雙參數(shù)譜的契合程度優(yōu)于二階濾波器，但是在實際問題的求解中應綜合考慮對模擬精度與計算復雜度的要求。

3)結(jié)合濾波器生成譜，利用Longuet-Higgins模型進行波浪時歷仿真，生成的時歷數(shù)據(jù)統(tǒng)計值與ITTC雙參數(shù)譜對應的時歷統(tǒng)計值之間的相似程度進一步驗證了由ARMA過程構(gòu)造的濾波器模型模擬隨機海浪的可行性。