馬 創(chuàng)，黃江濤，*，劉 剛，陳 憲，舒博文，，陳其盛，高正紅

(1.中國空氣動力研究與發(fā)展中心，綿陽 621000；2.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072)

0 引 言

聲爆是超聲速飛機(jī)飛行中特有的氣動聲學(xué)現(xiàn)象[1-2]。飛行器在超聲速飛行時，其近場產(chǎn)生復(fù)雜的激波-膨脹波波系，這些波系傳遞到地面，形成空間內(nèi)N形的聲壓分布[3]，這就是聲爆現(xiàn)象。聲爆現(xiàn)象會對生物體和建筑物產(chǎn)生不良影響，因此早期的超聲速客機(jī)只能在海洋上空進(jìn)行超聲速飛行。20世紀(jì)中后期，英法和蘇聯(lián)研制的以“協(xié)和”號和“圖-144”為代表的第一代超聲速客機(jī)都因聲爆過強(qiáng)而被許多國家禁止在境內(nèi)飛行[4]，嚴(yán)重影響了其商業(yè)運(yùn)營，最終以失敗告終。聲爆評估與抑制是超聲速民機(jī)發(fā)展必須解決的卡脖子問題[5]。

聲爆抑制涉及空氣動力學(xué)與聲學(xué)等多個學(xué)科，目前主要途徑有流動控制（靜音錐[5]、熱流注入[6]、矢量發(fā)動機(jī)[7]等）、新概念低聲爆布局（雙向飛翼[8-9]布局、可變前掠翼布局[10]、雙翼布局[11]等）、氣動外形優(yōu)化（代理優(yōu)化[12-14]、梯度優(yōu)化[15-17]、進(jìn)化算法[18]等）……研究表明，在各類方法中，飛機(jī)氣動外形設(shè)計是聲爆抑制的最有效途徑之一。其中，通過等效面積分布（equivalent area, EA）指導(dǎo)氣動外形優(yōu)化[19]是開展低聲爆氣動優(yōu)化設(shè)計的一項關(guān)鍵技術(shù)。

等效面積分布是沿機(jī)身軸線的體積截面積分布與升力分布的疊加，直接決定了超聲速飛行器的聲爆特性，而飛行器等效面積分布特征較大程度由近場過壓決定[20]。但是，通過近場過壓反設(shè)計抑制遠(yuǎn)場聲爆，缺乏直接的遠(yuǎn)場感知聲壓級作為指導(dǎo)。因此，探索有效的遠(yuǎn)場聲爆信號向近場過壓反演方法是開展聲爆抑制研究的關(guān)鍵環(huán)節(jié)。

傳統(tǒng)的低聲爆反設(shè)計主要基于JSGD（Jones-Seebass-George-Darden）方法開展。20世紀(jì)70年代末，Jones、Seebass、George、Darden 四人[21-26]首次提出基于超聲速線化理論的JSGD方法，為低聲爆反設(shè)計奠定了理論基礎(chǔ)。該方法在NASA的高速民用運(yùn)輸計劃（high speed civil transport, HSCT）中得到廣泛應(yīng)用[27]。國內(nèi)馮曉強(qiáng)等使用聲爆預(yù)測算法與遺傳算法對JSGD參數(shù)進(jìn)行了優(yōu)化[28]。但JSGD方法存在試湊的問題，且基于線化理論，對于復(fù)雜構(gòu)型的計算誤差較大，屬于反向等效面積分布設(shè)計方法[29]。因此，有必要探索一種正向等效面積分布設(shè)計方法。

本文擬對近場聲爆信號反演方法的“逆向傳播+POD（proper orthogonal decomposition）”和“逆向傳播+伴隨方程”兩種實現(xiàn)途徑開展研究，以期為等效面積分布指導(dǎo)的超聲速低聲爆飛行器氣動優(yōu)化設(shè)計提供技術(shù)支持。

1 聲爆信號的逆向傳播

以零為初始值開始近場信號反演，無論是POD反演方法還是伴隨方程反演方法，都將導(dǎo)致龐大的計算量，尤其對于POD方法，所需的樣本空間極大。因此獲取可靠的初始波形是進(jìn)一步準(zhǔn)確反演的關(guān)鍵。

