體自由度顫振飛行試驗(yàn)的邊界預(yù)測(cè)方法研究

雷鵬軒 呂彬彬 郭洪濤 余立 陳德華

摘要：體自由度顫振頻率低，參與顫振的模態(tài)頻率在亞臨界狀態(tài)往往已經(jīng)極為接近，加之基于大氣紊流激勵(lì)的顫振飛行試驗(yàn)數(shù)據(jù)信噪比通常較低，增加了體自由度顫振飛行試驗(yàn)?zāi)B(tài)辨識(shí)以及顫振預(yù)測(cè)的難度。對(duì)此提出了一種基于Matrix Pencil模態(tài)辨識(shí)方法的體自由度顫振預(yù)測(cè)方法。通過隨機(jī)減量技術(shù)對(duì)輸出響應(yīng)信號(hào)進(jìn)行系集平均，得到隨機(jī)衰減標(biāo)記;運(yùn)用Matrix Pencil方法擬合隨機(jī)衰減標(biāo)記獲取模態(tài)參數(shù)，并通過頻率和阻尼穩(wěn)定判據(jù)篩選真實(shí)模態(tài)，再通過阻尼比與顫振穩(wěn)定性判據(jù)變量外插獲取顫振點(diǎn)。通過對(duì)仿真數(shù)據(jù)與試驗(yàn)數(shù)據(jù)的應(yīng)用，可得到以下結(jié)論：Matrix Pencil模態(tài)辨識(shí)方法能有效辨識(shí)密集的顫振模態(tài)，并獲得清晰的模態(tài)辨識(shí)穩(wěn)態(tài)圖?；谧枘岜燃胺€(wěn)定性判據(jù)變量外插獲得的顫振預(yù)測(cè)結(jié)果較為合理，其中DTFM（ Discrete-Time Flutter Margin）判據(jù)變量的下降趨勢(shì)更明顯，外插結(jié)果與試驗(yàn)值更接近。該方法適用于體自由度顫振飛行試驗(yàn)的亞臨界預(yù)測(cè)。

關(guān)鍵詞：體自由度顫振;模態(tài)辨識(shí);顫振預(yù)測(cè);飛行試驗(yàn);Matrix Pencil方法

中圖分類號(hào)：V211.47

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：10044523（ 2022）01-020207

DOI： 10.16385/j .cnki.issn.10044523.2022.01.022

引 言

飛翼飛行器剛體短周期模態(tài)頻率高，易與一階彈性彎曲模態(tài)耦合發(fā)生一種特殊的顫振——體自由度顫振[1]。體自由度顫振的試驗(yàn)研究需要準(zhǔn)確模擬剛體自由度。受洞壁尺寸限制，在風(fēng)洞中開展動(dòng)力相似模型的自由飛試驗(yàn)是困難的。若僅模擬若干關(guān)鍵剛體自由度，也將不可避免地引入支撐結(jié)構(gòu)阻尼的不利影響，因此飛行試驗(yàn)成為其最有效的研究手段[2-3]。目前國(guó)外已將飛行試驗(yàn)列為體自由度顫振研究的主要手段，如2010年AFRL聯(lián)合洛·馬公司提出了建立X-56A飛行演示平臺(tái)（Multi-UtilityTechnology Testbed，MUTT）計(jì)劃[4-5]，其目的之一就是研究體自由度顫振。2015年美國(guó)明尼蘇達(dá)大學(xué)開展了飛翼外形（ mAEWingl，mAEWing2）體自由度顫振飛行試驗(yàn)研究[6-10]。

