一種基于狀態空間模型的結構振動模態軌跡提取方法

周 胡，錢忠文，王 濱，高樹健，劉福順

（1.浙江省深遠海風電技術研究重點實驗室，杭州311122；2.中國電建集團 華東勘測設計研究院有限公司，杭州 311122；3.中國國電集團有限公司 浙江分公司，杭州 310000；4.中國海洋大學 山東省海洋工程重點實驗室，山東 青島 266100）

近年來，由于海上風能的穩定性和可靠性，海上風力發電機已經得到了持續的發展。自2010年以來，裝機容量以每年30%的平均速率不斷增長。然而，由于復雜的海洋環境作用，海上風力結構極易發生損壞，在過去的十年中，由于風機結構的損壞導致了多起事故的發生，僅在2014年就發生了約1 000起事故。由于高昂的安裝和維護成本，對海上風力發電機的結構可靠性進行評估已經變得至關重要[1]。

海上風力發電機的軌跡是指其轉子軸心相對于風機中心的運動軌跡，包含有風機運行的故障征兆信息，可以直觀的反應風機的瞬時運動狀態，是評估風機運行的重要狀態特征參量。目前，對軸心軌跡的研究主要包括軸心軌跡提純和辨識法兩大類[2]。由于實際的振動信號會受到噪聲的干擾，軸心軌跡十分復雜，難以從中獲得清晰的特征，因此需要對軸心軌跡進行提純，從而還原結構的軸心軌跡。軸心軌跡的辨識則是研究軸心軌跡的特征提取和識別。目前對于軸心軌跡的提純和辨識主要是基于數字濾波、短時傅里葉變化和小波變換等[3]。這些方法都是采用一定的頻段劃分后，對原始信號進行濾波提純[4]。但是在基于傅里葉變換的方法中，由于實測信號往往不滿足周期性假設的前提，會導致分解結果不準確。而在實際應用中，不準確的分解將直接影響軸心軌跡的形狀，使軸心軌跡在工程中的應用價值被大大削弱[5]。Prony方法將信號表示為復指數形式，避免了信號的周期性假設，實現了信號的高精度分解[6]。但是其對噪聲極為敏感，在應用于實測信號時分解結果存在較大的誤差[7]。Hu等[8]通過引入低階狀態空間模型，將高階微分方程轉換為低級矩陣微分方程，避免了傳統Prony方法多項式求根的病態問題，顯著提高了Prony方法的抗噪性能，實現了其在工程實測信號中的應用。

除了軸心軌跡的分析方法外，海上風電結構軸心軌跡計算中還存在另一個難點，就是在海洋環境中，很難找到固定的參考點固定位移計[9]，而使用GPS測量位移時，除了會帶來高額的測量費用，其測量精度較低，僅能達到1 cm，也會給計算帶來誤差[10]。由于加速度的測量方便快捷且精度較高，已經被廣泛地應用于海上風電結構的振動測試。雖然理論上可以對加速度積分得到位移，但是由于初始條件的未知，積分得到的位移會發生不真實的漂移[11]。目前解決此類問題的方法主要分為基于頻域積分的方法[12]和基于高通濾波調整基線的方法[13]。基于頻域積分的方法由于Fourier變換的周期性假設，在處理實測信號時往往會導致結果的不準確。雖然基于高通濾波調整基線的方法已經得到了廣泛的應用，但是對于低頻為主導的海上風電結構，高通濾波除了會濾出漂移項，還會改變振動信號中的低頻項，產生不必要的誤差[14]。

針對以上兩個問題，本文根據實測海上風機振動信號的特點，通過使用基于低階狀態空間模型的復指數分解方法，發展了一種海上風力發電機模態軌跡的計算分析方法，并分別使用數值算例和海上風機實測振動數據驗證了發展方法的有效性。首先使用五自由度數值系統驗證了由加速度信號獲得位移中模態信息的可行性，然后通過實測的單樁海上風力發電機驗證了所提出方法在處理實測信號時的有效性。