遠(yuǎn)場聲爆信號的逆向傳播可以得到近場聲爆信號大致波形，為進(jìn)一步細(xì)節(jié)反演提供有力技術(shù)支撐。

1.1 聲信號逆向傳播模型

本文通過逆向增廣Burgers方程來建立聲爆信號的逆向傳播模型。推導(dǎo)過程如下。

經(jīng)典增廣Burgers方程的無量綱形式為：

等號右端項依次對應(yīng)聲信號在大氣中傳播的非線性效應(yīng)、經(jīng)典耗散、非均勻介質(zhì)和分子弛豫現(xiàn)象。式中：τ為無量綱時間，參考時間為1 / ω0，由近場聲爆信號的采樣頻率決定；Γ為無量綱氣體耗散系數(shù)； θv為無量綱分子松弛時間；Cv為無量綱松弛系數(shù)。

在經(jīng)典增廣Burgers方程中，幾何擴(kuò)散項和大氣分層項取決于聲管傳播路徑與面積，與時間相關(guān)的項只有非線性效應(yīng)項、經(jīng)典耗散項和分子弛豫項。建立逆向Burgers方程即對這幾項進(jìn)行修改，令時間反向流動[30]。

描述聲信號逆向傳播的逆向增廣Burgers方程為：

1.2 逆向增廣Burgers方程求解

逆向增廣Burgers方程數(shù)值求解時，對過壓信號進(jìn)行時間離散，對聲管傳播路徑進(jìn)行空間離散，并對相關(guān)參量進(jìn)行無量綱化處理。使用算子分裂法[31]對式（2）進(jìn)行求解，可分裂成以下五項：

求解過程中，在當(dāng)前高度依次更新無量綱壓力P，空間推進(jìn)直至目標(biāo)高度，完成求解[32]。

式（3）可整理為如下形式：

式（9）～式（11）給出了式（8）的離散方法：

式（4）經(jīng)典耗散項和式（5）分子弛豫項均為病態(tài)方程。為穩(wěn)定求解，引入偽拋物線型正則化方法進(jìn)行修正，修正形式如下：

依照工程經(jīng)驗選擇 ε =1×10-4。引入式（12）所示的偽拋物型方程改寫經(jīng)典耗散項和分子弛豫項。

經(jīng)典耗散項可以寫成：

分子弛豫項可以寫成：

其中，α取0.5。

改寫后的經(jīng)典耗散項和分子弛豫項均為三對角方程。使用TDMA方法[33]進(jìn)行求解。

式（6）和式（7）為守恒方程，求解方法如下：

2 基于POD理論的聲爆信號細(xì)節(jié)反演

由于正向傳播與逆向傳播均存在一定耗散，導(dǎo)致逆向傳播結(jié)果中局部激波信號丟失，需要進(jìn)一步對聲爆信號進(jìn)行精細(xì)化反演計算，以盡可能還原真實波形的細(xì)節(jié)。

本征正交分解是一種基于快照來描述高維空間的方法，其衍生的Gappy POD方法可進(jìn)一步修復(fù)此高維空間內(nèi)的任一殘缺快照，因而可以用于反演近場聲爆信號細(xì)節(jié)。

將設(shè)計好的目標(biāo)遠(yuǎn)場波形作為殘缺快照中的已知部分，待求的近場聲爆信號作為殘缺快照中的未知部分，通過對殘缺快照的修復(fù)來實現(xiàn)對遠(yuǎn)場聲爆信號的反演。

2.1 本征正交分解（POD）方法

POD本質(zhì)上是一種降維方法。對于近遠(yuǎn)場樣本組成的快照集 {Yi，i=(1,2,···,N)}(Yi∈Rn)，POD可以從中尋找主要模態(tài)作為基模態(tài)使所有快照在各階基模態(tài)組成的子空間上正交投影最大，表達(dá)式為：

然后，計算快照的脈動矩陣P：

采用奇異值分解（singular value decomposition, SVD）方法進(jìn)行特征值分解。奇異值分解計算效率高于傳統(tǒng)特征值分解方法，且在計算高階模態(tài)時更精確。對矩陣P進(jìn)行奇異值分解可得：