在顫振試驗(yàn)（飛行、風(fēng)洞）中，紊流激勵(lì)是常用的激勵(lì)方式之一。該方式不需要外部激勵(lì)部件，易于實(shí)現(xiàn)且更貼近于飛行器真實(shí)狀態(tài)。Scott等[11]、Gu等[12].劉基海等[13]以及黃超等[14-15]都在體自由度顫振風(fēng)洞試驗(yàn)中采用了紊流激勵(lì)的試驗(yàn)方案，其中Huang將ARMA模態(tài)辨識(shí)方法應(yīng)用于體自由度顫振預(yù)測(cè)當(dāng)中。但由于紊流激勵(lì)不充分，僅使用了來流速度高于30 m/s的幾個(gè)階梯，同時(shí)辨識(shí)結(jié)果散度也較大，且缺乏試驗(yàn)直接獲得的顫振點(diǎn)進(jìn)行驗(yàn)證。紊流激勵(lì)雖然易于實(shí)施，但試驗(yàn)數(shù)據(jù)信噪比往往較低，對(duì)于飛行試驗(yàn)還存在一定的非平穩(wěn)特性，且數(shù)據(jù)量有限。這些都對(duì)于基于飛行試驗(yàn)數(shù)據(jù)的顫振亞臨界預(yù)測(cè)提出了較高要求[16]。此外，現(xiàn)有的體自由度顫振試驗(yàn)數(shù)據(jù)反映出，即使在距離顫振邊界較遠(yuǎn)的亞臨界狀態(tài)，參與顫振的模態(tài)頻率也極為接近。這進(jìn)一步加大了模態(tài)辨識(shí)與顫振預(yù)測(cè)的難度。

針對(duì)體自由度顫振飛行試驗(yàn)中數(shù)據(jù)信噪比低、模態(tài)密集的特點(diǎn)，本文提出了一種適用于紊流激勵(lì)的體自由度顫振亞臨界預(yù)測(cè)方法，并基于已開展的體自由度顫振飛行試驗(yàn)數(shù)據(jù)，對(duì)顫振預(yù)測(cè)方法進(jìn)行了驗(yàn)證。

1 亞臨界預(yù)測(cè)方法

本文提出了一種適用于紊流激勵(lì)的顫振亞臨界預(yù)測(cè)方法，并將其應(yīng)用于飛行試驗(yàn)下的體自由度顫振預(yù)測(cè)。首先，通過隨機(jī)減量分析對(duì)輸出響應(yīng)信號(hào)進(jìn)行系集平均，得到隨機(jī)減量標(biāo)記;其次，運(yùn)用Ma-trix Pencil方法擬合隨機(jī)減量標(biāo)記獲取模態(tài)參數(shù)[17];最后，通過阻尼比以及穩(wěn)定性判據(jù)變量外插顫振臨界點(diǎn)。

1.1 隨機(jī)減量分析

隨機(jī)減量法是以一定的初始電平Uo（初始位移）為條件，采集經(jīng)過濾波的一組模型響應(yīng)子樣數(shù)據(jù)Ul（t），U2（t），…，UN（t）（如圖1（a）所示），經(jīng)過集系平均，消除由風(fēng)洞噪聲、電磁干擾等產(chǎn)生的隨機(jī)響應(yīng)R（t）（圖1（c）），保留因?yàn)槌跏紬l件產(chǎn)生的有效物理振動(dòng)響應(yīng)S（t）（如圖1（b））。平均后由噪聲引起的強(qiáng)迫隨機(jī)響應(yīng)的均方根為：

當(dāng)規(guī)定初始位移作為子樣截取的初始條件時(shí)，經(jīng)集系平均后，因?yàn)檎?fù)初始速度出現(xiàn)的概率相等，所以由初始速度引起的脈沖響應(yīng)相互抵消。因此S（t）中實(shí)際只包含由初始位移引起的氣動(dòng)力激勵(lì)響應(yīng)。可見，隨著平均數(shù)N的加大，R（t）將趨近于零，系集平均結(jié)果U（t）（稱為隨機(jī)衰減標(biāo)記）趨近于瞬時(shí)響應(yīng)S（t）。瞬時(shí)響應(yīng)可能是單型態(tài)的;也可能有兩個(gè)頻率接近的型態(tài)（如圖2所示）。

1.2

Matrix Pencil方法

對(duì)于一個(gè)具有M自由度的彈性系統(tǒng)來說，其振動(dòng)運(yùn)動(dòng)微分程為：

1.3 穩(wěn)定性判據(jù)

