李 斌，張玉杰，楊智春

(西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)

抖振是現(xiàn)代雙垂尾戰(zhàn)斗機［1］(大攻角機動引起旋渦破裂)和特種電子作戰(zhàn)飛機(外加電線類凸起物引起氣流分離)在服役過程中經(jīng)常要遭遇的復(fù)雜動載荷環(huán)境，屬于典型氣彈響應(yīng)問題。因此對于這類飛機在設(shè)計過程必須精細考慮其抖振載荷對結(jié)構(gòu)安全的影響。在現(xiàn)代飛機結(jié)構(gòu)設(shè)計中，振動嚴重部位的動態(tài)設(shè)計載荷的確定要求合理確定振動過程中的可能最大響應(yīng)、并估計最大響應(yīng)對應(yīng)的最大內(nèi)力載荷。尤其對于具有不確定特性的隨機振動過程，如何根據(jù)有限的測試或計算樣本數(shù)據(jù)，恰當?shù)毓烙嬜畲髣討B(tài)設(shè)計載荷，則顯得更為重要。

飛機設(shè)計的工程實踐中，在隨機振動工況下，確定動態(tài)設(shè)計載荷的傳統(tǒng)習(xí)慣方法是采用所謂的3σ法則，其中σ指振動響應(yīng)的均方根值。該法則的數(shù)學(xué)內(nèi)涵為，Gauss隨機變量在(－3σ，3σ)區(qū)間外取值的概率為0.27%。但實際飛機設(shè)計中對限制載荷的定義是服役中預(yù)期的最大載荷［2］。很顯然3σ法則確定的極值載荷一般是不滿足限制載荷定義的。但是在飛機沒有服役之前，缺乏安全使用壽命飛行數(shù)據(jù)的條件下，難以直接得到最大載荷峰值，只能根據(jù)有限的計算數(shù)據(jù)或風洞試驗數(shù)據(jù)，采用預(yù)估方法進行估計可能出現(xiàn)的載荷最大值。目前在飛機設(shè)計過程，可供參考的一類確定動態(tài)載荷峰值的方法是根據(jù)指定的超越次數(shù)指標來截取極限載荷。如對突風載荷設(shè)計，Hoblit［3］推薦的限制設(shè)計超越次數(shù)是2×10－5次/飛行小時。這種基于設(shè)計超越次數(shù)指標的抖振載荷設(shè)計方法在F－22飛機完成試飛測試，獲得大量抖振試飛數(shù)據(jù)后也得到了應(yīng)用。文獻［4］給出的確定極限載荷的可接受超越指標為全壽命周期內(nèi)107次飛行發(fā)生一次，或超越概率為1×10－4。但這種方法應(yīng)用的前提是必須先建立飛機全壽命周期的抖振載荷超越曲線，這對于設(shè)計階段的抖振載荷預(yù)估確定來說難以滿足。

確定設(shè)計階段的抖振最大設(shè)計載荷，可以利用的數(shù)據(jù)有:縮比模型的風洞試驗數(shù)據(jù)和抖振響應(yīng)計算數(shù)據(jù)。這兩種參考數(shù)據(jù)均只是有限長度，不可能覆蓋全壽命周期。因此，有必要建立一種可靠的統(tǒng)計預(yù)測方法，利用短時間段內(nèi)的響應(yīng)數(shù)據(jù)，來預(yù)測飛機全服役時間內(nèi)可能發(fā)生的最大抖振載荷。極值估計理論是解決這類問題的有效手段之一［5］，目前已經(jīng)在土木結(jié)構(gòu)的風工程研究領(lǐng)域得到較為廣泛的應(yīng)用［6］。本文結(jié)合某型飛機尾部結(jié)構(gòu)的抖振風洞試驗數(shù)據(jù)，運用次序統(tǒng)計和極值估計理論，基于極值I型分布和III型分布，分別應(yīng)用線性的最小二乘估計和最大似然估計進行極值模型的參數(shù)估計，并根據(jù)預(yù)設(shè)極值重現(xiàn)期指標進行抖振響應(yīng)數(shù)據(jù)的極值估計。此外，通過不同分析方法的對比，闡明了傳統(tǒng)的3σ法則的局限性，并給出了應(yīng)用極值估計理論確定動態(tài)設(shè)計載荷時的若干建議。

1 極值估計理論

用X1，X2，…，Xn代表n個獨立的樣本，x是任意樣本Xi的最大絕對值，在大樣本數(shù)目的條件下，Gumbel給出了x的三種漸近分布。這三種分布可以統(tǒng)一用廣義極值分布函數(shù)表示:

