基于改進局部遞歸率的脈動風速非平穩度分析方法

李利孝， 肖儀清， 鄭 斌， 宋麗莉

（1．哈爾濱工業大學深圳研究生院，廣東深圳518055；2．中國氣象局公共氣象服務中心，北京100081）

基于改進局部遞歸率的脈動風速非平穩度分析方法

李利孝1， 肖儀清1， 鄭 斌1， 宋麗莉2

（1．哈爾濱工業大學深圳研究生院，廣東深圳518055；2．中國氣象局公共氣象服務中心，北京100081）

為了克服遞歸趨勢（recurrence trend，RT）指標對不同信號非平穩度估計存在誤判的缺陷，分別采用互信息法和偽臨近法確定了遞歸量化分析的最佳延遲時間和最小嵌入維數，然后在遞歸量化分析基礎上，提出了歸一化局部遞歸率標準差（standard deviation of normalized local recurrence rate，SDNLRR）作為信號非平穩度量化指標．利用該指標，通過遞歸量化分析方法分析了白噪聲信號、正弦信號、調幅信號、線性調頻信號4個基本信號和2個實測臺風場脈動風速信號的非平穩特性，并與傳統的遞歸趨勢指標分析結果進行了對比．研究結果表明：利用SDNLRR指標對6個信號的非平穩度的量化比較準確率達100%，比RT指標的準確率提高了33．33%，消除了RT指標對正弦信號和平穩脈動風速信號的錯誤估計．

非平穩度；遞歸量化分析；局部遞歸率；脈動風速

天氣系統是一種復雜的非線性動力系統，它的發展演化同時受到大氣層結構狀況、近地面摩擦作用以及湍流驅動機制等多種相互作用且連續變化因素的影響，從而導致實測湍流風速信號具有一定程度的非平穩特征和非高斯特性．現代信號處理的主要分析手段（如譜分析等）以及基于隨機振動理論的結構動力響應計算方法均建立在平穩隨機過程假設基礎上，對于平坦開闊場地，常態風作用下結構動力抗風設計基本滿足該假定．然而對于復雜地形條件，或者處于臺風、雷暴等極端天氣過程條件下的動力抗風設計，則很難滿足平穩性假設［1］．

脈動風速非平穩度的強弱決定了其分析方法的選取以及結構風振響應計算結果的精度和可靠性．通常所用的非平穩檢驗方法（例如輪次檢驗、逆序檢驗）從時間序列的均值與方差等低階統計量出發［2］，僅能給出定性的判別結果，不能對脈動風速的非平穩度進行量化分析．

本文從遞歸分析出發，通過定義新的遞歸量化分析指標——歸一化局部遞歸率標準差法（standard deviation of normalized local recurrence rate，SDNLRR）對時間序列的非平穩度進行刻畫．通過對高斯白噪聲信號、正弦信號、調幅信號、線性調頻信號以及兩個實測脈動風速信號進行非平穩度對比分析來驗證SDNLRR指標的有效性．

1 遞歸圖

遞歸現象在自然界中十分普遍，例如潮汐漲落、日夜更替等．在動力系統理論中，遞歸現象是指系統中的某一特定狀態在不同時刻具有相似的特性．盡管遞歸現象很早就被人們所認識，但是直到1987年Eckmann等提出遞歸圖［3］概念之后，人們才真正掌握了對其進行定性分析和定量描述的手段．遞歸分析通常包括3部分：（1）相空間重構；（2）遞歸圖繪制；（3）遞歸量化分析．

1．1 相空間重構

根據Takens的相空間重構理論［4］，系統中的任一分量的演化和發展過程均是該系統中其他分量對其進行作用和影響的結果，因此，系統中與某一分量有關的信息都將在該分量的時間序列中得到體現．通過時間延遲可將單變量時間序列i＝1，2，…，N｝嵌入到由嵌入維數m和延遲時間τ重構的m維相空間中，即：

式中：M為m維重構相空間的向量數，