能夠用于指示顫振臨界點(diǎn)的穩(wěn)定性判據(jù)有很多，比如功率譜峰值倒數(shù)、阻尼比、Routh判據(jù)和DTFM判據(jù)等。其中阻尼比的物理概念清晰，應(yīng)用最為普遍，由于顫振發(fā)生時(shí)，至少有某一階模態(tài)的阻尼變?yōu)榱?，從而造成自激振?dòng)引發(fā)顫振，因此可以以阻尼作為顫振是否發(fā)生的判據(jù)，但是阻尼外推法往往存在著辨識(shí)誤差較大，且在亞臨界階段下降趨勢(shì)可能不明顯的問題。顫振基于彎扭二自由度運(yùn)動(dòng)方程，利用Routh判據(jù)確定系統(tǒng)的穩(wěn)定條件。通常來說對(duì)于頻率的辨識(shí)精度要高于阻尼，通過將發(fā)生耦合的模態(tài)的頻率信息引入預(yù)測(cè)判據(jù)，可以顯著提高預(yù)測(cè)精度，因此顫振裕度法可以提高顫振預(yù)測(cè)精度，但事先要知道發(fā)生耦合的模態(tài)。顫振裕度法以彎扭二自由度顫振方程為基礎(chǔ)，其理論同樣適用于俯仰模態(tài)與一彎模態(tài)耦合所致的體自由度顫振。對(duì)于二階自由度的彈性系統(tǒng)來說，其特征方程為：

當(dāng)Fx<0時(shí)，系統(tǒng)失穩(wěn)。畫出穩(wěn)定性參數(shù)隨著飛行速度變化曲線即可擬合外插出顫振臨界點(diǎn)，大量高速風(fēng)洞顫振試驗(yàn)表明，DTFM方法預(yù)測(cè)精度高、魯棒性強(qiáng)，是目前高速風(fēng)洞顫振試驗(yàn)常用的主流預(yù)測(cè)方法。gzslib202204041820

2 仿真信號(hào)應(yīng)用

本節(jié)以包含兩個(gè)密集模態(tài)的脈沖響應(yīng)信號(hào)為例檢驗(yàn)Matrix Pencil方法對(duì)密集模態(tài)的辨識(shí)能力。考慮到體自由度顫振頻率往往較低，系統(tǒng)的兩階模態(tài)參數(shù)分別為：f1=5 Hz，ξ1=0.05; f2=5.5 Hz，ξ2=0.05。同時(shí)在信號(hào)中加入白噪聲信號(hào)，白噪聲最大幅值與有效信號(hào)最大幅值比為1/10。帶噪聲的脈沖響應(yīng)信號(hào)如圖3所示，信號(hào)采樣頻率為100 Hz，長(zhǎng)度為400。

圖4和表1給出了Matrix Pencil方法的辨識(shí)結(jié)果，可見采用該方法可辨識(shí)獲得規(guī)律清晰的穩(wěn)態(tài)圖。采用6～20階模型進(jìn)行辨識(shí)，均獲得了穩(wěn)定的兩個(gè)極點(diǎn)，且兩個(gè)極點(diǎn)的模態(tài)頻率與真實(shí)值極為接近，模態(tài)阻尼辨識(shí)精度略差，但也基本滿足顫振試驗(yàn)中對(duì)密集模態(tài)辨識(shí)的需求。

3 飛行試驗(yàn)的應(yīng)用

本節(jié)采用實(shí)測(cè)的體自由度顫振飛行試驗(yàn)數(shù)據(jù)對(duì)上述顫振預(yù)測(cè)方法進(jìn)行檢驗(yàn)。首先對(duì)實(shí)驗(yàn)?zāi)Ｐ湍B(tài)特性做簡(jiǎn)要說明：在飛行試驗(yàn)開始前，在停機(jī)狀態(tài)下開展了地面振動(dòng)試驗(yàn)。表2給出了模態(tài)頻率實(shí)驗(yàn)結(jié)果，可見試驗(yàn)結(jié)果與設(shè)計(jì)目標(biāo)吻合良好。

飛行試驗(yàn)總共有3個(gè)亞臨界速度階梯，分別為28，30，26 m/s（如圖5所示）。此外在階梯1與階梯2之間以及階梯3之后，都通過飛行試驗(yàn)復(fù)現(xiàn)了體自由度顫振現(xiàn)象，直接獲得了顫振點(diǎn)。兩次出現(xiàn)顫振時(shí)的實(shí)時(shí)飛行速度分別為33.91和33.77 m/s，顫振頻率均為3.2 Hz。圖6給出了顫振發(fā)生時(shí)連續(xù)的6幀畫面，可以看出一個(gè)完整的振動(dòng)周期，對(duì)比遠(yuǎn)處的地平線變化可發(fā)現(xiàn)，機(jī)翼上反時(shí)機(jī)身抬頭;機(jī)翼下反時(shí)機(jī)身低頭，剛體運(yùn)動(dòng)與彈性振動(dòng)接近同相位。兩次顫振直飛試驗(yàn)結(jié)果非常接近。該結(jié)果可作為試驗(yàn)直接獲取的體自由度顫振速度對(duì)顫振預(yù)測(cè)結(jié)果進(jìn)行檢驗(yàn)。