式中，β=0稱為極值Ⅰ型分布，實際是一種雙參數(shù)分布函數(shù);β＜0表示極值Ⅱ型分布;β＞0表示極值Ⅲ型分布。后兩種都屬于三參數(shù)分布函數(shù)，分布函數(shù)的尾部較長。參數(shù)α，β，v具有確定的物理含義。β稱為形狀參數(shù)，確定極值分布函數(shù)的形狀;v為位置參數(shù)，表示最頻繁發(fā)生的極值水平，對應(yīng)分布密度函數(shù)的中心峰值位置;α為尺度參數(shù)，表示極值水平隨對數(shù)時間變量的增長率;α/β則表示可能發(fā)生的極值。實際應(yīng)用過程中，極值Ⅱ型或Ⅲ型分布的適用性由估計形狀參數(shù)β決定。

在極值估計中，重現(xiàn)期是一個非常重要的概念，假定單個樣本數(shù)據(jù)占據(jù)的物理時間長度為Δt，則重現(xiàn)期的數(shù)學(xué)表達式為:

式中，t稱為對數(shù)減縮變量，t= －lnln［1/F(x)］。重現(xiàn)期tN的物理意義是觀察到等于或大于x的極值所需的操作時間長度。1－F(x)表示超越概率。

將式(4)代入式(1)可得I型分布函數(shù)的表達式為:

比較式(1)和式(5)，可以發(fā)現(xiàn)，Ⅰ型分布屬于雙參數(shù)分布，估計方便，計算簡單，而且對取值無上下限要求。Ⅲ型分布屬于三參數(shù)分布，可估計的數(shù)據(jù)類型較廣，能夠有效地彌補Ⅰ型分布的缺陷。

2 參數(shù)估計方法

2.1 基于經(jīng)驗分布假設(shè)的最小二乘估計法

對于Ⅰ型分布，x的總體分布函數(shù)如式(5)所示。引入變量:

則式(5)可化為:

將式(7)代入式(3)，推導(dǎo)出x－T關(guān)系式:

或x－t關(guān)系式:

對于Ⅲ型分布，x的總體分布函數(shù)如式(1)所示。引入變量:

則式(1)可化為:

將式(11)代入式(3)推導(dǎo)出x－T關(guān)系式:

或x－t關(guān)系式:

應(yīng)用最小二乘法進行參數(shù)估計的一般步驟為:(1)將選定時間段內(nèi)的數(shù)據(jù)分為n份，從每份中選取最大絕對值x;

(2)將這n個x值按從小到大順序排列;

(3)引入經(jīng)驗分布函數(shù)，經(jīng)過次序排列后的極值分布函數(shù)可以表示為:F(x)=i/n，或F(x)=i/(n+1)，i=1，2，…n。為避免F(x)=1時t無解的情況出現(xiàn)，一般取后者。引入t的定義，則t=－lnln［(n+1)/i］;

(4)對于Ⅰ型分布，考慮到式(9)，引入x關(guān)于t的線性表達式:x=at+b;對于Ⅲ型分布，考慮到式(13)，引入x關(guān)于t的指數(shù)表達式:x=DeEt+F;

(5)依據(jù)經(jīng)驗分布函數(shù)，對應(yīng)每一個級別x，可以計算得到對應(yīng)的t，將它們代入到步驟(4)的參數(shù)方程，可以建立a，b或D，E，F(xiàn)為未知量的矛盾方程組，然后求解該矛盾方程組的最小二乘解。其中求解D，E，F(xiàn)的具體做法是，令d=eEt，根據(jù)d和x相關(guān)系數(shù)最大的原則確定E。當E確定后，就可由線性最小二乘法擬合D和F，詳細算法參考文獻［7］;

(6)將 a，b或 D，E，F(xiàn)的估計值分別回代到式(6)、式(10)中，可得參數(shù) α，β，v的估計值。

2.2 最大似然估計法

下面用最大似然估計法分別導(dǎo)出Ⅰ型分布和Ⅲ型分布參數(shù)的計算公式。

對于如式(7)所示的Ⅰ型分布函數(shù)，分布密度為:

似然函數(shù)為:

對式(15)兩邊取自然對數(shù)可得:

lnL(m，w)取極值的條件為:

可以利用Gauss-Newton迭代法求解式(17)、(18)組合而成的以m，w為未知數(shù)的非線性方程組。本文具體計算時應(yīng)用MATLAB自帶的非線性求解函數(shù)“l(fā)sqnonlin”進行求解。