1 基本原理

1.1 復指數序列

對于非周期信號x（t），其復指數序列可以表達為

式中:n為級數個數；γ和λ分別為各成分的復指數序列的參數，分別可以表示為

1.2 Prony方法

Prony方法使用復指數序列擬合信號，可以概括為以下三個步驟:

步驟1使用高階微分方程擬合數據點求解系數

式中:K為采樣點數，a0，a1，…，an-1為方程組的系數。

步驟2在由步驟1求解得到多項式系數后，計算方程組的特征根（表示為zi（i=1，2，…，n）），從而求解復指數參數λi，并通過參數求解各成分的頻率和衰減因子。

步驟3將求解得到的特征根表示為如下形式，之后可以通過最小二乘法求解復指數參數γi

2 基于低階狀態空間模型的模態軌跡分析

2.1 加速度信號的復指數表達

在同一平面內對風機結構的x向和y向加速度信號進行測試，由于在實測風機信號中，不僅僅包含系統的結構信息，同時還包含有噪聲信號、機械振動信號等。結構信息中包含了結構的模態頻率和阻尼因子，噪聲成分通常較為復雜，包括由機械運轉和其他因素引起的非衰減成分。下面以x向信號為例進行算法的展示。首先，采用復指數序列式（1）的形式對測試信號x（t）進行表征

式中，n為信號中復指數序列中成分的數量。

式（5）可以視為系數為bi的Nn階常微分方程的通解，即

式（6）的特征多項式可以表示為

式（7）的根即為 λi（i=1，2，3，…，n）

2.2 基于低階狀態空間模型的模態分解

引入

假設bi=1，式（3）可以轉換為一階矩陣微分方程

為了求解式（9），引入Hankel矩陣

式中，ξ和η分別為Hankel矩陣的行數和列數。

分別令式（6）中的k=0和k=1，并對H（0）和H（1）進行奇異值分解，可以得到

式中:S1為奇異值矩陣；U1和V1為正交矩陣；A為狀態矩陣，其特征根為zn。通過特征根zn，就可以計算得到對應的復指數參數λi=ln zn/Δt。同時由式（13）可以通過最小二乘法求得γi

式中，N為采樣點數。更詳細的求解過程可以參考包興先的研究。

2.3 模態軌跡重構

通過求得的γ和λ，即可以重構測試得到的x向加速度信號中的模態成分

通過對式（14）進行二次積分，可以得到x向位移中對應模態的時域曲線

同樣，可以對y向加速度進行上述過程可以得到y（t）中相應模態的時域曲線

3 五自由度數值算例

在本數值算例中，將使用五自由度數值模型驗證所提出方法的正確性。該模型如圖1所示，其中mi=60 kg，ki=3×107N/m，ci=600 Ns/m（i=1，2，…，5）。

圖1 五自由度數值模型Fig.1 The five-degree-of freedom system

根據系統的特征方程可以得到，該系統的模態頻率分別為32.032 Hz，93.501 Hz，147.395 Hz，189.349 Hz和215.962 Hz；其五階模態阻尼比分別為0.002 012 6，0.005 874 9，0.009 261 1，0.011 897 0和0.013 569 0。

首先，構造一個數值信號作為系統的輸入外荷載

式（17）中各參數如表1所示，合成的外荷載如圖2所示。將該合成的外荷載作用于五自由度系統的前三個自由度上，即可以得到系統的動力響應。

圖2 合成外荷載時域曲線Fig.2 The simulated external force in the time domain

表1 式（17）中外荷載各成分的參數Tab.1 The parameters used in Eq.17

為了獲得該五自由度系統的動力響應，使用Newmark-β進行計算。 如圖3（a）和圖3（b）分別展示了使用Newmark-β方法計算得到的系統第5個自由度的加速度和位移的時域結果。在計算中，設置時間間隔Δt=0.01 s，計算點數N=1 000。