為了保證重構的相空間能夠重現原系統的特性，必須合理選擇嵌入維數m與延遲時間τ．Takens證明了當嵌入維數m≥2d＋1（d為吸引子關聯維數）時，重構相空間與原系統在拓撲意義上等價，該關系式提供了無限長無噪聲序列下嵌入維數的選取準則，在該種情況下延遲時間可以任意選取．然而，通常時間序列為有限長度且伴有噪聲，因此嵌入維數m和延遲時間τ的選取在相空間重構的過程中十分重要．

對于m和τ的選取主要分為兩種觀點：第一種認為，嵌入維數m和延遲時間τ之間并無直接聯系，可以分別求取，例如，用G-P算法、偽臨近法、改進偽臨近法及求取最佳延時的自相關函數法和互信息法等求取嵌入維數m；另一種觀點認為兩者之間存在相關關系，其中比較具有代表性的是C-C算法．

本文分別選取了應用較廣的互信息法（mutual information function，MIF）［5］和偽臨近法（false nearest neighbours，FNN）［6］確定最佳延遲時間τ和最小嵌入維數m．

在給定系統A情況下可得到系統B信息量的大小，稱為系統A與B之間的互信息，通過計算原始時間序列之間的互信息，可求得在不同延時τ時兩者之間的關聯程度．由于延遲時間τ選取的基本原則是使序列中兩點之間的相互影響最弱，可以通過改變延遲時間τ，選取其首次得到最小τ值時作為最優延遲時間．

偽臨近點是指兩個在m維重構相空間中接近的點在m＋1維空間發生了遠離，否則稱作真臨近點．當嵌入維數小于系統真實維數時容易產生偽臨近點，當嵌入維數m達到或超過系統真實維數時偽臨近點將消失．因此，通過計算不同嵌入維數時偽臨近點所占比例，找出該比例降至比較穩定的低水平時所對應的嵌入維數即為最小嵌入維數．

圖1為分別用互信息法以及偽臨近法在某脈動風速時程相空間重構過程中確定延遲時間以及嵌入維數的示意圖，從圖1中可以得到嵌入維數m＝5，延遲時間τ＝5．

1．2 遞歸圖

用文獻［3］的方法構造遞歸矩陣R，與延遲后的序列

式中：Rij為遞歸矩陣的元素；

Θ（·）為單位階躍函數；

ε為預先給定的鄰域半徑．

圖1 相空間重構參數Fig．1 Parameters for the phase space reconstruction

遞歸圖是遞歸矩陣的一種可視化表達形式，從式（2）可見，矩陣R是一個由0和1組成的M×M階方陣．計算點在整個區域填充，當Rij＝1時，在遞歸圖上用黑點表示，意味著從狀態Xi出發，經過一段時間后，在狀態Xj發生了遞歸；當Rij＝0時，在遞歸圖上用白色表示，意味著兩種狀態遠離．由于任意狀態均與其本身遞歸，因此遞歸圖的主對角線（line of identity，LOI）為一黑實線．鄰域半徑ε的選取依據目前尚無定論，若過小會增大噪聲影響，而過大則會引起虛假遞歸．本文以保證遞歸率在1%左右的方法來選取的鄰域半徑ε值［7］．

2 遞歸量化分析

遞歸圖的主要分析方法包括定性遞歸分析及量化遞歸分析［8］．前者利用總結歸納得到的遞歸圖模式對一些簡單信號的特征作出判別；后者使用量化分析手段，更適用于湍流等復雜信號．遞歸量化分析方法［8］是由Webber和Zbilut提出的，根據遞歸點在遞歸圖上的分布特征，遞歸量化指標主要包括遞歸率γRR、確定性DET、測度熵ENT、遞歸趨勢RT、層狀度LAM和平均捕獲時間TT等．在這些量化指標中，與非平穩特征相關的主要是遞歸率RR和遞歸趨勢RT．

遞歸率RR用來度量遞歸點密度的大小，可以反映相點在m維空間中的聚集程度以及遞歸頻率的大小，其定義為

遞歸趨勢RT用來表征遞歸率隨著延時的變化趨勢，其定義為

式中：K為將遞歸圖副對角線下方的三角形區域按平行副對角線方向等距分割成的區塊數，其取值通常應當滿足N－K＞10；

γRRk為區塊k（k＝1，2，…，K）內的遞歸率，稱之為局部遞歸率．