對于如式(11)所示的Ⅲ型分布函數(shù)，它的分布密度為:

似然函數(shù)為:

對式(20)兩邊取自然對數(shù)可得:

對 lnL(A，B，C)取極大值，得:

同樣可以利用Gauss-Newton迭代法求解式(22)、式(23)、式(24)組合而成的以A，B，C為未知數(shù)的非線性方程組。對于Ⅲ型分布，初值的選取顯著影響計算結(jié)果。參照文獻［7］，A，B，C 初值分別取為 1.5xn，5，1.5xn－0.4，其中，xn為數(shù)量為 n 的極值樣本中的最大值。

3 應(yīng)用實例與分析

應(yīng)用長期極值估計理論，可以建立振動響應(yīng)極值大小與樣本數(shù)量(即重現(xiàn)期)之間的關(guān)系。實際設(shè)計中則可以通過給定重現(xiàn)期來求取對應(yīng)的設(shè)計極值。例如在土木工程結(jié)構(gòu)的設(shè)計規(guī)范中，確定結(jié)構(gòu)風載荷時通常定義“幾十年”一遇的重現(xiàn)期。

對于飛機抖振設(shè)計載荷確定來說，重現(xiàn)期的確定可以遵循如圖1所示的流程框圖來預(yù)設(shè)。首先影響確定影響抖振響應(yīng)烈度的敏感參數(shù)，并用敏感參數(shù)來進行典型狀態(tài)的劃分。然后按照飛機使用習(xí)慣，確定飛機的典型任務(wù)剖面。要求每個典型任務(wù)剖面內(nèi)，每飛行小時各典型狀態(tài)的可能累積發(fā)生時間相對穩(wěn)定。最后根據(jù)飛機的全壽命周期使用統(tǒng)計情況，將整個飛行壽命周期內(nèi)的典型狀態(tài)可能發(fā)生的時間累加起來則可以得到對應(yīng)的“典型狀態(tài)”的設(shè)計重現(xiàn)期。依據(jù)該重現(xiàn)期，則可進行設(shè)計載荷極值的估計。

圖1 抖振設(shè)計載荷確定框圖Fig.1 Flow chart to determine buffet design load

下面以某型飛機尾部結(jié)構(gòu)的抖振響應(yīng)風洞試驗測試數(shù)據(jù)為例進行極值估計方法的實例討論。如圖2所示為通過風洞試驗獲得的某型飛機尾部結(jié)構(gòu)的典型加速度振動響應(yīng)數(shù)據(jù)，數(shù)據(jù)采樣頻率為2 048 Hz，采樣時間24 s。為了有利于統(tǒng)計模型的參數(shù)估計，對原始數(shù)據(jù)幅值進行了無量綱規(guī)范化處理，即將測試數(shù)據(jù)統(tǒng)一處理到±0.8的變化范圍內(nèi)，Ca=對該數(shù)據(jù)進行正態(tài)分布檢驗，統(tǒng)計和擬合得到的概率密度曲線和分布曲線分別如圖3所示，可見響應(yīng)數(shù)據(jù)符合正態(tài)分布假設(shè)，響應(yīng)屬平穩(wěn)隨機過程。

圖2 測試的抖振響應(yīng)樣本Fig.2 Sample of measured buffet response data

圖3 抖振數(shù)據(jù)正態(tài)分布檢驗Fig.3 Probability distribution check of buffet data

3.1 關(guān)于極值模型和估計方法的討論

為了檢驗極值估計模型和估計方法的適用性，以圖2所示的24 s長度的數(shù)據(jù)作為數(shù)據(jù)總體，根據(jù)前文對數(shù)據(jù)的規(guī)范化處理方法，可知其理想極值為0.8。以每128個時域測試數(shù)據(jù)作為一個樣本，則每個子樣本對應(yīng)的時間長度Δt=0.062 5 s。提取各個子樣本內(nèi)的極值(絕對值)，進行極值分布參數(shù)估計分析。應(yīng)用第2節(jié)所述的最小二乘法和最大似然法分別進行極值Ⅰ型分布和Ⅲ型分布的參數(shù)估計，結(jié)果如表1和圖4所示。

表1 Ⅰ型和Ⅲ型分布模型對比Tab.1 Comparison betweenⅠ andⅢ distribution model

