基于隨機子空間法的海洋平臺模態特征實時提取方法研究

黃 焱，陳 濤，朱本瑞

(1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300350；2.天津大學 建筑工程學院，天津 300350)

海洋平臺作為海洋油氣資源開發的基礎性設施，是海上生產作業和生活的基地。為了保障海洋平臺結構的安全服役，避免重大的惡性事故發生，必須對海洋平臺結構在服役期內進行定期或不定期的檢測和安全評估。海洋平臺結構模態特征實時提取就是以此目標出發而衍生出的一項重要手段。目前，常用的模態特征提取方法包括峰值拾取法、自然激勵技術(Natural Excitation Technique,NET)、ITD(Ibrahim Time Domain,ITD)法、隨機子空間方法(Stochastic Subspace Identification,SSI)等。其中，SSI因其計算模型參數少，無迭代運算，收斂性好，運算快捷，魯棒性好等特點[1-3]，受到廣泛關注。

隨機子空間方法根據識別過程中相關性計算方式的不同，可以分為數據驅動隨機子空間方法和協方差驅動隨機子空間方法[4-5]。在同樣有限數據的情況下，兩種方法識別精度差別不大，但協方差驅動隨機子空間方法表現出更高的計算效率[6]。對于長期處于復雜激勵下的大型結構物，協方差驅動隨機子空間方法在結構模態實時提取中更具優勢。Peeters等[7]分析了隨機子空間方法的性質，并利用該方法識別受風載作用下的鋼桿的振動特性，取得了良好效果。樊可清等[8]推導了隨機子空間辨識算法中的協方差矩陣構造方式，通過減少辨識過程中的重復運算量達到縮短運算時間的目的，從而提高了方法的實用性，并通過對香港汀九大橋狀態監測系統的數據處理結果加以驗證。Reyndersat等[9]研究了隨機子空間識別模態參數的不確定性，提出用穩定圖方法剔除虛假模態，并通過一個橋梁仿真模型驗證了該方法在辨識密集模態時的有效性。

目前，隨機子空間方法多運用于陸上鋼結構、橋梁、建筑物等工程結構，而對于難以監測水下結構的海洋平臺涉獵甚少。辛俊峰等[10]針對一個四樁腿單層甲板海洋平臺，研究了冰激振動下的結構模態提取工作，并取得了較理想的效果。但其研究計算的結構比較簡單，環境工況的設置單一，忽略了環境條件的多變性及隨機性，缺乏對復雜工況下所帶來的虛假模態的考慮。針對復雜環境下的海上平臺而言，隨機子空間方法計算效率不高，有限測點下無法反映結構整體特征等問題尚未有很好的解決辦法，因此，如何能夠有效獲得海洋平臺等大型結構模態參數的研究仍有待深入。為此，本文基于海洋平臺的結構特點，提出了一種改進的協方差驅動隨機子空間模態實時提取方法，為結構健康評估提供參數化指標。

1 理論基礎

1.1 離散時間狀態空間模型

對于N個自由度的系統，在外界環境激勵下，其振動微分方程可用下述方程表示：

(1)

式中:M表示N×N維質量矩陣；C表示N×N維阻尼矩陣；K表示N×N維剛度矩陣；q(t)表示t時刻下的N維位移列向量；f(t)表示N維外部激勵向量。

實際測試中，考慮到信號采集的時間離散性和測量時存在的噪聲干擾現象，對式(1)進行采樣離散化處理以及微分求解后，得到系統離散時間狀態空間模型如下：

(2)

式中:yk和xk分別表示第k個采樣時刻系統的輸出量和狀態向量；wk和vk分別表示環境干擾噪聲序列和傳感器自身測量誤差所帶來的噪聲序列；A∈R2N×2N，表示系統狀態矩陣；C∈Rm×2N，表示系統輸出矩陣；N為系統階數，m是所采用的測點個數。

此空間模型是由狀態方程和輸出方程組合，構成一個系統完整的動態描述。其基本假設前提為：

(1)噪聲項為均值為零的平穩隨機過程，即

