基于改進(jìn)多目標(biāo)遺傳算法的艦尾紊流模擬方法

陶楊，韓維

(1.海軍航空工程學(xué)院 飛行器工程系，煙臺(tái)264001;2.海軍裝備研究院，上海200436)

艦尾紊流的特性分析一直以來(lái)都是艦載機(jī)飛行仿真研究中的重點(diǎn)問(wèn)題.艦載機(jī)回收時(shí)，最大著艦誤差源自艦尾紊流對(duì)艦載機(jī)的擾動(dòng)，美軍軍用規(guī)范MIL-F-8785C［1］通過(guò)對(duì)大量實(shí)驗(yàn)和實(shí)測(cè)數(shù)據(jù)的總結(jié)、統(tǒng)計(jì)和分析后，把這種擾動(dòng)視為4種成分的合成并對(duì)它們進(jìn)行了定量描述，本文將其中的一個(gè)分量——自由大氣紊流分量作為研究對(duì)象.

目前多數(shù)紊流處理方法都是將Gauss白噪聲作為輸入信號(hào)，經(jīng)過(guò)紊流頻譜函數(shù)分解成的傳遞函數(shù)濾波后得到相應(yīng)的紊流序列［2-4］.這種方法使用簡(jiǎn)單方便，但生成紊流序列的相關(guān)特性不能很好符合紊流模型相關(guān)函數(shù)的理論表達(dá)式，且它們之間的誤差與采樣步長(zhǎng)大小成反比，采樣步長(zhǎng)越小誤差越大.

針對(duì)這一問(wèn)題，本文基于文獻(xiàn)［5-6］將所求大氣紊流序列先離散化再通過(guò)數(shù)值方法求解的思想對(duì)艦尾紊流進(jìn)行了研究，提出了一種艦尾紊流數(shù)值模擬的新方法.該方法將3個(gè)方向的艦尾紊流離散成Euler前向差分格式進(jìn)行求解，同時(shí)針對(duì)上述方法存在計(jì)算精度低的問(wèn)題，引入了修正系數(shù)，并將紊流相關(guān)性檢驗(yàn)的準(zhǔn)則看作待優(yōu)化目標(biāo)函數(shù)，使用改進(jìn)的多目標(biāo)遺傳算法進(jìn)行修正系數(shù)求解，進(jìn)而實(shí)現(xiàn)在小步長(zhǎng)情況下獲得較精準(zhǔn)的紊流序列.

1 自由大氣紊流生成方法

式(1)所示為MIL-F-8785C根據(jù)統(tǒng)計(jì)學(xué)分析給出的自由大氣紊流分量的頻譜函數(shù):

式中Ω為角頻率.將式(1)的頻譜函數(shù)根據(jù)分子分母的不同形式分為兩類(lèi):第1類(lèi)包括Φu(Ω)和Φw(Ω);第2類(lèi)為Φv(Ω).兩者的處理方法略為不同，下面分別進(jìn)行討論.

1.1 第1類(lèi)頻譜函數(shù)

因Φu(Ω)和Φw(Ω)之間存在線性關(guān)系，故以Φw(Ω)為例說(shuō)明進(jìn)行.Φw(Ω)共軛分解后得到的傳遞函數(shù)，由 Euler前向差分格式可得紊流速度:

式中{r}表示均值為0、方差為1的Gauss白噪聲序列.這里借鑒文獻(xiàn)［7］中的思想加入修正系數(shù)，差分方程變?yōu)?/p>

式中A，B為待定修正系數(shù).由遞推公式(3)可得

因?yàn)閣均值為0［5］，所以其方差為

根據(jù)文獻(xiàn)［5］，可知 AP＜1，則式(6)的等比級(jí)數(shù)和收斂:

在隨機(jī)過(guò)程中，頻譜函數(shù)與相關(guān)函數(shù)之間存在Fourier變換對(duì)的關(guān)系.對(duì)Φw(Ω)作Fourier逆變換可得到該方向紊流序列以距離差ξ為自變量的距離相關(guān)函數(shù):

