葉靜嫻，呂中榮，汪利

(中山大學航空航天學院，廣東 廣州 510000)

海洋資源和能源越來越受到人們的重視[1-2]。與之相對應的海上結構，如海上采油平臺、風力發(fā)電機等，因其所處的環(huán)境相當復雜，往往同時受到波浪荷載、風荷載、冰荷載等荷載的共同作用[3]。其中，波浪荷載是最主要的荷載[4]?？紤]到波浪荷載的隨機性和結構所處海洋環(huán)境的復雜性，波浪荷載難以直接實時測量，這將對海洋結構的設計和狀態(tài)評估帶來困擾。目前，結構響應的測量手段很成熟，如：加速度傳感器可以直接提供加速度信息，應變片可以直接測量結構的應變。因此，可以海洋結構的實測響應為基礎，反演得到波浪荷載。這一工作不僅能為海上結構的設計提供必要的荷載信息，同時也能為海上結構的運行監(jiān)測和疲勞監(jiān)測等創(chuàng)造條件。

根據響應反演結構所受荷載，是典型的反問題。國內外學者提出了一些波浪荷載的反演方法。Maes 等[5]基于一個風力渦輪機相似模型的加速度時程和動態(tài)模型，以遞歸聯(lián)合輸入狀態(tài)估計算法反演了作用在其上的破波荷載、以及其他荷載，并進行了實驗驗證。Jensen等[6]利用ARMA校正模型，針對一個墩柱系統(tǒng)，以加速度作為輸入，進行了系統(tǒng)參數識別和隨機波浪荷載反演，并通過實驗驗證了該方法的有效性。王子健[7]第一次用全量補償復合反演算法進行了橋墩的參數識別與波浪荷載反演，并結合概率論中的統(tǒng)計平均法與消去法，提出了改進的復合反演算法。以上方法大多將波浪荷載反演問題轉化成一個狀態(tài)估計問題，識別過程中將波浪荷載視為狀態(tài)變量，然后采用濾波類方法進行狀態(tài)變量估計。該類方法對于隨時間變化緩慢的荷載識別問題較為有效，但如果荷載隨時間變化較為劇烈，如：隨機波浪荷載，該類方法將會產生較大的識別誤差。本文提出了一種基于逐步優(yōu)化算法的隨機波浪荷載反演方法，該方法針對每兩個相鄰時間點建立目標函數，以加速度時程為輸入，逐步識別一系列時間點上的波浪荷載，計算代價較低。該方法將荷載當作外部輸入，而非狀態(tài)變量，因此對時間變化劇烈的荷載如波浪荷載，仍能得到滿意的識別結果。

1 基于波浪譜生成隨機波浪荷載

Longuet-Higgins模型是常用的描述隨機波浪的模型[8]。該模型認為，固定點的波面位置η(t)可看做是由許多不同振幅、不同周期和不同隨機初位相的諧波迭加而成。

(1)

其中，t為時間；an為各組成波的振幅；ωn為圓頻率；εn為隨機初相位，它是均布于[0,2π]的隨機變量?？紤]到波面位置的隨機性，通常使用功率譜密度函數Sηη(ω)描述波面位置的統(tǒng)計信息，即稱為波面譜。本文考慮一種常用的P-M不規(guī)則波面譜[9]，其表達式為：

(2)

其中，HS為有義波高?；诓ɡ撕奢d與波面位置的相關性，推導波浪荷載與波面譜的關系。考慮如圖1所示小尺寸樁柱所受的波浪荷載，以海底為坐標原點，建立如下所示坐標系，根據Morison方程，z處的波浪荷載為：

f(z,t)=fD(z,t)+fM(z,t)

(3)

圖1 圓管柱所受隨機波浪荷載Fig.1 Random wave load on small pile column

根據波面函數(1)，由線性波浪理論可得水質點速度和加速度的表達式如下：

(4)

其中，Tn為周期，滿足Tn=2π/ωn，d為水深，kn為波數，滿足以下色散關系：

ωn=gkntanhknd

(5)

進一步把Morison方程中的阻力項線性化處理，可得：

(6)

其中，σu(z)為速度的均方差，表達式如下：

(7)

由于k是ω的函數，它們滿足式(5)中的關系，故式(7)不能通過直接積分求得。為簡化起見，僅對上式波譜部分做積分，結合P-M譜，可得：

(8)

波浪荷載在柱高方向呈現(xiàn)出空間相關性，由式(3)和(4)可求得t時刻波浪荷載沿水深方向的分布，故可定出相應的空間分布系數?？紤]水面下z1至z2處所受的合力，將波浪荷載沿著結構方向積分即可，再乘上相應的空間分布系數矩陣，即可得到各個自由度上的波浪荷載。通過積分得到總速度合力FD(t)和總慣性合力FM(t)如下：