E[wk]=0E[vk]=0

(3)

(2)噪聲信號是均值為零的白噪聲序列，且與結構真實狀態無關，即

(4)

(3)系統為線性時不變系統(既滿足疊加原理又具有時不變特性)，狀態序列為平穩隨機過程，即

(5)

式中：E表示數學期望；狀態協方差矩陣Σ與時間k無關。為此，定義輸出協方差矩陣Ri如下：

(6)

式中：G表示狀態—輸出協方差矩陣，其表達式為：

(7)

在實際測試中，通常只選取一段數據量(S個有限數據點)參與計算，因此，輸出協方差矩陣可估計為：

(8)

1.2 協方差驅動隨機子空間方法計算流程

協方差驅動隨機子空間法的計算主要分為五個步驟：首先由響應數據y按反向對角線元素相等規則構成Hankel矩陣，并按矩陣的行列數將Hankel矩陣分為過去部分Yp與將來部分Yf；其次計算Yp與Yf的相關性；然后對計算后的矩陣進行奇異值分解(SVD分解)；再依據分解后結果獲得系統狀態矩陣A和輸出矩陣C；最后通過矩陣A和C求得結構模態參數。

(1)構造Hankel矩陣，并分為過去部分Yp與將來部分Yf[11]：

(9)

式中:yi∈Rm×1，表示第i時刻下所有測點(m個)的響應數據。

(2)構造T型矩陣(Toeplitz矩陣)，計算Yp與Yf的相關性，即

(10)

式中:Ri為輸出協方差矩陣。

將式(8)代入式(10)，并對T型矩陣分解可得到：

(11)

式中：表示協方差矩陣T1|i可表示為擴展可觀測矩陣Oi∈Rmi×2N與擴展可控矩陣Γi∈R2N×mi的乘積。

(3)T型矩陣奇異值分解(SVD)

在對T矩陣進行分解之前，通常先對其進行加權處理[12]，再對加權過后T型矩陣進行奇異值分解，可得：

(12)

式中：S1∈RN×N為對角矩陣，U、V∈Rmi×2N，是正交矩陣，即滿足UTU=VTV=I

(4)系統狀態矩陣A和輸出矩陣C

比較式(9)和(10)可知：

(13)

根據T型矩陣T1|i構造方式，構造下一時間延遲協方差矩陣T2|i+1

(14)

式中：矩陣T1|i和T2|i+1具有相同的結構，只是包含的協方差矩陣Rk時延變為從2到i+1。

觀察矩陣T1|i和T2|i+1可知，輸出矩陣C為拓展可觀測矩陣Oi的前m行，系統狀態矩陣A可用矩陣T2|i+1矩陣表示(*+表示廣義逆)：

(15)

(5)系統模態識別

對矩陣A進行特征值分解：

A=ψΛψ-1

(16)

式中:Λ=diag(λ)∈RN×N，是一個對角矩陣，該矩陣中的元素為系統的離散時間復特征值；ψ為離散時間特征向量組成的矩陣。

通過系統矩陣的特征值及特征向量可以得到系統模態參數為

(17)

式中:ξ為阻尼比，ω為圓頻率，Φ為振型。

隨機子空間方法在識別結構模態信息時，因輸入信號不滿足白噪聲的假定前提或輸出信號受到環境干擾，計算時容易出現虛假模態現象。且由于其算法采用狀態空間模型，需要預先人為確定系統階次n，而為了保證計算時避免出現模態遺漏問題，通常需要對系統階次進行過估計，因而，亦不可避免的產生虛假模態[13]。

對于虛假模態的消除，目前廣泛采用穩定圖[14]來實現。理由是虛假模態信息不是反映系統本身特征，是由于各種干擾產生，必然不如結構真實模態的計算結果穩定。穩定圖的做法是先選取一定常數據段，求不同系統階數下識別的模態結果反映到一個二維坐標系中(坐標圖的橫坐標為頻率值，縱坐標為階次)，最后呈現穩定的頻率信號值便被認定為結構的真實模態信息。目前，隨機子空間辨識法與穩定圖相結合的做法在工程應用上取得了較好的效果。