圖3 計算得到的加速度和位移Fig.3 The calculated accelerations and displacements

在使用本文方法進行結構模態軌跡求解時，首先需要對加速度進行復指數分解。分解時首先需要假設信號中的模態階次，即振動加速度中假設的信號成分的數量的兩倍。在分解時，由于荷載中共包含五個成分，且該系統共有五個自由度，因此其振動信號中應該有十個頻率成分，所以這里設置分解的模態階次為20。

為了檢查分解模態階次的設置是否合適，可以通過對分解后的信號進行重構，并將重構結構與原始信號進行對比來判斷，重構結果如圖4所示。從圖中可以看到，分解結果與重構結果吻合較好，這也證明了分解結果的正確性。從加速度的分解結果中選擇對應于結構一階模態的γ和λ進行重構，同時，可以利用式（16）對位移中的一階模態信息進行重構，重構結果如圖5所示。

圖4 分解重構結果Fig.4 The reconstruction result after decomposition

同時，為了驗證重構結果的正確性，同樣使用本文提到的基于低階狀態空間模型的復指數分解方法對Newmark-β計算的結構位移進行分解，并提取出其中的一階模態，結果如圖5所示。可以看到，使用加速度重構得到的一階模態位移信息與位移中的位移信息吻合較好，這證明了本文方法中通過加速度信息求解一階模態位移信息的正確性，也意味著在結構位移無法測量的情況下，可以使用本文方法對結構位移中的模態信息進行重構。

圖5 位移中的一階模態時域圖和由加速度重構的一階位移模態Fig.5 The first-mode in calculated displacements and the reconstructed first-mode displacements using accelerations

4 單樁式海上風電結構實測數據分析

為了驗證本文提出的模態軌跡分析方法在處理實測數據時的有效性，使用某單樁式海上風電結構實測數據進行模態軌跡分析。

4.1 現場測試情況介紹

測試風機結構為某4 MW單樁式海上風力發電機，位于江蘇省如東縣的黃海區域，塔筒高度為86 m，總質量為2 500 t。該風機結構的基礎為直徑5.5 m的單樁式結構，固定于在床上，海床下的穿透深度為50.35 m，從海床表面到法蘭的高度為22.65 m。該風機共具有三個葉片，每個葉片長度為88 m。在測試中，在風機結構內壁自下而上安裝五個加速度傳感器，分別記為1～5號，以法蘭為零點，五個加速度傳感器的高度分別為0，16.3 m，42.7 m，59.4 m和73.4 m。 該風機結構的一階設計頻率為0.305 6 Hz，二階設計頻率為1.729 6 Hz。測試的單樁式海上風機結構如圖6所示。

圖6 測試的風機結構Fig.6 The test offshore wind turbine

測試時，使用的三向加速度傳感器的型號為4630-005-060-C，其測量范圍為0～5g，電壓靈敏度為400 mV/g。測試中，同時對三個方向的振動加速度信號進行采集，所使用的數據采集儀為IMC-Access公司生產CRONOSflex400。測試中使用的加速度傳感器及采集儀，如圖7所示。

進行測試時，設置采樣頻率為200 Hz，采樣時間設置為1 h，5號傳感器記錄得到的x向和y向加速度信號時域圖，如圖8所示。從時域圖中可以看出，由于風荷載和波浪荷載隨時間的變化，加速度振動的幅值在不同時間發生變化。

圖8 5號傳感器測試的x向和y向加速度信號Fig.8 The measured accelerations by the 5th acceleration sensors in x direction and y direction