(9)

以z1=0和z2=d代入上式進行積分, 并考慮色散關系以及σu的近似計算，由譜分析可得總速度力譜SFD(ω)和總慣性力譜SFM(ω)[10]。

(10)

相加即得總波浪荷載譜：

SFF(ω)=SFD(ω)+SFM(ω)

(11)

根據平穩(wěn)隨機過程的諧波合成法，波浪荷載F(t)可按如下方式生成[11]：

(12)

(1)劃分頻率區(qū)間的方法有等分能量法和等分頻率法。在多數情況下宜采用前者，但對結構進行動力分析宜用后者。則上式中，Δω=(ωm+1-ω1)/m，m一般可取50～100，本文取50。

(3)頻率范圍ω1～ωm+1的選取, 在高頻率和低頻率段各舍去總能量的μ部分，這里μ取0.002。對于P-M譜有：

(13)

2 時程分析

多自由度系統(tǒng)的運動方程為：

(14)

(15)

其中， Δt是時間步長，參數γ和β在本文中分別取0.50和0.25。將式(15)代入(14)可得：

(16)

通過式(16)，可由ti時刻的響應求得ti+1時刻的位移，將位移解代入(15)可求得ti+1時刻的加速度和速度如下：

(17)

循環(huán)以上步驟，可以求得所有離散時間點的位移、速度和加速度。

3 基于逐步優(yōu)化算法的荷載反演方法

重新整合公式(16)和(17)，可得狀態(tài)方程如下：

xi+1=Hxi+G1Pi+1+G2Pi

(18)

(19)

其中，I為相應維度的單位矩陣。系統(tǒng)的觀測方程如下：

yi+1=Lxi+1+M-1Pi+1

(20)

其中，yi+1表示ti+1時刻的測量數據，L為與實測結構響應相關的矩陣。本文主要考慮加速度測量數據，由式(14)可得L=[-M-1Κ,-M-1C]。

(21)

其中，λ是權值，其取值與測量誤差水平有關。一般地，測量數據的誤差越小，λ應越大；反之，λ則應越小。目標函數(21)是一個二次函數，其優(yōu)化問題可直接求解，可得如下結果：

(22)

其中，矩陣A、b的表達式如下：

(23)

根據迭代式(23)，最終可以逐個時間點識別波浪荷載?？紤]波浪荷載具有一定的空間分布規(guī)律，將波浪荷載向量Pi用總波浪荷載Fi與空間分布系數矩陣Q的乘積來表示，即Pi=QFi。此時，迭代計算須用如式(24)-(25)替代。

(24)

(25)

在波浪荷載反演過程中，考慮波浪荷載的空間分布，能有效地縮減識別參數的個數，從而提高識別的精度。

4 算例驗證

由于測量不可避免地存在誤差，為了考察誤差對反演算法的影響，在用Newmark-β法求得結構加速度時程后，添加一定水平的高斯白噪聲，如下所示。

(26)

將含有不同程度噪聲干擾的加速度時程作為輸入，采用本文提出的荷載反演算法進行反演，識別出了兩個自由度的荷載時程、速度時程以及位移時程，并進行誤差分析。圖3給出了無噪聲情況下的反演結果?？梢钥闯觯寒敓o噪聲干擾時，波浪荷載的測量值與識別值幾乎完全重合，這說明該反演算法在無噪聲情況下是切實可行的。

圖2 固定式海洋平臺及其簡化模型Fig.2 Fixed offshore platform and its’ simplified model

圖4-5分別是噪聲水平為5%和10%情況下的識別結果。由圖4可見，在5%的噪聲水平下，波浪荷載的識別值與測量值吻合很好，剔除幾個極值點，其余點的識別結果的相對誤差均小于1%。 由圖5可見，在10% 的噪聲水平下，剔除幾個極值點，其余點的識別結果的相對誤差均小于5%，滿足工程需要。因此，該反演算法具有較強的抗噪聲能力。

圖3 無噪聲情況下第二個自由度波浪荷載識別結果Fig.3 Wave force identified results of the second degree of freedom in the case of no noise interference

圖4 噪聲水平為5%時第二個自由度波浪荷載識別結果Fig.4 Wave force identified results of the second degree of freedom in the case of 5% noise level

圖5 噪聲水平為10%時第二個自由度波浪荷載識別結果Fig.5 Wave force identified results of the second degree of freedom in the case of 10% noise level

5 結 論

本文提出了一種波浪荷載反演的逐步優(yōu)化算法，該算法通過Newmark-β數值積分法將動力問題表達成狀態(tài)方程，然后再基于每一個時間點建立目標函數，逐步識別波浪荷載。數值算例表明，即使在10%的噪聲水平下，反演得到的波浪荷載仍具有很高的精度。