其中：U∈RN×N，V∈RN×N。 Λ ∈RN×n且僅有對角線元素Λii= σi(σ1≥ σ2≥ ···≥ σr≥0)非零，并得到各階模態(tài)：

任意一個快照（如第i個）都可以由各階模態(tài)的線性疊加得到，

定義各模態(tài)的能量為其對應(yīng)的奇異值平方，

按照能量原則來選取基模態(tài)，定義能量下限為99.9%，即滿足式（27）的前p階模態(tài)被選擇作為基模態(tài)。

2.2 Gappy POD方法

對于一個殘缺樣本，在POD的基礎(chǔ)上發(fā)展出的Gappy POD方法可對其進(jìn)行修復(fù)。對于一個如下形式的殘缺樣本

其中：K為樣本中的已知元素，在本文中為遠(yuǎn)場聲爆信號值；U為樣本中的未知元素，在本文中為待求的相應(yīng)近場聲爆信號值。

POD方法提取出前p階模態(tài)作為主模態(tài)，按以下形式給出：

基于最小二乘法，使用主模態(tài)對殘缺樣本中的殘缺部分進(jìn)行修復(fù)：

其中， Γ =(γ1,γ2···γL)T為殘缺樣本I在主模態(tài)張成的新的p維空間中的坐標(biāo)。

修復(fù)得到的殘缺樣本中的缺失數(shù)據(jù)為：

2.3 聲爆信號的“逆向傳播+POD”反演方法

前文對POD理論的介紹也表明POD方法對樣本空間的依賴性。由于超聲速飛行器近場激波-膨脹波波系的存在，近場聲爆信號波形大多呈現(xiàn)“N”形起伏。直接使用POD方法進(jìn)行反演，對主模態(tài)的要求過于嚴(yán)苛，難以得到理想結(jié)果。

本文提出了逆向傳播與POD結(jié)合的反演方法。首先進(jìn)行聲爆信號的逆向傳播，得到近場聲爆信號的波形主特征。以此為初始值進(jìn)行POD計算，大大減少了POD的計算成本，僅需構(gòu)造波形細(xì)節(jié)所需的樣本空間，再對逆向傳播的近場波形進(jìn)行擾動取樣即可。

采用“擾動取樣—POD”迭代計算的方法進(jìn)行反演，盡可能準(zhǔn)確地還原聲爆信號的局部激波。

3 基于伴隨方程的聲爆信號細(xì)節(jié)反演

基于聲爆伴隨方程求解梯度信息與差分方法相比，避免了大量重復(fù)進(jìn)行Burgers方程求解的正計算過程，是一種高效的梯度求解方法，可基于逆向傳播結(jié)果反演聲爆信號細(xì)節(jié)。

將目標(biāo)遠(yuǎn)場聲爆信號視為目標(biāo)函數(shù)，近場過壓信號的每個時間離散點(diǎn)對應(yīng)的過壓值作為設(shè)計變量，以逆向傳播得到的近場聲爆信號為初始值進(jìn)行梯度尋優(yōu)，以得到更精細(xì)化的反演結(jié)果。

3.1 聲爆伴隨方程推導(dǎo)與梯度求解

將增廣Burgers方程式（1）進(jìn)行算子分裂后寫成矩陣形式：

式（32）～式（35）分別對應(yīng)氧氣分子（O2）弛豫效應(yīng)、氮?dú)夥肿樱∟2）弛豫效應(yīng)、經(jīng)典耗散、非線性扭曲這四個聲爆信號傳播的物理環(huán)節(jié)。kn為增廣Burgers方程式（1）右端第四項與第五項的乘積，表征幾何擴(kuò)散和大氣分層對聲爆信號幅值的影響。P、q、r、t分別表示求解過程中聲壓的四個中間值。

由上述矩陣形式，構(gòu)造式（36）所示的拉格朗日量，推導(dǎo)聲爆伴隨方程的矩陣形式。

其中：l為目標(biāo)函數(shù)；N為垂直空間網(wǎng)格數(shù)；D表示設(shè)計變量，即每個時間離散點(diǎn)對應(yīng)的過壓值； γ0、 γ1、β、λ分別為式（32）～式（35）表示的四個物理環(huán)節(jié)對應(yīng)的伴隨變量。