取 ξ=nh，有

根據(jù)相關(guān)函數(shù)的定義，可知:

對(duì)比式(9)和式(10)可得

將式(11)代入紊流速度方差表達(dá)式(7)可得

1.2 第2類(lèi)頻譜函數(shù)

對(duì)于第2類(lèi)頻譜函數(shù)，可將其因式分解表示為兩個(gè)獨(dú)立頻譜函數(shù)和的形式，之后采用同1.1節(jié)類(lèi)似的方法分別處理.

由此得到傳遞函數(shù):

總的紊流序列{v}可先分別求出各自紊流序列{v1}，{v2}，再求和得到{v}={v1}+{v2}.類(lèi)似式(3)，有

式中，P1=e－h(huán)/1000/A1;P2=e－3h/400/A2;A1，B1，A2，B2為待定修正系數(shù).

取ξ=nh，總的相關(guān)函數(shù):

紊流序列{v}的方差:

該頻譜函數(shù)的相關(guān)函數(shù)為

若式(16)和式(18)兩者表示同樣的衰減率，則比較可得

式(19)必然滿(mǎn)足:

可解出

式中

另外由式(17)可得

這樣，若待定修正系數(shù)確定，即可根據(jù)不同形式的紊流頻譜函數(shù)不斷遞推生成所需的紊流序列.

2 待修正系數(shù)確定方法

艦尾紊流作為大氣紊流的一種特殊存在形式，它擁有大氣紊流的所有特征.相關(guān)性函數(shù)因?yàn)槟軌蚍从吵鑫闪鞯娜啃誀?，所以一直以?lái)多用為檢驗(yàn)和分析紊流的主要指標(biāo)［7-9］.

另外，艦尾紊流也是一個(gè)隨機(jī)過(guò)程，需要用統(tǒng)計(jì)特性對(duì)其進(jìn)行描述.而均方差作為隨機(jī)過(guò)程重要的統(tǒng)計(jì)特性，它反映了隨機(jī)過(guò)程的能量，同時(shí)也是反映數(shù)據(jù)離散程度最常用的一種量化形式，是檢驗(yàn)仿真結(jié)果精確與否的重要指標(biāo).

因此，為檢驗(yàn)生成的紊流序列{x}是否正確，這里采取了兩個(gè)驗(yàn)證準(zhǔn)則:一是讓{x}的均方差σx趨近理論值 σx，th，即 f1(x)二是通過(guò)對(duì){x}進(jìn)行相關(guān)性檢驗(yàn)使相關(guān)函數(shù)Rx(ξ)貼近理論相關(guān)函數(shù) Rx，th(ξ)，即 f2(x)=

在兩個(gè)準(zhǔn)則之間函數(shù)關(guān)系不明確的情況下，不能通過(guò)權(quán)重關(guān)系將之合并為一個(gè)待優(yōu)化函數(shù)用常規(guī)的單目標(biāo)優(yōu)化尋得最優(yōu)解，因此采用多目標(biāo)尋優(yōu)的思想進(jìn)行求解.則問(wèn)題可描述為f(x)=min［f1(x)，f2(x)］，兩個(gè)約束準(zhǔn)則作為待優(yōu)化的目標(biāo)函數(shù)，待修正系數(shù)為通過(guò)尋優(yōu)最終得到的決策變量.

2.1 改進(jìn)的NSGA2算法