2 隨機子空間方法在海洋平臺結構的應用

根據隨機子空間方法基本原理可知，在識別過程中需要的人工干預少，算法本身抗噪能力強，識別結果的準確度高。但目前隨機子空間方法的研究多限于測點能夠全覆蓋的陸上橋梁等工程結構，對于海洋工程結構，該方法的應用存在諸多阻礙。主要原因在于海洋工程結構既無法做到測點全覆蓋布置，又時刻受到各種復雜海況的干擾，導致采集到的信號干擾成分過多并且很難保留結構完整信息，極大地影響了算法的識別效果。因此，隨機子空間方法在海洋工程結構上的應用仍有待進一步研究。

本節以渤海某導管架平臺為研究對象，來分析隨機子空間方法在海洋平臺模態參數的識別效果。該平臺由六樁腿的導管架鋼結構作為支撐，上部組塊分為3層，作業水深12.6 m，平臺總質量為9 610.4 t，見圖1。采用有限元軟件ANSYS建立該平臺結構有限元模型，見圖2。建模時對井口區隔水套管進行模擬，并充分考慮平臺歷次改造的結構與設備布置變化，保證了各層甲板重心以及模型的總重量重心與平臺設計及改造資料一致，使模型在各方面盡可能接近平臺結構的真實狀態。采用隨機波浪理論模擬現實海況，對模型平臺施加隨機波浪激勵，模擬各種海況下平臺所受的波浪沖擊。根據該平臺所處海域水文資料統計結果，90°波浪入射方向為該平臺的主浪向，其它出現概率較大的波浪方向為0°、30°，波浪入射方向見圖3。為盡量覆蓋平臺可能遭遇的波浪環境，針對這3種波浪入射方向均設置了3種出現概率最大的波高條件，有效波高分別為0.5 m，1.0 m，1.5 m，不同工況對應編號,如表1所示。

圖1 渤海某導管架平臺示意圖

圖2 平臺有限元模型及測點位置示意圖

圖3 平臺水線面處波浪入射方向示意圖

表1 有限元模型計算工況及其編號

現實監測中，考慮到傳感器的適用性與施工可行性，測點通常布置在平臺上部甲板上。本文選取的目標平臺共有三層主甲板六根樁腿，按照一般的監測原則，每層甲板可布置6個測點，共計18個測點。對應于現實結構，數模計算時，選取平臺結構上部三層甲板與六根樁腿相交的18個節點作為加速度信號采集點，即測點數目m=18，采樣頻率為10 Hz，采樣時間為9 000 s。

通過ANSYS模態分析得出平臺前三階頻率信息，其中包括x方向、y方向以及扭轉方向(z)的模態。如表2所示。

表2 導管架平臺有限元模型模態信息

面對目標平臺如此復雜的結構，經過前期的估算，發現至少要形成18個測點的布置方案，才有可能對系統形成較好的識別。針對該布置方案，采用傳統的隨機子空間方法的處理方式是：選用監測結構響應的所有18個測點按順序組成一列向量構建Hankel矩陣，再利用隨機子空間算法進行模態識別，最后結合穩定圖提取結構穩定模態。九種工況下系統的識別結果如圖4所示。需要說明的是，由于海洋工程結構的復雜性與測點布置的局限性，可以預判此時所形成的模態識別結果中必然會含有一定數量的虛假模態，因此，為與上面給出的結構真實模態計算結果相區分，將這些識別結果稱作準模態。考慮到均值可以表征準確度，而方差則反映了穩定性水平，因此識別結果中分別給出了每組數據的均值μ與方差σ，如表3所示。

圖4 傳統方法不同工況穩定圖

如表3數據所示，其中準五階，準六階，準七階頻率與結構前三階真實頻率信息有很好的對標性，但其余的準一至四階與準八階模態信息并不是結構本身所擁有的，即傳統方法下出現了無法消除的虛假模態信息。