觀察表1可以發(fā)現(xiàn)，無論采用最小二乘法還是最大似然估計法，Ⅰ型分布估計出來的極值偏差都比實際值大14%以上，過于保守。Ⅲ型分布估計得到的最大抖振載荷更接近實際值，圖4中兩種模型得到的極值估計曲線與試驗數(shù)據(jù)統(tǒng)計曲線的對比也說明了這一點，Ⅰ型分布估計在初始階段與實際統(tǒng)計曲線吻合較好，但隨著對數(shù)減縮時間變量的增長，曲線后段的偏離越來越大。導(dǎo)致這種現(xiàn)象的主要原因是雙參數(shù)模型較難描述長尾分布。如圖5所示為超越概率分布曲線，Ⅲ型分布模型明顯改善了曲線后段的描述精度，擬合后的總體趨勢曲線與實際統(tǒng)計數(shù)據(jù)之間的吻合性非常好。

圖4 Ⅰ型和Ⅲ型分布估計值和試驗值的對比Fig.4 Comparison of measured and predicted extreme value using distributionⅠandⅢ

圖5 超越概率分布曲線對比Fig.5 Comparison of measured and predicted exceedance probability using distributionⅠandⅢ

表2進一步采用抖振風洞試驗中獲得的10組不同的24 s長度的測試數(shù)據(jù)，基于Ⅲ型分布，分別應(yīng)用最小二乘法和最大似然估計法進行了載荷極值估計，估計結(jié)果如表2所示。對比分析各組數(shù)據(jù)表明，應(yīng)用最小二乘法估計出來的最大抖振響應(yīng)值較試驗值略低，屬于偏冒險的情況;應(yīng)用最大似然估計法，估計出的最大抖振響應(yīng)值較試驗值略高，屬于偏保守的設(shè)計。兩種參數(shù)估計算法得到極值估計精度都在±5%以內(nèi)，適于工程應(yīng)用，且整體來說最大似然估計法的精度要略高于最小二乘法。

表2 兩種參數(shù)估計方法的結(jié)果比較Tab.2 Comparison between results obtained by two parameter estimation methods

3.2 關(guān)于測試數(shù)據(jù)樣本量的討論

由于在實際的動態(tài)信號測量過程中，穩(wěn)定狀態(tài)持續(xù)時間和數(shù)據(jù)采集時間總是有限的。尤其對于飛機的飛行測試試驗，飛行參數(shù)及姿態(tài)持續(xù)穩(wěn)定的時間非常短。本節(jié)將從進行長期極值有效估計的角度，討論測試數(shù)據(jù)樣本量對估計精度的影響，以便確定合理的采樣長度。

以圖2所示24 s秒長度的測試數(shù)據(jù)為例，分別截取時間長度為 0.5 s，1 s，1.5 s，2 s，2.5 s，5 s，24 s 的數(shù)據(jù)，按Ⅲ型分布模型應(yīng)用最大似然估計法進行參數(shù)估計，并進行極值預(yù)測。同時對不同長度的樣本數(shù)據(jù)求相應(yīng)的均方根值，并按照傳統(tǒng)的3σ法則估計設(shè)計載荷，以便對比。

計算結(jié)果如表3所示，應(yīng)用2 s長數(shù)據(jù)，32個統(tǒng)計樣本，進行參數(shù)估計，即可較為準確的識別分布模型參數(shù)，對應(yīng)24 s重現(xiàn)期估計得到的極值與實際極值之間的誤差小于5%。并且由識別參數(shù)隨樣本數(shù)量變化趨勢可見，參數(shù)識別結(jié)果隨樣本數(shù)量的增加趨于穩(wěn)定。而簡單由3σ法估計的設(shè)計載荷明顯低于24 s內(nèi)的極值載荷。如果采用4σ，則基本逼近24s內(nèi)的極值載荷。更進一步，由不同時間長度數(shù)據(jù)得到的均方根值變化可知，對于本文所分析的平穩(wěn)過程數(shù)據(jù)，僅0.5 s數(shù)據(jù)即可準確獲得振動響應(yīng)的基本統(tǒng)計特征。這說明平穩(wěn)過程的數(shù)據(jù)樣本長度對均方根值計算的影響不大。但是如要根據(jù)均方根值來估計長期極值，則需要借用合理的判據(jù)來確定σ的放大倍數(shù)，而本文所提的極值估計法恰好可以用來實現(xiàn)合理倍數(shù)的估計。

3.3 長期極值估計

假定圖2所示的數(shù)據(jù)表示飛機在以下飛行條件下的典型抖振響應(yīng)數(shù)據(jù)。3 000 m高度巡航，馬赫數(shù)Ma=0.6，機翼迎角 αw=4°，方向舵偏角 15°(尾部結(jié)構(gòu)抖振響應(yīng)對方向舵偏角比較敏感，所以選用該參數(shù)進行狀態(tài)劃分)。以該數(shù)據(jù)段作為一個劃分后的典型狀態(tài)代表。