式（36）使用鏈?zhǔn)角髮?dǎo)法則對設(shè)計變量D求導(dǎo)，可得：

令所有狀態(tài)變量對設(shè)計變量的全導(dǎo)數(shù)為0，可以得到聲爆伴隨方程。

式（39）為聲爆伴隨方程，可進(jìn)行迭代求解。將迭代結(jié)果代入式（38），可求得目標(biāo)遠(yuǎn)場聲爆信號對當(dāng)前近場聲爆信號設(shè)計變量的梯度信息。

3.2 聲爆信號的“逆向傳播+伴隨方程”反演方法

絕大多數(shù)氣動設(shè)計問題都可抽象為優(yōu)化問題。本文的對聲爆信號反演的研究，可視為如下所示的單變量最優(yōu)化問題：

其中：Ftarget為目標(biāo)遠(yuǎn)場過壓信號；Ni為當(dāng)前近場過壓信號，以逆向傳播所得的過壓信號為初始值；FNi為 相應(yīng)的遠(yuǎn)場過壓信號；EM為歐氏距離；Ii為優(yōu)化目標(biāo)函數(shù)，只與當(dāng)前一輪正計算的遠(yuǎn)場聲爆信號有關(guān)。

通過求解聲爆伴隨方程更新目標(biāo)函數(shù)對設(shè)計變量的梯度信息，基于序列二次規(guī)劃算法進(jìn)行梯度尋優(yōu)，可得近場過壓信號的反演結(jié)果。

4 基于LM1021標(biāo)模的聲爆信號反演技術(shù)驗證

采用第二屆聲爆會議SBPW（the second sonic boom prediction workshop）提供的LM1021標(biāo)模算例[34]對提出的方法進(jìn)行可信度驗證。LM1021外形示意于圖1。

4.1 近場聲爆信號與遠(yuǎn)場預(yù)測

在飛行器流場中，機(jī)身周向角0°的近場區(qū)域提取壓力值，并計算得到過壓值。近場過壓信號在空間內(nèi)的分布如圖2所示。

圖2 近場聲爆信號參考值Fig.2 Reference of near-field sonic boom signal

圖3給出了聲爆預(yù)測程序得到的遠(yuǎn)場聲爆信號預(yù)測結(jié)果。由于預(yù)測程序中進(jìn)行了時空變換，遠(yuǎn)場聲爆信號為地面某一空間位置固定的觀測點(diǎn)處的壓力值隨時間的變化。

圖3 遠(yuǎn)場聲爆信號參考值Fig.3 Reference of far-field sonic boom signal

本文提出的反演方法是對圖3所示的遠(yuǎn)場聲爆信號進(jìn)行反演計算，使反演結(jié)果具有與圖2所示的近場聲爆信號等效的可信度。

4.2 遠(yuǎn)場聲爆信號的逆向傳播

基于逆向增廣Burgers方程求解波形B的逆向傳播結(jié)果，目標(biāo)高度為飛行高度H= 16 747 m，聲爆信號時間離散點(diǎn)數(shù) ω0=50000。為了減少離散格式的數(shù)值耗散，垂直方向網(wǎng)格數(shù)為 5 0000。正則化參數(shù)ξ=1×10-6。

圖4給出了逆向傳播結(jié)果與真實近場過壓分布的對比。經(jīng)過逆向傳播得到的波形，較為準(zhǔn)確地反映了真實近場過壓分布的主要特征。從圖5中可以看出，僅經(jīng)過逆向傳播計算，近場聲爆信號仍然存在較大誤差，精度未達(dá)到設(shè)計要求。逆向傳播誤差較大部分主要由病態(tài)方程正則化造成，需要進(jìn)一步對聲爆反演信號進(jìn)行精細(xì)化計算，以還原被耗散掉的細(xì)節(jié)。

圖4 逆向傳播結(jié)果Fig.4 Result of reverse propagation

圖5 逆向傳播結(jié)果相應(yīng)的遠(yuǎn)場聲爆信號Fig.5 Far-field wave corresponding to near-field wave of reverse propagation

4.3 聲爆信號細(xì)節(jié)反演