由表3中數據和圖4可知，在各種工況下，傳統隨機子空間法不僅在真實模態頻率上出現了穩定點，在許多非系統模態頻率上也表現出穩定點集中現象。且穩定圖中數據點波動較大，即各種工況下的識別結果存在較大波動。這主要是由于測點無法實現結構全覆蓋，導致算法對結構信息識別的“發散式”推測，以及對結構系統階數的過估計所致。

表3 不同工況頻率識別結果

由此可見，傳統方法往往會導致虛假的模態提取。因此，如何基于隨機子空間方法基本原理，發展一套能在有限測點下增強反映結構整體特征，并且有效消除虛假模態的改進方法仍需進一步討論。

3 改進隨機子空間方法

海洋平臺工作環境的限制，決定了其不可能和陸地結構物一樣盡可能采用測點全覆蓋的方式來進行監測信號獲取。因此，如何根據平臺水上結構有限測點數據信息推測整體模態特征是隨機子空間法運用的關鍵技術難點。由此目標出發，本文提出了兩點改進：Hankel矩陣元素重構，即每個元素重新嵌套子矩陣的方法和基于信號累積特征的虛假模態辨識方法與剔除準則。

3.1 Hankel矩陣元素重構

隨機子空間方法的核心是Hankel矩陣的建立，其描述的是動態響應的時間連續繼承關系，在數學上體現為過往響應與未來響應間的相關性。因此，當Hankel矩陣中描述的相關性僅是有限的結構片段時，就會在預先設定的模態階數所控制的滿秩矩陣運算中，產生對結構剩余部分的“發散式”推測。由隨機子空間方法的基礎理論可知，數字矩陣間的協方差關系是形成結構信息描述的基本手段。那么可以推知，如果將局部位置(測點)的數據作為獨立的Hankel矩陣元素進行運算，實際上是抹除了局部位置間本來可能具備的關聯性而進行全新的關聯性重構，這也是造成“發散式”推測的根本原因。因此，對于現實中監測數據僅為結構的片段信息時，就必須在Hankel矩陣的構建中盡可能保留甚至增強反映結構整體特征的數字信息量，即以某種測點數據排列方式構建一個關聯子矩陣作為Hankel矩陣運算的基本元素。

在結構有限測點的情況下，為了在Hankel矩陣中盡可能保留結構整體特征的數字信息量，可以在選擇測點數據排列方式時盡可能的對標測點在實際結構中的空間排布，這里主要考慮兩方面對標性：

(1)子矩陣元素排布對標實際測點空間排布：Hankel矩陣中每個元素yi表示在第i個時刻下各個測點響應信號值，那么只需要把各個測點在某個平面內的空間排布位置與子矩陣中各元素位置對應，這些測點便可以按一定的排列方式組成一個子矩陣，用來取代Hankel矩陣中對應元素。為了能夠充分收集到結構各個方向的模態信息，可以選取結構各個橫剖面和縱剖面內的測點組成多種子矩陣交替計算形式。為了能更清楚表明結構測點位置信息，將上述導管架平臺測點位置表示成如圖5的形式。

圖5 測點位置示意圖

工程實際中，不同的測點位置對結構整體信息敏感程度也不一樣，如果僅選取一個面內的測點響應數據不足以表達結構的整體信息。為了盡可能方便提取結構所有方向模態信息，把能表達結構不同方向信息的橫剖面和縱剖面測點數據分別構成多個子矩陣形式，具體如圖6所示。圖中，yki表示第k個傳感器在第i時刻下的響應數據。

圖6 子矩陣形式

(2)子矩陣元素之間關聯性的保持對標響應數據獲取方向：按照上述子矩陣構建方式，已經完成了測點位置信息的表達，但子矩陣內部仍然缺少元素方向關聯性(運算時各子矩陣內沿同一方向振動的測點數據不被打散)的表達。為了補全這種方向關聯性，子矩陣內部元素還需進行進一步的向量重構，即把子矩陣元素沿振動方法合并成多個行向量。所以子矩陣構建可進一步表達為：

(18)

其中:

(19)

(20)