3.1 模態(tài)辨識(shí)結(jié)果

基于三個(gè)亞臨界階梯的振動(dòng)信號(hào)開展顫振模態(tài)辨識(shí)。每個(gè)階梯選取信號(hào)段長(zhǎng)20 s，數(shù)據(jù)長(zhǎng)度2000。首先采用隨機(jī)減量法對(duì)數(shù)據(jù)段進(jìn)行降噪處理，提取隨機(jī)衰減標(biāo)記。信號(hào)處理結(jié)果如圖7所示。

針對(duì)所提取的隨機(jī)衰減標(biāo)記，采用Matrix Pen-cil方法進(jìn)行模態(tài)參數(shù)辨識(shí)。為了有效區(qū)分虛假模態(tài)，采用穩(wěn)態(tài)圖方法確定模型階次，應(yīng)用阻尼穩(wěn)定和頻率穩(wěn)定判據(jù)對(duì)辨識(shí)結(jié)果進(jìn)行了篩選。

首先采用較為簡(jiǎn)單的頻率穩(wěn)定、阻尼穩(wěn)定以及模態(tài)貢獻(xiàn)度判據(jù)對(duì)辨識(shí)結(jié)果進(jìn)行了初步篩選，初級(jí)判據(jù)如下式所示：式中模態(tài)貢獻(xiàn)度Ai如下式所示，定義為模態(tài)幅值ai占整體全部模態(tài)幅值之和的百分比：

根據(jù)體自由度顫振頻率，及顫振問題所關(guān)心的阻尼比范圍，設(shè)定頻率篩選范圍為0～10 Hz，阻尼比篩選范圍為0～0.3，模態(tài)貢獻(xiàn)度大于5%。

經(jīng)過初步篩選后，依據(jù)連續(xù)兩個(gè)模型階次所辨識(shí)的頻率與阻尼相對(duì)誤差作進(jìn)一步穩(wěn)定模態(tài)的篩選，高級(jí)判據(jù)如下式所示：

一般來說，頻率的辨識(shí)精度較高，阻尼的辨識(shí)精度較差，因此選取εf和εξ分別為5%和10%。

圖8給出了采用6～20階模型辨識(shí)獲得的穩(wěn)態(tài)圖。圖中黑色點(diǎn)表示滿足初級(jí)判穩(wěn)判據(jù)，但不滿足高級(jí)判據(jù)的極點(diǎn)。紅色加號(hào)代表滿足頻率穩(wěn)定的極點(diǎn)，綠色圓圈代表阻尼穩(wěn)定的極點(diǎn)。

首先從圖8所示的亞臨界狀態(tài)功率譜曲線可見，對(duì)于體自由度顫振，即使在亞臨界狀態(tài)參與顫振的模態(tài)間頻率也非常接近，從頻譜上幾乎難以觀察出明顯的多個(gè)頻率峰值，尤其是飛行速度為30 m/s時(shí)的功率譜，幾乎呈現(xiàn)出單頻振動(dòng)。若無法準(zhǔn)確辨識(shí)出參與顫振的模態(tài)參數(shù)，將會(huì)直接影響顫振預(yù)測(cè)結(jié)果。

其次，從辨識(shí)出的極點(diǎn)穩(wěn)定性可見，采用MatrixPencil方法可以準(zhǔn)確辨識(shí)出多個(gè)穩(wěn)定模態(tài)。從貢獻(xiàn)度分布圖中可見辨識(shí)出的模態(tài)分支總貢獻(xiàn)度占到了整段振動(dòng)信號(hào)的80%以上。其中從飛行速度26 m/s時(shí)的功率譜中辨識(shí)出三個(gè)穩(wěn)定極點(diǎn)，參考模態(tài)貢獻(xiàn)度的大小，選取較大的兩階模態(tài)進(jìn)行顫振預(yù)測(cè)。