聲爆信號的細(xì)節(jié)反演采用兩種方案，分別為POD反演方法與伴隨方程反演方法，均以逆向傳播的結(jié)果為初始值進(jìn)行。

4.3.1 基于POD方法的聲爆信號細(xì)節(jié)反演

首先將逆向傳播得到的近場過壓信號參數(shù)化，以定義設(shè)計域和控制變量。精細(xì)化反演計算的設(shè)計域為強(qiáng)激波區(qū)域，基于三次非均勻有理B樣條（non uniform rational B-spline, NURBS）將設(shè)計域參數(shù)化，定義設(shè)計變量共80個，如圖6所示。

圖6 設(shè)計域與控制點(diǎn)分布Fig.6 Design area and control points distribution

逐個擾動控制變量，擾動量為 ± 0.001，擾動生成的每個近場樣本都通過增廣Burgers方程預(yù)測其對應(yīng)的遠(yuǎn)場波形樣本。一對相應(yīng)的近、遠(yuǎn)場聲爆信號樣本組成完整的快照。圖7給出了擾動生成的快照集。

圖7 取樣得到的快照集Fig.7 Sampling snapshots

通過POD方法可以得到從大到小排列的若干特征值和與之對應(yīng)的模態(tài)。圖8給出了前3階模態(tài)形態(tài)。設(shè)定目標(biāo)能量值為99.9%，經(jīng)計算，前21個模態(tài)可以滿足目標(biāo)能量值的要求，可作為主模態(tài)來描述樣本空間。建立包含目標(biāo)遠(yuǎn)場聲爆信號的殘缺快照，使用GPOD方法修復(fù)。對上述“擾動取樣-POD”過程迭代計算，殘差定義同式（41）中的目標(biāo)函數(shù)Ii，收斂條件為相鄰迭代步的殘差變化不超過1×10-5。共迭代400步收斂。

圖8 POD得到的前3階模態(tài)Fig.8 The first three POD modes

4.3.2 基于伴隨方程的聲爆信號細(xì)節(jié)反演

以逆向傳播結(jié)果作為初始值，其時間離散點(diǎn)對應(yīng)的過壓值為設(shè)計變量，基于聲爆伴隨方程求解目標(biāo)遠(yuǎn)場聲爆信號對設(shè)計變量的梯度信息，并更新設(shè)計變量。使用序列二次規(guī)劃（sequential quadratic programming,SQP）算法尋優(yōu)，收斂條件為相鄰迭代步的殘差變化不超過 1×10-5。

4.3.3 結(jié)果分析

圖9給出了兩種反演方法得到的近場聲爆信號。兩種方法對主膨脹波前的聲爆信號反演精度相當(dāng)；對于主膨脹波后的聲爆信號，“逆向傳播+伴隨方程梯度尋優(yōu)”方法得到的結(jié)果有更真實的物理細(xì)節(jié)。

圖9 反演結(jié)果Fig.9 Reverse results

經(jīng)過細(xì)節(jié)反演得到的近場過壓信號由增廣Burgers方程可得相應(yīng)的遠(yuǎn)場聲爆信號，如圖10所示。兩種方法的可信度良好，在尾部激波區(qū)，“逆向傳播+伴隨方程”反演方法的結(jié)果與參考值高度一致，而POD方法存在一定偏差。

圖10 遠(yuǎn)場聲爆信號驗證Fig.10 Far-field sonic boom signal verfication

遠(yuǎn)場聲爆信號通過快速傅里葉變換（fast Fourier transform, FFT）得到頻域內(nèi)的聲壓級、響度級分布，如圖11所示。圖11表明，僅通過逆向傳播得到的結(jié)果，相應(yīng)的遠(yuǎn)場聲壓級與響度級信號在中高頻區(qū)域誤差較大，原因是聲信號逆向傳播中激波信號耗散嚴(yán)重，表明了使用POD方法與伴隨方法還原波形細(xì)節(jié)的必要性。

圖11 遠(yuǎn)場聲爆信號的頻域特征Fig.11 Frequency domain characteristics of far-field sonic boom signal