遺傳算法基于群體操作的特點(diǎn)使其能在一次運(yùn)算后并行得到多個(gè)解，因此非常適合作為有多個(gè)解的多目標(biāo)優(yōu)化問(wèn)題［10］.目前，比較有代表性的多目標(biāo)遺傳算法包括向量估計(jì)遺傳算法VEGA［11］、多目標(biāo)進(jìn)化算法 MOGA［12］、小生境Pareto遺傳算法 NPGA［13］、非劣排序遺傳算法NSGA［14］、強(qiáng) 度 Pareto 進(jìn) 化 算 法 SPEA［15］、NSGA2［16］等，其中 NSGA2 算法因其低時(shí)間復(fù)雜性、高魯棒性等特點(diǎn)應(yīng)用較為廣泛.本文采用的尋優(yōu)方法就是一種改進(jìn)的NSGA2算法，在此將改進(jìn)的部分加以說(shuō)明.

2.1.1 基于序值的交叉、變異操作

將支配關(guān)系引入遺傳機(jī)制中，位于不同前端的個(gè)體，序值越小表明該個(gè)體越優(yōu)秀，也越有希望能遺傳到下一代.因此在進(jìn)行個(gè)體之間的交叉操作時(shí)，希望能有更多的低序值個(gè)體的基因參與，而自身變異時(shí)則希望能更多地保留低序值個(gè)體的基因.基于以上論斷，標(biāo)準(zhǔn)模擬二進(jìn)制交叉操作和多項(xiàng)式變異操作［17］可改進(jìn)為以下形式.

交叉操作:

式中rank為序值.

式中，r為［0，1］之間均勻分布的隨機(jī)數(shù);ηc為交叉分布指數(shù).

變異操作:

2.1.2 基于擁擠距離的修剪方法

NSGA2中，采用的是精英選擇的策略，低序值的個(gè)體會(huì)被優(yōu)先保留.這樣做存在一定的不足，一旦算法陷入了局部最優(yōu)，因?yàn)樾藜舻舻亩际歉咝蛑档膫€(gè)體，而自身種群沒(méi)有新個(gè)體更新，因此很難跳出.為了保持種群的多樣性，避免算法早熟現(xiàn)象，這里對(duì)修剪方法進(jìn)行改進(jìn).假設(shè)每一前端的個(gè)體由兩部分構(gòu)成，其中原前端按擁擠距離保留一定比例的個(gè)體，其余的由下一前端的個(gè)體補(bǔ)上，候補(bǔ)的個(gè)體同樣遵循按擁擠距離從大到小的選取方法.

2.2 基于改進(jìn)NSGA2算法的參數(shù)選擇

算法流程如圖1所示，具體步驟如下.

步驟1 隨機(jī)初始化種群，種群中每個(gè)個(gè)體由待修正系數(shù)xi(i=1，2，…，n)和其對(duì)應(yīng)的待優(yōu)化目標(biāo)函數(shù)值fj(j=1，2)兩部分組成.

步驟2 根據(jù)目標(biāo)函數(shù)值fj判斷種群個(gè)體間的支配關(guān)系，通過(guò)非支配排序得到每個(gè)個(gè)體的序值和擁擠距離.

步驟3 依據(jù)個(gè)體序值和擁擠距離，采用錦標(biāo)賽選擇方法得到父代種群.

步驟4 保存父代種群到臨時(shí)記憶庫(kù)，同時(shí)開(kāi)始遺傳進(jìn)化過(guò)程，對(duì)父代種群進(jìn)行基于序值的交叉和變異操作，得到子代種群.

步驟5 將臨時(shí)記憶庫(kù)中的父代與步驟4產(chǎn)生的子代合并為中間種群.

步驟6 重復(fù)步驟2的操作，得到該中間種群個(gè)體的序值和擁擠距離.

步驟7 采用改進(jìn)的修剪方法整理中間種群并得到最終的種群.

步驟8 判斷是否達(dá)到最大遺傳代數(shù)，若滿(mǎn)足，結(jié)束優(yōu)化，并根據(jù)需求選擇修正系數(shù);否則，返回步驟2，同時(shí)清空臨時(shí)記憶庫(kù).

圖1 修正系數(shù)選擇流程圖Fig.1 Flow chart of correction coefficient selection

3 計(jì)算結(jié)果及分析