遞歸趨勢表征局部遞歸率γRRk關于k的線性回歸斜率值．對于一個完全平穩的時間序列，其局部遞歸率應是一個常數，不隨時間間隔發生變化，遞歸趨勢為0；而對于一些時變的非平穩系統，其局部遞歸率將隨著時間間隔發生變化．目前這種方法得到了廣泛的應用［9-11］．

然而，由于遞歸趨勢是通過最小二乘法線性回歸的斜率值，因此，比較適用于具有漂移模式特征的信號，而對于一些局部非平穩性較強但總體趨勢穩定的信號，可能會出現誤判．湍流信號成分復雜，除了漂移模式外可能還同時存在突變模式等，使用遞歸趨勢來衡量湍流信號可能會存在問題．為了克服遞歸趨勢RT的這一局限性，在這里定義歸一化局部遞歸率標準差SDNLRR來衡量信號的非平穩程度，其定義為

歸一化局部遞歸率標準差σSDNLRR的實質是局部遞歸率的變異系數．非平穩度越強，σSDNLRR值越大，完全平穩時σSDNLRR值為0．與ηRT相比，σSDNLRR對于遞歸率的局部突變更加敏感，對一些模式復雜的非平穩信號分析更有效．另外，由于ηRT是遞歸率的線性擬合斜率，而遞歸率的大小又受到鄰域半徑ε的影響，因此，同一信號的鄰域半徑ε不同時，其計算得到的遞歸趨勢也不同，這對同一個時間序列來說極不合理．而σSDNLRR通過局部遞歸率均值對不同區塊的局部遞歸率歸一化，在一定程度上削弱了鄰域半徑選取對非平穩程度刻畫的影響．

食物有保質期，信息資源也會有。過了期的信息資源就如同昨日黃花一文不值。作為工程造價而言，是為工程做鋪墊的一項工作，對于時間肯定有著嚴格的要求。因為工程造價對于工程而言是一個必要且充分的條件，但又不能因為它的拖延而推遲工程的期限[4]。所以，為了造價方案在預計時間段內做出對于信息資源就必須以最快的速度獲取，并且是最新的、具有時效性的信息資源。而大數據下的計算機和各類系統恰如其分地解決了這個問題。

3 數值模擬及實測驗證

3．1 信號來源

為了驗證本文定義的非平穩度量化指標σSDNLRR的有效性，選取信號S1～S6進行驗證．

（1）理論上平穩的高斯白噪聲信號S1由MATLAB的randn函數生成．

（2）正弦信號S2為x（t）＝sin（2π×0．02t）．

（3）線性調頻信號S3為x（t）＝sin（2πf1t），ft＝5×10－5t＋0．02，可以產生線性變化的頻率為0．02～0．08 Hz正弦調頻信號．

（4）調幅信號S4為x（t）＝A［1＋M× cos（2πfmt）］sin（2πfct），式中：載波幅值A＝1；載波基頻fc＝0．02 Hz；調制頻率fm＝0．005 Hz；調制幅度為50%．

（5）同時通過逆序檢驗與輪次檢驗的“平穩”脈動風速時程信號S5選自臺風黑格比過程的實測脈動風速時程．

（6）未通過逆序檢驗和輪次檢驗的非平穩脈動風速時程S6選自了臺風黑格比過程的實測脈動風速時程．

在驗證的過程中，同時計算了傳統的ηRT值作為對比，信號S1～S6的時程曲線如圖2所示，各信號采樣頻率均為1 Hz，樣本長度N均為600．

圖2 信號時程圖Fig．2 Time history of the signals

3．2 遞歸圖及遞歸分析

采用第1節所述方法對上述6個信號進行相空間重構和遞歸圖繪制．為了統一標準，K值均取為10，按照控制總體遞歸率為1%左右的準則選取鄰域半徑ε，嵌入維數m和延遲時間τ分別采用偽臨近法（FNN）和互信息法（MIF）計算．根據上述參數，計算了信號S1～S6的遞歸圖及相應的局部遞歸率，見圖3和圖4．利用上述遞歸圖及相應的局部遞歸率結果，根據式（4）和（5）分別計算了各信號的ηRT和σSDNLRR，見表1．