或

(21)

可以看出，相比較傳統隨機子空間方法把所有m個測點信號作為一個列向量輸入，改進后的Hankel矩陣有明顯的兩個優勢：子矩陣中盡可能的包含了結構整體特征的數字信息，使得計算結果更能反映結構的真實模態信息。另外，新的輸出協方差矩陣由原來的m×m維變為現在的m/t×(m/t)維(t為Hankel矩陣子矩陣列數)，極大地縮減了后續矩陣運算維度，明顯的提高了隨機子空間方法的計算效率。

為驗證元素重構后的識別效果，將上述五種子矩陣組合的構造方式運用于模型數據計算，分別采用隨機子空間算法，對各自識別結果進行穩定圖處理，得到五組穩定圖數據。將這些穩定圖數據進行均化合并，得到一組綜合穩定圖，九種工況下識別結果,如表4所示。

表4 不同工況頻率識別結果

對比表3和表4可知，相比較傳統方法識別結果，Hankel矩陣元素重構后，表現出多方面的優點。①改進后方法中每一階次數據的方差下降十倍以上，識別結果數據穩定程度大幅提升。②改進后方法識別的準五階，準六階，準七階頻率信息，其結果也更加接近結構的前三階模態真實值，增加了識別結果的準確度。③經過Hankel矩陣元素重構，新的Toeplitz矩陣由原來的mi×mi維變為現在的mi/t×mi/t維(t為Hankel矩陣子矩陣列數)，即新的協方差矩陣維度縮減了3×3倍，極大地縮減了后續矩陣運算維度。在相同計算條件下(戴爾DELL T630塔式服務器，32 GB內存，雙英特爾處理器)對兩種算法分別完成一組模態識別所需的計算時間進行比較，結果表明：傳統隨機子空間算法用時34.8 min，而改進算法僅用10.5 min，計算效率提高了三倍。綜合而言，通過Hankel矩陣元素重構后，改進算法既提高了模態識別的準確度與穩定性，又極大地的提升了計算效率，具有良好的優越性。

應注意的是，不同響應數據經過計算處理后，很多非系統真實特征的信息在穩定圖上亦呈現出很穩定的結果，如表4中的準一階，準二階，準三階，準四階，準八階模態。此外，在識別過程中，也很容易使得結構真實模態信息在不同的系統階數下表現為不同的階次，即出現模態越階現象。那么這時候，傳統的穩定圖方法會把結構真實模態信息消除，從而出現誤判。為此，必須引入一種新的能夠有效消除虛假模態的準則，來輔助提取結構真實特征信息。

3.2 基于信號累積特征的虛假模態辨識方法與剔除準則

由于虛假模態的出現，是由于計算過程中存在各種干擾，從而引起的對結構模態信息的發散性推測。它反映的不是系統本身特征，且不同時段下結構所受到的干擾不同，因此，隨著時間的推移，虛假模態的不穩定現象會更加突出，而真實模態則具有模態魯棒性的時間累積增強效應。因此，可以通過采用不同時間段的數據計算結果來統計結構信息的穩定程度，從而甄選出結構的真實模態信息。

基于此原理，本文提出一種基于信號累積特征的虛假模態辨識方法與剔除準則，主要是通過不同時段數據結合不同子矩陣的形式依次做隨機子空間方法計算，最后把識別結果在各個小區間下做分段統計。那么，結構的真實模態必然會在一個區間范圍內保持較好穩定性，而虛假模態則會相對雜亂的排布在其他小區間，其表現出來的穩定性自然也會越差。其基本步驟如下：

步驟1：選取一較長時間段數據S并按數據長度均分成多段數據Si。

步驟2：選取結構不同截面測點組成j種子矩陣構建形式，并把把每段數據Si有時間交叉的分割成j小段數據Sij(數據Sij間有80%重疊率，使得結構真實模態信號更穩定而虛假模態產生較大變化，數據Sij表示在數據Si下的第j小段數據)。