為方便比較，使用本文的改進(jìn)NSGA2算法計(jì)算時(shí)所取的各向紊流采樣步長(zhǎng)為0.1 m，產(chǎn)生白噪聲的隨機(jī)數(shù)種子為23 341，修正系數(shù)限定在［0.5，2.5］之間，優(yōu)化前的原始修正系數(shù)均為 1.通過(guò)上述方法求出的Pareto非劣解分布如圖2所示，可以看出兩個(gè)約束準(zhǔn)則存在著此消彼長(zhǎng)的關(guān)系，并不能同時(shí)達(dá)到最優(yōu)，因此需根據(jù)實(shí)際工程需要選擇其中的一個(gè)解作為該問(wèn)題的最終解.

圖2 x，y，z向 Pareto解分布圖Fig.2 Distribution image of Pareto solution in x，y，z direction

本文以相關(guān)性誤差f2為主需求同時(shí)兼顧紊流方差誤差f1為原則選取的優(yōu)化值及對(duì)應(yīng)的修正系數(shù)如表1所示，圖3為相應(yīng)的各向紊流相關(guān)性檢驗(yàn)結(jié)果.可見(jiàn)，優(yōu)化后的目標(biāo)函數(shù)值均比未優(yōu)化的結(jié)果小一個(gè)數(shù)量級(jí)以上，并且優(yōu)化后紊流相關(guān)性函數(shù)與理論值基本吻合，達(dá)到了良好的效果.采用相應(yīng)修正系數(shù)得出的各方向紊流速度如圖4所示，因?yàn)閤向和z向頻譜函數(shù)存在線性關(guān)系，所以結(jié)果也較類(lèi)似.

表1 優(yōu)化結(jié)果Table1 Optimization results

圖3 x，y，z向相關(guān)性檢驗(yàn)Fig.3 Correlation tests in x，y，z direction

圖4 x，y，z向紊流序列Fig.4 Turbulence sequence in x，y，z direction

4 結(jié)論

本文對(duì)標(biāo)準(zhǔn)NSGA2算法進(jìn)行改進(jìn)，提出了一種基于改進(jìn)NSGA2算法的艦尾紊流數(shù)值模擬方法，通過(guò)仿真實(shí)驗(yàn)，得到了以下結(jié)論:

1)本文算法得到的Pareto非劣解集在解空間的分布均勻且連續(xù)，算法的分布性較好.

2)仿真得到的艦尾紊流序列的均方差誤差和相關(guān)函數(shù)誤差均在要求誤差范圍之內(nèi)，滿(mǎn)足該模型頻譜函數(shù)的理論表達(dá)式，表明通過(guò)該方法生成的艦尾流紊流序列的準(zhǔn)確性高，而且在小步長(zhǎng)情況下也能夠得到較精確的修正系數(shù)，驗(yàn)證了本文方法的正確性和合理性.

References)

［1］ Moorhouse D J，Woodcock R J.Background information and user guide for MIL-F-8785C，military specification-flying qualities of piloted airplanes，AD-A11942［R］.1982.

［2］ 胡國(guó)才，王奇，劉湘一，等.艦尾流對(duì)艦載機(jī)著艦軌跡和動(dòng)態(tài)響應(yīng)的影響研究［J］.飛行力學(xué)，2009，27(6):18-21.Hu G C，Wang Q，Liu X Y，et al.Influence of carrier airwake on carrier-based aircraft landing trajectory and dynamic response［J］.Flight Dynamics，2009，27(6):18-21(in Chinese).

［3］ 蔣康博，劉超，袁東.近艦區(qū)風(fēng)場(chǎng)建模與著艦仿真分析［J］.飛行力學(xué)，2010，28(6):11-15.Jiang K B，Liu C Yuan D.Close-carrier-area wind field modeling and carrier-landing simulation analysis［J］.Flight Dynamics，2010，28(6):11-15(in Chinese).