此外，圖8還給出了模態(tài)阻尼比的聚類分布，并圈出了聚類分布的分散程度△ξ，以阻尼比的絕對(duì)誤差表示?？梢钥闯?，各模態(tài)分支的阻尼比分布較為緊湊。隨著飛行速度增大，紊流擾動(dòng)加強(qiáng)，貢獻(xiàn)度較高的模態(tài)分支阻尼比辨識(shí)精度顯著增強(qiáng)，從±2.6%提高至±0.4%。

最終從整體上看，基于Matrix Pencil方法的模態(tài)辨識(shí)結(jié)果合理可靠，能實(shí)現(xiàn)密集模態(tài)的有效辨識(shí)，獲取真實(shí)的亞臨界顫振模態(tài)。表3給出了三個(gè)階梯的模態(tài)辨識(shí)結(jié)果。

3.2 顫振預(yù)測(cè)結(jié)果

基于上一節(jié)的模態(tài)辨識(shí)結(jié)果，首先采用阻尼比作為判據(jù)進(jìn)行顫振預(yù)測(cè)。圖9給出了模態(tài)阻尼隨飛行速度的變化曲線。從圖中可以看出，最終的發(fā)散模態(tài)為模態(tài)1（俯仰模態(tài)）。由于模態(tài)阻尼隨飛行速度的下降趨勢(shì)并不明顯，采用了二次多項(xiàng)式擬合和基于最后兩個(gè)階梯的線性擬合進(jìn)行外插，預(yù)測(cè)的顫振速度分別為30.43和30.67 m/s。顫振頻率采用發(fā)散模態(tài)分支的頻率變化曲線外插獲取，預(yù)測(cè)的顫振頻率為2.82 Hz。

由于阻尼比無法反映出顫振發(fā)生前模態(tài)頻率耦合的過程，采用既考慮頻率影響又考慮阻尼影響的穩(wěn)定性判據(jù)進(jìn)行顫振預(yù)測(cè)。從圖10和11可見，穩(wěn)定性判據(jù)隨飛行速度的增大下降趨勢(shì)更明顯，且基本呈線性關(guān)系。兩種方法的顫振外插結(jié)果分別為30.90和33.42 m/s，較阻尼比外插結(jié)果都更貼近于真實(shí)值。表4給出了三種方法的顫振預(yù)測(cè)結(jié)果，可見顫振預(yù)測(cè)誤差均在10%以內(nèi)，且均為保守。其中DT-FM方法預(yù)測(cè)結(jié)果最準(zhǔn)確，誤差僅為 1.12%。

4 結(jié) 論

針對(duì)體自由度顫振飛行試驗(yàn)信噪比低、顫振模態(tài)密集的特點(diǎn)，本文提出了一種基于Matrix Pencil模態(tài)辨識(shí)方法的適用于紊流激勵(lì)的顫振預(yù)測(cè)方法。經(jīng)密集模態(tài)測(cè)試信號(hào)和飛行試驗(yàn)數(shù)據(jù)檢驗(yàn)，所采用的模態(tài)辨識(shí)方法能實(shí)現(xiàn)密集模態(tài)的有效辨識(shí)，獲得模態(tài)關(guān)系清晰的穩(wěn)態(tài)圖。通過阻尼比和穩(wěn)定性判據(jù)變量外插能獲得較為合理的體自由度顫振預(yù)測(cè)結(jié)果。其中，DTFM方法構(gòu)造的穩(wěn)定性判據(jù)變量Fz下降趨勢(shì)明顯，且預(yù)測(cè)結(jié)果與飛行試驗(yàn)直接獲得的顫振點(diǎn)最為接近，表明本文提出的顫振預(yù)測(cè)方法有效可靠，能實(shí)現(xiàn)體自由度顫振飛行試驗(yàn)的數(shù)據(jù)處理和顫振預(yù)測(cè)。gzslib202204041820

參考文獻(xiàn)：

[1]楊超，黃超，吳志剛，等.氣動(dòng)伺服彈性研究的進(jìn)展與挑戰(zhàn)[J].航空學(xué)報(bào)，2015，36（4）：1011-1033.

Yang Chao， Huang Chao， Wu Zhigang， et al.Progressand challenges for aeroservoelasticity research[J].ActaAeronautica et Astronautica Sinica， 2015， 36（4）：10111033.