從圖3、4及表1可見，高斯白噪聲信號S1的遞歸點整體分布均勻，其ηRT值及σSDNLRR值均為所測試的6種信號中最小，因此，兩種方法都能較好地表征白噪聲信號的非平穩度．

正弦信號S2是一種平穩性較好的周期信號，雖然其遞歸圖不如信號S1分布均勻，但整體分布規則，表現出許多與主對角線平行的直線．從S2的局部遞歸率分布圖可見，除第10區塊由于面積過小引起誤差偏大外，其余區塊的局部遞歸率波動較小．但正是由于γRR10的影響使得S2的ηRT＝－0．081 2，高于調幅信號S4，很不合理；其σSDNLRR值則為0．196 0，略大于信號S1，小于調頻信號S3與調幅信號S4，與實際情況相符．

調頻信號S3的遞歸圖呈現了漂移模式，其局部遞歸率具有單調變化的趨勢，因此存在著比較明顯的非平穩性．從圖3可見，遞歸點在主對角線附近集中，γRR1約為2．2，明顯大于平均值；信號S3的

兩者均反映出信號的非平穩特性，相對而言，σSDNLRR更顯著．

圖3 信號遞歸圖Fig．3 Recurrence plots of the signals

圖4 信號局部遞歸率Fig．4 Local recurrence rate of the signals

理論上調幅信號S4較正弦信號S2的平穩性更差，其遞歸圖的也證明了這一點．S4的遞歸圖僅由幾條與主對角線平行的直線構成，整體均勻程度較S2差，非平穩性比S2更顯著．從S4的局部遞歸率圖上可以看到，除γRR4和γRR8外，其他區塊遞歸率為0．然而，由于S4的局部遞歸率變化比較對稱，其ηRT值僅為0．068 3，小于信號S2，未能反映出該信號應有的非平穩特征；相對而言，σSDNLRR＝2．121 7則更加合理．

分別計算信號S5和S6的指標ηRT與σSDNLRR，從圖3（e）和圖3（f）可見，信號S5的遞歸點分布較S6更均勻，因此，S5比S6的平穩性更好，符合實際情況．S5的ηRT＝0．040 6，σSDNLRR＝0．282 7，與平穩信號S1和S2接近，表明信號S5較為平穩．信號S6遞歸圖上的遞歸點集中于左下角和右上角，分布不均勻．對照信號S6的時程可以發現，在約350 s風速發生了突變，因此S6的非平穩性明顯，其ηRT＝－0．168 2，σSDNLRR＝1．428 3．因此，采用指標σSDNLRR對脈動風速時程進行非平穩程度量化描述是合理有效的．

綜上所述可見，傳統的遞歸趨勢TND對上述6個信號中的正弦信號和調幅信號的非平穩度進行量化分析時存在估計不準確的問題，而采用本文定義的歸一化局部遞歸率標準差SDNLRR可以準確地對上述6個信號的非平穩度進行量化分析，適用于實測脈動風速信號的非平穩特性分析．

表1 平穩度量化指標比較Tab．1 Comparison between quantitative non-stability indicatorsηTNDandσSDNLRR

4 結束語

脈動風速的非平穩性度量及其合理設計是結構抗風設計中亟待解決的問題，本文基于遞歸量化分析方法，對傳統的非平穩性度量指標進行了改進，提出了歸一化局部遞歸率標準差σSDNLRR對時間序列的非平穩性進行表征，并使用了4種通用信號和兩個實測脈動風速信號進行了比較驗證．結果表明，與傳統的非平穩性指標（遞歸趨勢指標ηRT等）相比，本文提出的歸一化局部遞歸率標準差σSDNLRR對信號非平穩性的表征更合理．

本文方法通過對局部遞歸率進行歸一化，在一定程度上削弱了由不同樣本之間總體遞歸率的差別引起的非平穩性差異．通過計算局部遞歸率的標準差而不是對其斜率進行線性擬合，擴展了非平穩性指標的適用范圍．最后，對實測“平穩”脈動風速和非平穩脈動風速的驗證結果表明，指標σSDNLRR能較好地表征不同實測脈動風速的非平穩性．

［1］ CHEN J，XU Y L．On modelling of typhoon-induced non-stationary wind speed for tall buildings［J］．The Structural Design of Tall and Special Buildings，2004，13（2）：145-163．