［4］ 呂開(kāi)東，李新飛，姜邁，等.艦載機(jī)著艦過(guò)程的艦尾氣流場(chǎng)數(shù)值仿真分析［J］.飛行力學(xué)，2013，31(1):18-23.Lü K D，Li X F，Jiang M，et al.Simulation analysis on carrier landing disturbance model［J］.Flight Dynamics，2013，31(1):18-23(in Chinese).

［5］ 趙震炎，肖業(yè)倫，施毅堅(jiān).Dryden大氣紊流模型的數(shù)字仿真技術(shù)［J］.航空學(xué)報(bào)，1986，7(5):433-443.Zhao Z Y，Xiao Y L，Shi Y J.A digital simulation technique for Dryden atmospheric turbulence model［J］.Acta Aeronautica et Astronautica Sinica，1986，7(5):433-443(in Chinese).

［6］ 屈香菊，李勇.一種改進(jìn)的紊流風(fēng)模型及其仿真算法［J］.系統(tǒng)仿真學(xué)報(bào)，2004，16(1):10-13.Qu X J，Li Y.An improved model of atmospheric turbulence and its simulation algorithm［J］.Journal of System Simulation，2004，16(1):10-13(in Chinese).

［7］ 馬東立.大氣紊流數(shù)字仿真的改進(jìn)方法［J］.北京航空航天大學(xué)學(xué)報(bào)，1990，16(3):57-63.Ma D L.An improvement of the digital simulation method for atmospheric turbulence［J］.Journal of Beijing University of Aero-nautics and Astronautics，1990，16(3):57-63(in Chinese).

［8］ 吳揚(yáng)，姜守達(dá).非質(zhì)點(diǎn)飛行器模型的大氣紊流仿真［J］.沈陽(yáng)工業(yè)大學(xué)學(xué)報(bào)，2010，32(1):22-26.Wu Y，Jiang S D.Simulation of atmospheric turbulence for nonparticle aircraft model［J］.Journal of Shenyang University of Technology，2010，32(1):22-26(in Chinese).

［9］ Reid L D.Correlation model for turbulence along the glide path［J］.Journal of Aircraft，1978，15(1):13-20.

［10］ 申曉寧，郭毓，陳慶偉，等.一種保持群體多樣性的多目標(biāo)遺傳算法［J］.控制與決策，2008，23(12):1435-1440.Shen X N，Guo Y，Chen Q W，et al.Multi objective optimization genetic algorithm keeping diversity of population［J］.Control and Decision，2008，23(12):1435-1440(in Chinese).

［11］ Schaffer J D.Multiple objective optimizations with vector evaluated genetic algorithms［C］//Proceedings of the 1st International Conference on Genetic Algorithms.New Jersey:Lawrence Erlbaum Associates，1985:93-100.

［12］ Fonseca C M.Multi-objective genetic algorithms with application to control engineering problems［D］.UK:The University of Shefield，1995.

［13］ Horn J，Nafpliotis N，Goldberg D E.A niched pareto genetic algorithm for multi-objective optimization［C］//Proceedings of the 1st IEEE Conference on Evolutionary Computation.New York:IEEE，1994:82-87.

［14］ Srinivas N，Deb K.Multi-objective function optimization using non-dominated sorting genetic algorithms［J］.Evolutionary Computation，1995，2(2):221-248.

［15］ Zitzler E，Thiele L.Multi-objective evolutionary algorithms:a comparative case study and the strength Pareto approach［J］.IEEE Trans on Evolutionary Computation，1999，3(4):257-271.

［16］ Deb K，Pratap A，Agarwal R B.A fast and elitist multi-objective genetic algorithms:NSGA2［J］.IEEE Transactions on Evolutionary Computation，2002，6(2):182-197.

［17］ Deb K，Agrawal R B.Simulated binary crossover for continuous search space［J］.Complex Systems，1995，9(3):1-15.