步驟3：每一段數據Sij對應一種子矩陣規則進行計算，得到j種Hankel矩陣，依次做子空間運算識別結構模態，并把整體數據段S下的頻率計算結果在各個區間下做分段統計，即每組識別結果按數值大小劃分在不同區間(為使各區間足夠把結構各階次模態區分開，定義每個區間的長度為0.02 Hz)。

步驟4：把包含數據量少的區間作為不穩定區間(定義每個區間內包含P種子矩陣的識別結果以及每個區間內包含Q個數據，取P、Q閾值為p、q，則各區間數據結果小于p、q的定義為不穩定區間)，拋除不穩定區間，留下穩定區間。

步驟5：統計穩定區間內前幾組連續區間，其中每組連續區間內數據期望值便作為結構的真實模態信息模態。

針對文中算例，采用上述步驟進行虛假模態提取，其流程圖如圖7所示。取數據S為9 000 s，均分成五組S1，S2，S3，S4，S5，并把每組數據Si在80%數據重疊率下分成五段數據Sij，每一段數據Sij各選取一種子矩陣規則進行運算，交替做隨機子空間運算，識別出多組結構模態。最后把整體數據S下的頻率計算結果在各個區間下做分段統計，把包含數據量少的區間作為不穩定區間。定義不穩定分段時，經過多次驗證，P、Q閾值分別取4和6。拋除不穩定區間，留下穩定區間，其中每組連續穩定區間內數據期望值便作為結構的真實模態信息。以各連續穩定區間內數據數學期望值為橫坐標，各連續穩定區間內數據總量為縱坐標，統計各工況模態計算結果，如圖8所示。

圖7 模態提取流程圖

圖8 不同工況頻率識別結果

依據上述方法，對結構9種工況下響應數據識別結果進行分析。得出目標平臺的前三階頻率識別結果。如表5和圖9所示。

觀察表5和圖9可以發現：經過Hankel矩陣重構與虛假模態消除處理后，最終識別出的準一階，準二階，準三階模態與結構前三階真實模態有很好的對標性，有效剔除了傳統方法下的眾多虛假模態，識別結果準確度也保持在95%以上。對于有限測點下的海洋平臺結構，改進后的隨機子空間方法有很好的優越性。

圖9 不同工況頻率識別結果

表5 不同工況頻率識別結果

(1)改進后的隨機子空間算法能有效識別結構各階次模態結果，并把識別結果集中在一定動態范圍之內，有效排除了眾多虛假模態的干擾，準確提取出系統模態所在區間。

(2)相對于Hankel矩陣元素重構后的識別結果而言(不考慮虛假模態的情況下)，穩定段內識別出的模態信息與系統真實模態非常接近，而且識別結果的波動也比較平穩。也就是說，在不損失計算結果穩定程度的基礎上，識別結果準確度有了明顯的提升。

總而言之，改進后的隨機子空間方法識別結果與理論分析保持很好的一致性。采用子矩陣與信號累計的形式，使計算結果盡可能的包含了結構整體特征的數字信息，極大地縮減了矩陣運算維度，明顯的提高了隨機子空間方法的識別準確度和計算效率。相比原重復性計算同一組數據下做穩定圖，能更加明顯的排除虛假模態干擾以及有效提取結構真實模態信息。

4 結 論

針對海洋平臺模態識別問題，本文在傳統隨機子空間方法理論基礎上，提出了新的Hankel矩陣元素重構方式以及虛假模態辨識方法與剔除準則。為了更好反映結構各個方向的整體數字特征信息，在原有Hankel矩陣元素上重新嵌套新的子矩陣。對標測點實際空間位置建立多種子矩陣組合方式，并結合有時間重疊的響應數據做交替計算。為消除識別結果中的眾多虛假模態以及修正各組次下相近模態計算的越階現象，把所有識別結果采用各區間分段統計的方式來處理。通過對目標平臺有限元模型響應信號進行計算，在有限測點下情況下，改進的方法能更加明顯的排除虛假模態干擾并有效提取結構前三階真實模態信息，也極大地提高了原有算法計算效率，具有良好的適用性和魯棒性。