“逆向傳播+伴隨方程”方法的反演結(jié)果，其遠(yuǎn)場聲壓級與響度級高頻信號精度優(yōu)于“逆向傳播+POD”方法。聲爆信號的高頻部分為激波，表明伴隨方法對激波信號的反演能力優(yōu)于POD方法。

本文提出的聲爆反演方法旨在對人為設(shè)計的、達(dá)到適航要求的遠(yuǎn)場聲爆信號進(jìn)行反演，這一遠(yuǎn)場聲爆信號主要參考感覺噪聲級（preceived noise level, PNL）進(jìn)行設(shè)計。感覺噪聲級是根據(jù)三分之一倍頻程聲壓級計算出來的表示人耳感覺到的噪聲吵鬧程度的單一數(shù)字[35]，直接表征了聲爆信號對人耳的噪聲影響。反演方法的感覺噪聲級精度驗證如表1所示。

表1 感覺噪聲級對比Table 1 Comparison of PNL

兩種反演方法的感覺噪聲級精度較逆向傳播結(jié)果均有大幅提升，其中，“逆向傳播+伴隨方程”反演方法的精度高于“逆向傳播+POD”反演方法。

4.4 等效面積分布分析

等效面積分布是飛行器聲爆特性的決定性因素，因此也成為低聲爆綜合優(yōu)化設(shè)計的指導(dǎo)函數(shù)。Abel算法可將近場聲爆信號轉(zhuǎn)化為飛行器的等效面積分布（EA），指導(dǎo)飛行器低聲爆氣動優(yōu)化設(shè)計。其數(shù)學(xué)表達(dá)如式（42）所示。

圖12給出了兩種反演方法對應(yīng)的等效面積分布，縱坐標(biāo)為歸一化處理的等效面積。相對于逆向傳播結(jié)果，基于伴隨方程反演得到的等效面積分布均更接近參考值，表明局部激波信號丟失對等效面積分布精度有一定不利影響，本文提出的聲爆信號細(xì)節(jié)反演方法可有效降低這一影響。

圖12 等效面積分布對比Fig.12 Comparison of equivalent area distribution

5 結(jié) 論

本文開展了“逆向傳播+POD”與“逆向傳播+伴隨方程”兩種超聲速飛行器近場聲爆信號反演方法的研究工作，驗證了其在低聲爆反設(shè)計優(yōu)化中的可行性。該方法可進(jìn)一步應(yīng)用于未來超聲速民機(jī)的低聲爆設(shè)計，指導(dǎo)該類飛機(jī)的氣動設(shè)計和優(yōu)化。相關(guān)結(jié)論如下：

1）遠(yuǎn)場聲爆信號頻域內(nèi)的聲壓級、響度級、感知聲壓級分析表明，僅經(jīng)過逆向傳播得到的近場聲爆反演結(jié)果中局部激波信號丟失。而本文提出的兩種聲爆信號細(xì)節(jié)反演方法在逆向傳播結(jié)果的基礎(chǔ)上較好地還原了激波信號，感知聲壓級更精準(zhǔn)。

2）盡管近場聲爆信號波形存在一定誤差，從遠(yuǎn)場聲爆信號頻域分析與人耳感知聲壓級分析結(jié)果來看，基于本文提出的方法開展遠(yuǎn)場聲爆信號的設(shè)計、反演與飛行器氣動外形優(yōu)化設(shè)計是可行的。

3）聲壓級頻域分析結(jié)果表明，伴隨方程反演方法對聲壓中高頻信號（即激波信號）的反演更精準(zhǔn)，因此“逆向傳播+伴隨方程”聲爆信號反演方法比“逆向傳播+POD”聲爆反演方法存在優(yōu)勢。

4）由反演結(jié)果得到的等效面積分布與參考值吻合度較好，表明兩種反演方法能夠為后續(xù)基于等效面積分布指導(dǎo)的優(yōu)化設(shè)計工作提供有力支撐。

下一步將根據(jù)遠(yuǎn)場感知聲壓級設(shè)計聲特性良好的遠(yuǎn)場聲爆信號，基于本文提出的反演方法得到相應(yīng)的近場聲爆信號，進(jìn)一步得到等效面積分布，開展低聲爆伴隨綜合優(yōu)化設(shè)計。