理論上結構的軸心軌跡可以通過結構平面內同時刻的x向和y向位移數據直接獲得，但是其前提條件必須首先獲得結構的位移。但是在結構的位移測量中，需要有固定的參考點固定位移計。但是在海洋環境中進行測試時，很難找到固定的參考點，因此目前的海上結構物的振動測試，都是進行加速度測試，然后通過使用測試的加速度數據轉換成位移數據。目前已經有很多的研究去除加速度積分后的積分項來獲得結構的真實位移數據。大多數的加速度轉位移的方法都是基于Trifunac等提出的方法發展期出來的，因此本文首先使用Trifunac等提出的方法對測試得到的加速度信號進行分析從而重構位移信號。對圖7中1 000～1 050 s的加速度信號進行分析，在計算中，設置高通濾波中的截止頻率為0.1 Hz，轉換的位移結果如圖9所示。

圖7 測試中使用的加速度傳感器和采集儀Fig.7 The used acceleration sensors and IMC collectors in the measurement

圖10繪制了由圖9中重構位移直接得到的結構在1 000～1 050 s時的軸心軌跡圖，可以看到由于受到復雜的環境噪聲干擾，結構的運動狀況較為復雜，軸心軌跡非常混亂，無法對結構的運行狀態進行正確的監測和檢測。

圖9 由測試的加速度信號轉為的位移信號x向和y向Fig.9 The displacements reconstructed by measured accelerations in x direction and y direction

圖10 由重構的位移繪制的軸心軌跡圖Fig.10 The axis trajectory drawn by the reconstructed displacements

4.2 加速度信號的復指數表達

使用本文方法繪制風機軸心軌跡時，首先要對測試得到的加速度信號進行復指數分解。為了展示更加清晰地分解效果，這里僅選取1 000～1 005 s的信號進行展示。在分解中設置模態階次為600，信號分解后的重構結果與測試得到的加速度信號如圖11所示，可以看到二者具有較好的一致性，這也證明本文使用的基于低階狀態空間模型的復指數分解方法在處理實測信號時，具有較高的分解精度，這也證明分解后的復指數序列可以代表實際的測試信號。

4.3 風機模態軌跡

從圖11中分解后的信號中提取出一階模態后，可以利用式（15）重構x向位移一階模態，同樣也對測試的y向加速度信號進行分解，并利用重構的得到的兩個方向的一階位移模態繪制軸心一階模態軌跡圖，如圖12所示。可以看到相比于圖10中結構的軸心軌跡圖來說，圖12中結構的模態軌跡圖更加清晰，可以直觀的反應出結構的運行狀態，實現對運行風機的監測和檢測。

圖11 分解重構結果Fig.11 The reconstruction result after decomposition

圖12 1 000～1 050 s的結構模態軌跡Fig.12 The axis model trajectory from 1 000-1 050 s

同時，重復上述過程對圖7中2 000～2 050 s的加速度信號進行分析，結構的軸心一階模態軌跡如圖13所示。可以看到本文方法重構的海上風機結構模態軌跡穩定清晰，避免了現場測試中環境噪聲的干擾，可以有效地對結構進行故障特征監測。

圖13 2 000～2 050 s的結構模態軌跡Fig.13 The axis model trajectory from 2 000-2 050 s

5 結 論

本文通過引入低階狀態空間模型，對加速度信號進行高精度的復指數分解，發展了一種海上風機結構的模態軌跡監測新方法。主要結論如下:

（1）該方法使用復指數的形式對信號進行表征，避免了傳統方法中信號周期性的假設，同時通過引入低階狀態空間模型，提高了實測數據的分解精度，避免了傳統Prony方法由于對噪聲敏感而無法處理實測數據的病態問題。

（2）相比于傳統的軸心軌跡分析方法，由于負載的海洋環境，直接用位移得到的結果往往無法實現結構的監測和監測。而本文提出的方法即克服了海上位移測試的難點，也為海上結構結構模態軌跡分析提供了新的思路。

（3）通過五自由度數值算例和實測的海上風電結構加速度信號進行分析，檢驗了本文方法在處理工程實際問題時的有效性。結果表明，本文提出的方法有效地克服了環境噪聲的干擾，得到的海上風電結構模態軌跡穩定清晰，對于結構的健康監測具有較強的工程應用價值。