設(shè)每個飛行小時內(nèi)，該狀態(tài)的累計發(fā)生時間為300 s，飛機的設(shè)計壽命為6 000飛行小時。據(jù)此可以推斷在飛機的整個壽命周期內(nèi)該狀態(tài)可能發(fā)生的總持續(xù)時間為 1 800 000 s。取 Δt=0.062 5 s，按照式(2)式(3)設(shè)定tN=1 800 000 s，可求得設(shè)計超越概率為:

將式(25)代入到式(1)，并且根據(jù)事先得到的估計參數(shù)，可求得對應(yīng)的設(shè)計極值為1.145。

4 結(jié)論

本文利用次序統(tǒng)計和長期極值估計理論，結(jié)合某型飛機尾部結(jié)構(gòu)的抖振風洞試驗數(shù)據(jù)，應(yīng)用極值Ⅰ型分布和Ⅲ型分布參數(shù)進行了抖振響應(yīng)極值估計研究，相關(guān)結(jié)論如下:

表3 樣本長度對估計精度的影響(Ⅲ型分布，最大似然法估計)Tab.3 Prediction accuracy varies with tN(distribution Ⅲ ，maximum likelihood method)

(1)對于平穩(wěn)隨機過程，Ⅰ型分布模型的極值估計精度明顯沒有用Ⅲ型分布模型估計的精度高。原因在于Ⅰ型分布為雙參數(shù)模型，難以準確描述長尾分布，三參數(shù)的Ⅲ型分布模型明顯改善了超越概率曲線的后段描述精度。

(2)在數(shù)據(jù)樣本足夠的情況下，最小二乘法估計法和最大似然估計法都可以較為準確地實現(xiàn)分布參數(shù)估計，極值估計精度在±5%范圍內(nèi)，適于工程應(yīng)用。但總體趨勢上最大似然估計法的估計精度略高于最小二乘法。

(3)通過數(shù)據(jù)對比分析可知，對平穩(wěn)隨機振動過程，傳統(tǒng)的3σ設(shè)計原則不能滿足飛機設(shè)計中對確定限制載荷的要求，根據(jù)3σ設(shè)計原則確定的限制載荷將偏危險。長期極值估計法更加符合飛機限制載荷的確定要求，值得在工程設(shè)計中進行推廣應(yīng)用。

(4)要實現(xiàn)飛機抖振載荷準確有效的極值估計，對采集樣本數(shù)量有一定的基本要求。該基本要求可以確定最小響應(yīng)采集時間，該指標對于振動響應(yīng)測試方案和數(shù)據(jù)處理方法的制定有一定的指導(dǎo)性意義。

［1］Lee B H K.Vertical tail buffeting of fighter aircraft［J］.Progress in Aerospace Sciences，2000，36:193 －279.

［2］中國民用航空管理局.中國民用航空規(guī)章第25部:運輸類飛機適航標準(CCAR－25－R4)［S］.北京:中國民用航空局，2011.

［3］HoblitF M. GustLoads on Aircraft:Concepts and Applications［M］.USA:AIAA，1988.

［4］Patel S R，Black C L，Anderson W D.F/A－22 Vertical tail buffet strength certification［C］.46th AIAA/ASME/ASCE/AHS/ASC Structures，StructuralDynamics ＆ Materials Conference，Austin，Texas，2005.

［5］Lee B H K.Statistical analysis of wing/fin buffeting response［J］.Progress in Aerospace Sciences，2002，38:305 －345.

［6］段忠東，歐進萍，周道武.極值風速的最優(yōu)概率模型［J］.土木工程學(xué)報，2002，35(5):11 －16.DUAN Zhong-dong，OU Jin-ping，ZHOU Dao-wu.The Optimal Probabilistic Distribution for Extreme Wind Speed［J］.Journal of China Civil Engineering，2002，35(5):11 －16.

［7］趙 偉，楊永增，于衛(wèi)東，等.長期極值統(tǒng)計理論及其在海洋環(huán)境參數(shù)統(tǒng)計分析中的應(yīng)用［J］.海洋科學(xué)進展，2003.24(4):471－476.ZHAO Wei，YANG Yong-zeng，YU Wei-dong，et al.Longterm extreme value statistics theory and its application to the marine environmental parameter statistics and analysis［J］.Advances in Marine Science，2003，24(4):471 －476.