[2] Love M H， Zink P S，Wieselmann P A， et al.Bodyfreedom flutter of high aspect ratio flying wings [C].Collection of Technical PapersAIAA/ASME/ASCE/AHS/ASC Structures ， Structural Dynamics and Mate-rials Conference ， 2005 ， 3 ： 1808-1830.

[3]Nicolai L， Hunten K， Zink P S， et al. System benefitsof active flutter suppression for a sensorcrafttype vehicle[C]. 13th AIAA/ISSMO Multidisciplinary AnalysisOptimization Conference ， 2010 ： 9349.

[4]Bumett E， Atkinson C， Beranek J， et al. Ndof simula-tion model for flight control development with flight testcorrelation [C]. AIAA Modeling and Simulation Tech-nologies Conference， 2010： 7780.

[5]Beranek J， Nicolai L， Buonanno M. et al. Conceptualdesign of a multiutility aeroelastic demonstrator [C].13th AIAA/ISSMO Multidisciplinary Analysis and Op-timization Conference ， 2010 ： 9350.

[6]Regan C D， Taylor B R. mAEWingl： Design， build，testinvited [C] . AIAA Atmospheric Flight MechanicsConference ， San Diego ， USA ， 2016 ： 1747.

[7]Gupta A， Seiler P J， Danowsky B P. Ground vibrationtests on a flexible flying wing aircraft [ C ] . AIAA Atmo-spheric Flight Mechanics Conference ， 2016 ： 1753.

[8]Kotikalpudi A. Pfifer H. Balas G J. Unsteady aerody-namics modeling for a flexible unmanned air vehicle[C]. AIAA Atmospheric Flight Mechanics Confer-ence . 2015 ： 2854.

[9]Schmidt D K， Zhao W， Kapania R K. Flightdynamicsand flutter modeling and analyses of a flexible flyingwing drone [C] . AIAA Atmospheric Flight MechanicsConference ， 2016 ： 1748.

[10]Danowsky B P， Kotikalpudi A， Schmidt D K， et al.Flight testing flutter suppression on a small flexible fly-ingwing aircraft [J]. 2018 Multidisciplinary Analysisand Optimization Conference， DOI： 10.2514/6.20183427.

[11]ScottRC，VetterT K，PenningKB，et al.AeroServ0～elastic testing of a sidewall mounted free nying windtunnel model[C]. Collection of Technical PapersAIAA Applied Aerodynamics Conference， 2008，DOI：10.2514，6.2008 7186.

[12]Gu Yingsong，Yang Z C，Marzocca P，et al.Body free—dom nutter of a nexible blended wing body like plateAn experimental study[C].17th Intemational Forumon Aeroelasticity and Structural Dynamics， IFASD2017，Como，Italy，2017，4（4）：2558 2573.

[13]劉基海，谷迎松，楊智春.支撐剛度對(duì)飛翼模型固有模態(tài)和體自由度顫振特性的影響[J].振動(dòng)工程學(xué)報(bào)，2018，31（5）：727 733.

Liu Jihai，Gu YingS0ng，Yang ZhiChun.InnuenCe ofsupport stiffness on modal and body freedom nuttercharacteristics of a nying wing model[J].Joumal ofVibration Engineering，2018，31（5）：727 733.

[14]黃超.柔性飛翼機(jī)顫振主動(dòng)抑制系統(tǒng)建模、設(shè)計(jì)與驗(yàn)證[D].北京：北京航空航天大學(xué)，2018.

Huang Chao.Modeling，design，and verification of ac—tive nutter suppression system acting on flexible nyingwing aircraft[D].Beijing：Beihang University，2018.

[15]Huang Cha0，Wu Zhigang，Yang Cha0，et al.Flutterboundary prediction for a nying wing model exhibitingbody freedom nutter[C].58th AIAA/ASCE/AHS/ASC Structures，Stmctural Dynamics，and MaterialsConference 2017.2017：1 11.

[16]唐煒，史忠科.飛機(jī)顫振試飛數(shù)據(jù)處理[M].北京：國(guó)防工業(yè)出版社，2018.

[17]郭洪濤，路波，呂彬彬，等.面向紊流激勵(lì)響應(yīng)的顫振亞臨界預(yù)測(cè)方法[C].第13屆中國(guó)系統(tǒng)仿真技術(shù)及其應(yīng)用學(xué)術(shù)年會(huì)，中國(guó)黃山，2011.