［2］ BECK T W，HOUSH T J，WEIR J P，et al．An examination of the runs test，reverse arrangements test，and modified reverse arrangements test for assessing surface EMG signal stationarity［J］．Journal of Neuroscience Methods，2006，156（1）：242-248．

［3］ ECKMANN J，KAMPHORST S O，RUELLE D．Recurrence plots of dynamical systems［J］．Europhysics Letters，1987，4（9）：973-977．

［4］ TAKENS F．Detecting strange attractors in turbulence［M］．［S．l．］：Springer，1981：366-381．

［5］ FRASER A M，SWINNEY H L．Independent coordinates for strange attractors from mutual information［J］．Physical Review A，1986，33（2）：1134-1140．

［6］ KENNEL M B，BROWN R，ABARBANEL H D．Determining embedding dimension for phase-space reconstruction using a geometrical construction［J］．Physical Review A，1992，45（6）：3403-3411．

［7］ ZBILUT JP，ZALDIVAR-COMENGES J，STROZZIF．Recurrence quantification based Liapunov exponents for monitoring divergence in experimental data［J］．Physics Letters A，2002，297（3）：173-181．

［8］ WEBBER C L，ZBILUT JP．Dynamical assessment of physiological systems and states using recurrence plot strategies［J］．Journal of Applied Physiology，1994，76（2）：965-973．

［9］ ZBILUT J P，THOMASSON N，WEBBER C L．Recurrence quantification analysis as a tool for nonlinear exploration of nonstationary cardiac signals［J］．Medical Engineering and Physics，2002，24（1）：53-60．

［10］ 閆潤強．語音信號動力學特性遞歸分析［D］．上海：上海交通大學，2006．

［11］ 楊棟，任偉新，李丹，等．基于局部遞歸率分析的振動信號非平穩評價［J］．中南大學學報：自然科學版，2013，44（7）：3024-3032．YANG Dong，REN Weixin，LI Dan，et al．Local recurrence rate analysis based non-stationarity measurement for operational vibration signal［J］．Journal of Central South University：Science and Technology，2013，44（7）：3024-3032．

（中文編輯：秦萍玲 英文編輯：蘭俊思）

Method for Analysis of Non-stationarity of Fluctuating Winds Based on Revised Local Recurrence Rate

LI Lixiao1， XIAO Yiqing1， ZHENG Bin1， SONG Lili2
（1．Shenzhen Graduate School，Harbin Institute of Technology，Shenzhen 518055，China；2．Public Meteorological Service Center，China Meteorological Administrator，Beijing 100081，China）

In order to overcome the misestimate of the non-stationarity of different signals by the recurrence trend（RT），the mutual information function and false nearest neighbors are employed to determine the time delay and minimum embedding dimension of the recurrence quantification analysis（RQA），respectively．Then a novel index，i．e．，the standard deviation of normalized local recurrence rate（SDNLRR），which is based on the RQA，was proposed to quantify the non-stationary degree of the signals．Utilizing the SDNLRR，the non-stationarity of four basic signals（white noise signal，sinusoidal signal，amplitude-modulated signal and linear frequency modulation signal）and two field observed fluctuating wind speed histories were analyzed and compared with the analysis results of RT．The results show that the proposed SDNLRR could offer a quantitative comparison of the nonstationarity of the above six signals with a 100%accuracy．The new method eliminates the misestimates of the sinusoidal signal and the stationary fluctuating wind speed signal in RT，and hence is more accurate than the RT estimation by 33．33%．

non-stationarity；recurrence quantification analysis；local recurrence rate；fluctuating winds

TU973．213

A

0258-2724（2016）01-0065-06

10．3969／j．issn．0258-2724．2016．01．010

2014-07-30

國家自然科學基金資助項目（51308168，51278161）；中國博士后科學基金資助項目（2013M531045，2014T70343）

李利孝（1984—），男，博士，研究方向為結構風工程，電話：18682013431，E-mail：lilixiao1984＠gmail．com

李利孝，肖儀清，鄭斌，等．基于改進局部遞歸率的脈動風速非平穩度分析方法［J］．西南交通大學學報，2016，51（1）：65-70．