組合優(yōu)化算法在沖擊溫度反演中的應(yīng)用

張寧超, 葉 鑫, 李 多, 解孟其, 王 鵬, 劉福生, 鈔紅曉

1. 西安工業(yè)大學(xué)電子信息工程學(xué)院, 陜西 西安 710021 2. 西南交通大學(xué)高溫高壓物理研究所, 四川 成都 610031 3. 西北機(jī)電工程研究所, 陜西 咸陽(yáng) 712099

引 言

沖擊波物理研究中, 溫度是描述沖擊壓縮下材料狀態(tài)和性質(zhì)不可或缺的物理量, 對(duì)高壓科學(xué)、 地球科學(xué)、 行星科學(xué)和材料科學(xué)具有特殊重要性[1-3]。 同時(shí), 特殊材料輻射光譜溫度反演研究, 對(duì)航空航天、 深空探測(cè)和工業(yè)制造等領(lǐng)域有重要應(yīng)用價(jià)值[4-5]。 在溫度求解過(guò)程中, 多光譜的n個(gè)通道, 構(gòu)成了n+1個(gè)未知量的欠定方程組, 目標(biāo)真實(shí)溫度的獲取, 旨在求解該欠定方程組。 當(dāng)前對(duì)光譜溫度求解方法主要包括發(fā)射率模型假設(shè)、 神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練等數(shù)據(jù)處理方法。 前者存在一定局限性, 只有當(dāng)假定的材料發(fā)射率模型與實(shí)際發(fā)射率相符時(shí)才能準(zhǔn)確的反演目標(biāo)溫度, 面對(duì)未知發(fā)射率的材料, 溫度反演將束手無(wú)策。 早期Flower等[6]利用模型假設(shè)法, 開(kāi)展了較多的實(shí)驗(yàn)與應(yīng)用研究, 雖然得到了一些材料的反演溫度, 但模型的不確定性限制了其推廣。 楊永軍[7]針對(duì)高溫材料表面溫度難以準(zhǔn)確測(cè)量的問(wèn)題, 利用光譜儀測(cè)量材料的能量信息, 通過(guò)多光譜測(cè)溫獲取其表面溫度進(jìn)而計(jì)算出光譜發(fā)射率, 并利用不銹鋼在1 100 K時(shí)的發(fā)射率進(jìn)行了驗(yàn)證。 孫紅勝等[8]針對(duì)彌散介質(zhì)下的材料溫度測(cè)量, 分析并總結(jié)了幾種主要的輻射測(cè)溫方法的優(yōu)缺點(diǎn), 為彌散介質(zhì)下的高溫測(cè)量提供了研究方向。 孫曉剛等[9]較早將遺傳算法與神經(jīng)網(wǎng)絡(luò)相結(jié)合, 將其應(yīng)用于目標(biāo)溫度的測(cè)量與反演, 該方法提高了反演的精度, 但需要大量樣本數(shù)據(jù)作為支撐, 求解溫度前需對(duì)樣本進(jìn)行訓(xùn)練, 所需要的數(shù)據(jù)集大且訓(xùn)練所需要的時(shí)間較長(zhǎng)。 然而在沖擊溫度測(cè)量實(shí)驗(yàn)中, 材料在高溫高壓狀態(tài)下的光譜發(fā)射率模型很難獲得, 尤其在面對(duì)發(fā)射率未知的材料或被測(cè)材料發(fā)生結(jié)構(gòu)相變時(shí), 這種方法將不再適用。 根據(jù)約束優(yōu)化理論, 如果將材料溫度求解轉(zhuǎn)化為目標(biāo)優(yōu)化問(wèn)題, 通過(guò)建立相應(yīng)的溫度反演模型, 利用優(yōu)化算法對(duì)模型進(jìn)行求解, 即可實(shí)現(xiàn)溫度的反演。 根據(jù)測(cè)量特點(diǎn)及應(yīng)用范圍, 戴景民等[10]對(duì)發(fā)射率測(cè)量方法進(jìn)行了分類, 分析了發(fā)射率測(cè)量方面的問(wèn)題, 將人工神經(jīng)網(wǎng)絡(luò)與二次細(xì)分方法相結(jié)合提高了測(cè)溫精度, 但該方法所需發(fā)射率樣本集過(guò)多, 存在一定局限性。 刑鍵等[11]利用約束優(yōu)化對(duì)溫度反演進(jìn)行了系統(tǒng)研究, 利用二次測(cè)量法結(jié)合迭代遞推算法, 提出了無(wú)需假設(shè)發(fā)射率模型的輻射溫度測(cè)量方法, 利用發(fā)射率約束偏差的思想提高了溫度反演的速度, 基于廣義逆-坐標(biāo)變換提高了溫度反演的精度。 根據(jù)應(yīng)用場(chǎng)景不同, 張福才等[12]使用了不同類型的優(yōu)化算法, 分別對(duì)單目標(biāo)和多目標(biāo)模型進(jìn)行了算法仿真, 并對(duì)火箭發(fā)動(dòng)機(jī)運(yùn)行時(shí)尾焰溫度值進(jìn)行了測(cè)量與反演, 溫度反演過(guò)程在效率和精度上都取得進(jìn)展。 沖擊溫度屬于材料內(nèi)部受到高強(qiáng)度壓縮能量在瞬間的釋放, 因此沖擊溫度測(cè)量過(guò)程更快, 對(duì)溫度反演的時(shí)間效率要求更高, 同時(shí)沖擊壓縮下材料的結(jié)構(gòu)極容易發(fā)生相變, 更不能保證發(fā)射模型的穩(wěn)定性與唯一性。 為此本工作基于多光譜測(cè)量沖擊溫度, 將測(cè)溫方法與優(yōu)化組合理論相結(jié)合, 避免目標(biāo)發(fā)射率模型對(duì)溫度反演的影響, 提出將平衡優(yōu)化器算法(equilibrium optimizer, EO)與序列二次規(guī)劃法(sequential quadratic programming, SQP)算法組合使用, 既提高了溫度反演計(jì)算的效率, 也確保了溫度反演的精度和結(jié)果的穩(wěn)定性。

1 實(shí)驗(yàn)部分

1.1 多光譜測(cè)溫

圖1為沖擊溫度的實(shí)驗(yàn)測(cè)量系統(tǒng)。 二級(jí)輕氣炮發(fā)射高速飛片, 產(chǎn)生沖擊波通過(guò)金屬基板, 由于窗口材料一直處于透明狀態(tài), 所以光纖束將記錄沖擊波壓縮下金屬后界面發(fā)光信息。

圖1 沖擊輻射實(shí)驗(yàn)測(cè)試系統(tǒng)Fig.1 Shock radiation test system

圖2是實(shí)驗(yàn)測(cè)量得到的金屬銅后界面發(fā)光強(qiáng)度信息。 飛片撞擊材料沖擊波傳到金屬后界面時(shí), 由于透明窗口與金屬基板間存在空氣間隙,t1時(shí)刻出現(xiàn)發(fā)光尖峰, 說(shuō)明此刻光纖開(kāi)始接收金屬基板后界面的輻射信息,t2時(shí)刻輻射強(qiáng)度出現(xiàn)明顯拐折, 這說(shuō)明輻射受邊側(cè)稀疏波或者追趕稀疏波的影響。

圖2 金屬銅輻射強(qiáng)度曲線Fig.2 Copper impact radiation intensity curve

1.2 多光譜測(cè)溫的溫度反演模型

由普朗克黑體輻射定律可知,n個(gè)通道的輻射高溫計(jì)中,當(dāng)材料真實(shí)溫度為T, 黑體參考溫度為T′時(shí), 材料溫度與參考溫度之間的關(guān)系如式(1)

(1)

由式(1)結(jié)合輻射測(cè)溫理論, 理想狀態(tài)下每個(gè)通道求解得到的溫度都應(yīng)相同, 因此各通道反演得到的目標(biāo)溫度理論偏差應(yīng)趨近于零, 由此可構(gòu)成目標(biāo)函數(shù), 如式(2)所示

(2)

式(2)中,Ti為各個(gè)通道反演獲取的目標(biāo)溫度,E(Ti)為所有通道實(shí)際反演溫度的平均值。 對(duì)式(1)進(jìn)行變形可得

(3)

(4)

每個(gè)通道獲取到的溫度均值如式(5)

(5)

理想狀態(tài)下, 各通道反演得到的溫度均值與第一個(gè)通道溫度間的差值也應(yīng)為0, 由此可以構(gòu)成等式約束優(yōu)化模型如式(6)所示

T1-E(Ti)=0

(6)

由發(fā)射率0≤ε(λi,T)≤1構(gòu)成不等式約束優(yōu)化模型如式(7)所示

0≤ε(λi,T)≤1

(7)

式(7)中,f(x)為各個(gè)通道溫度值的方差。

2 算法仿真與實(shí)例分析

溫度反演模型建立了黑體輻射溫度、 目標(biāo)溫度、 發(fā)射率與波長(zhǎng)之間的函數(shù)關(guān)系, 該模型可利用優(yōu)化算法求解得出材料的發(fā)射率與物理溫度值。 傳統(tǒng)優(yōu)化算法求解比較依賴初始點(diǎn)且計(jì)算速度較慢, 為了適應(yīng)沖擊溫度測(cè)量的需求, 首先選用平衡優(yōu)化器算法, 該算法具有參數(shù)少、 耗時(shí)少且尋優(yōu)能力強(qiáng)的優(yōu)點(diǎn)。

2.1 基于平衡優(yōu)化器的溫度反演

平衡優(yōu)化器算法(EO算法)是以動(dòng)態(tài)質(zhì)量平衡為靈感而提出的啟發(fā)式優(yōu)化算法。 質(zhì)量平衡方程可用一階微分方程表示為

(8)

式(8)中:V為控制容積;c為控制容積內(nèi)的濃度;Q為控制容積的容積流速;ceq為控制容積平衡狀態(tài)的濃度;G為控制容積內(nèi)的質(zhì)量生成速率。

通過(guò)求解式(8)可得

c=ceq+(c0-ceq)F+G(1-F)/λV

(9)

(10)

t=(1-iter/Maxiter)a2iter/Maxiter

(11)

EO中取α用來(lái)控制勘探能力,a1越大勘探能力越強(qiáng), 開(kāi)采能力將變?nèi)?a2控制開(kāi)采能力,a2越大開(kāi)采能力越強(qiáng), 勘探能力逐漸越弱。 通常a1取值范圍為1～3,a2范圍為0～2。 由于沒(méi)有通用的模型, 為了選取最佳的EO算法參數(shù), 選用了四類常見(jiàn)材料發(fā)射率模型[13-14]對(duì)算法進(jìn)行仿真驗(yàn)證:

模型A:ε(λ)=ea+bλ,

模型B:ε(λ)=1/[a+bln(λ)/λ],

模型C:ε(λ)=a+bλ,

模型D:ε(λ)=a+bλ+cλ2

式中:ε(λ)為發(fā)射率;a、b和c為系數(shù)。

上述四種模型涵蓋了絕大部分材料的發(fā)射率特性, 根據(jù)選用材料沖擊溫度范圍, 設(shè)材料物理溫度為3 000 K, 通過(guò)發(fā)射率數(shù)據(jù)對(duì)EO算法中的參數(shù)進(jìn)行優(yōu)化調(diào)整, 參考溫度選取1 600 K, 具體數(shù)據(jù)如表1所示。

表1 四種模型的發(fā)射率參數(shù)Table 1 Emissivity parameters of four models

設(shè)定a2不變, 在a1取不同數(shù)值時(shí), 利用EO算法對(duì)目標(biāo)溫度進(jìn)行反演, 結(jié)合常見(jiàn)材料的發(fā)射率特征, 將發(fā)射率的范圍限定為0.1≤ε(λ,T)≤0.9, 結(jié)果如表2所示。

表2 不同a1值的四種模型的溫度反演結(jié)果(K)Table 3 Temperature inversion results of four models with different a1 values (K)

對(duì)表2中的溫度數(shù)據(jù)進(jìn)行求解, 獲取不同a1取值下的發(fā)射率數(shù)值。 結(jié)果如圖3所示。 模型A中當(dāng)a1取值為1.8和2.0時(shí)發(fā)射率曲線與理論值最為接近, 當(dāng)a1取值為2.8時(shí), 由于算法的開(kāi)采能力大幅減弱, 導(dǎo)致反演的發(fā)射率與理論值誤差較大; 模型B中當(dāng)a1取值為1.4、 1.8、 2.0和2.4時(shí)得到的發(fā)射率曲線與理論值最為接近, 其中當(dāng)a1為1.4時(shí)的發(fā)射率低于理論值; 模型C中當(dāng)a1取值為2.0和2.4時(shí)得到的發(fā)射率曲線與理論值最為接近,a1為2.0時(shí)的發(fā)射率低于理論值,a1為2.4時(shí)的發(fā)射率高于理論值; 模型D中a1取值為1.8、 2.0、 2.2和2.4時(shí)的發(fā)射率與理論值更為接近, 其中a1為1.8和2.4時(shí)的發(fā)射率高于理論值。

圖3 不同a1值時(shí)四種模型的發(fā)射率趨勢(shì)Fig.3 Emissivity trend of four models at different a1 values

設(shè)定a1不變, 選取a2范圍為0.2～1.8, 步長(zhǎng)為0.2, 對(duì)四種模型的溫度進(jìn)行求解, 結(jié)果如表3所示。

表3 不同a2值的四種模型的溫度反演結(jié)果(K)Table 3 Temperature inversion results of four models with different a2 values (K)

對(duì)表3中的溫度數(shù)據(jù)進(jìn)行求解, 獲取不同a2取值下的發(fā)射率數(shù)值, 如圖4所示。

圖4 不同a2值時(shí)四種模型的發(fā)射率趨勢(shì)Fig.4 Emissivity trends of four models at different a2 values

結(jié)合圖4和表3的數(shù)據(jù), A模型中a2取值為0.8、 1.0與1.8時(shí), 發(fā)射率變化曲線與理論曲線接近; B模型中a2取值為0.4或1.0時(shí), 發(fā)射率值曲線與理論值接近; C類模型a2取值為1.0或1.4時(shí), 溫度反演效果較佳; D類模型a2取值為1.0或1.8時(shí)發(fā)射率與理論值更為貼近。 結(jié)合圖1可得, 當(dāng)a1取1.8,a2取1.0時(shí), EO算法的溫度反演效果最佳。 然而, 在多次求解尋找最優(yōu)值時(shí)發(fā)現(xiàn), EO算法得到的溫度值波動(dòng)較大, 穩(wěn)定性較低。 這是由于EO算法在優(yōu)化過(guò)程中, 每一次迭代會(huì)等概率的隨機(jī)從平衡池中選擇一個(gè)當(dāng)前的最優(yōu)解進(jìn)行迭代。 由于從平衡池中每次選取的當(dāng)前最優(yōu)解存在隨機(jī)性, 使得迭代的尋優(yōu)方向變得隨機(jī), 導(dǎo)致每次計(jì)算出的溫度值之間的浮動(dòng)較大, 穩(wěn)定性偏低。

2.2 基于平衡優(yōu)化器結(jié)合多起點(diǎn)序列二次規(guī)劃的溫度反演

為解決平衡優(yōu)化器求解溫度值穩(wěn)定性偏低, 且誤差較大的問(wèn)題, 將平衡優(yōu)化器算法(equilibrium optimizer, EO)與序列二次規(guī)劃法(sequential quadratic programming, SQP)相結(jié)合, 對(duì)溫度反演方法進(jìn)行完善, 解決EO算法反演求解穩(wěn)定性不足的問(wèn)題。 EO算法運(yùn)算速度快, 收斂性強(qiáng), 但最終求解結(jié)果的穩(wěn)定性較差, 解的準(zhǔn)確性較低, 易陷入局部最優(yōu)。 SQP法邊界搜索能力強(qiáng), 解的穩(wěn)定性高, 但初始點(diǎn)的不同, 會(huì)導(dǎo)致運(yùn)算量增大而影響收斂速度, 增加算法的計(jì)算效率。 結(jié)合兩種算法各自的優(yōu)點(diǎn), 用EO算法求得的解作為SQP算法的初始點(diǎn)之一, 利用SQP方法展開(kāi)多起點(diǎn)迭代尋優(yōu), 實(shí)現(xiàn)多起點(diǎn)尋優(yōu)的方式求解, 保證精度和穩(wěn)定性高于單一算法, 實(shí)現(xiàn)精度高、 速度快且穩(wěn)定性強(qiáng)的溫度反演。

利用SQP方法來(lái)完成溫度求解, 隨機(jī)在可行域范圍內(nèi)初始化多個(gè)起始點(diǎn), 將求得的解進(jìn)行對(duì)比并輸出滿足式(7)的最優(yōu)解, 即為目標(biāo)精度最高的溫度值。 仿真結(jié)果表明起始點(diǎn)個(gè)數(shù)n直接影響溫度值精度與計(jì)算時(shí)間。 分別取3, 6, 9, 12和15為起始點(diǎn), 對(duì)四種發(fā)射率模型進(jìn)行求解, 當(dāng)起始點(diǎn)的個(gè)數(shù)為3或15時(shí), 求得的四種模型對(duì)應(yīng)溫度值誤差都很高, 起始點(diǎn)個(gè)數(shù)取9時(shí), 四種模型誤差較低。 不同起始點(diǎn)個(gè)數(shù)四種模型溫度反演結(jié)果如表4所示。

表4 不同起始點(diǎn)數(shù)n的四種模型的溫度反演結(jié)果Table 4 Temperature inversion results of four models with different starting points n

為了驗(yàn)證本溫度反演方法的有效性, 分別利用EO-SQP組合算法、 基于內(nèi)點(diǎn)罰函數(shù)(interior penalty function, IPF)的對(duì)四種模型溫度進(jìn)行驗(yàn)證, 將兩種方法計(jì)算的精度和時(shí)間進(jìn)行對(duì)比, 表5中a1,a2起始點(diǎn)的個(gè)數(shù)取值分別為1.8, 1.0, 9。

表5 兩種方法溫度反演結(jié)果對(duì)比Table 5 Comparison of temperature inversion results between two methods

結(jié)果表明IPF方法在求解B類模型的溫度時(shí), 誤差較小, 僅為0.39%。 但在求解A、 C和D三類模型時(shí), 溫度誤差普遍較大。 尤其在求解C類模型時(shí), 誤差達(dá)到了2.37%, 且IPF在反演溫度時(shí), 程序迭代運(yùn)算所需要的時(shí)間較長(zhǎng)。 本文算法求解得到的四類模型溫度值, 誤差均小于1%, 且算法的運(yùn)行時(shí)間小于3 s, 相較于IPF, 在提高精度的前提下, 反演效率大幅提升, 求解精度和速度都要優(yōu)于IPF。

由亮度溫度模型, 各通道反演出的光譜發(fā)射率可表示為

(12)

根據(jù)式(12), 對(duì)表5中兩種方法反演出的溫度結(jié)果進(jìn)行計(jì)算, 獲取其對(duì)應(yīng)的發(fā)射率曲線, 并與目標(biāo)發(fā)射率理論值進(jìn)行對(duì)比, 如圖5所示。

圖5 3 000 K時(shí)四種模型發(fā)射率變化曲線Fig.5 Emissivity curves of four models at 3 000 K

圖5中IPF方法和本組合反演算法對(duì)模型A、 B和D都較適用且發(fā)射率反演誤差較小, 但I(xiàn)PF應(yīng)用于C類材料模型的溫度反演時(shí), 發(fā)射率曲線與目標(biāo)發(fā)射率理論值偏差較大, 結(jié)合表5中數(shù)據(jù)可知, 本文提出的組合優(yōu)化算法適用范圍廣。

2.3 實(shí)例分析

利用典型的輻射高溫計(jì), 獲取8波長(zhǎng)目標(biāo)沖擊發(fā)光強(qiáng)度值。 測(cè)試系統(tǒng)各通道的中心波長(zhǎng)分別為: 0.809、 0.779、 0.702、 0.650、 0.589、 0.533、 0.509和0.488 μm。 分別利用最小二乘法(least square method, LSM)、 內(nèi)點(diǎn)罰函數(shù)法(IPF)以及本算法(EO-SQP)對(duì)沖擊過(guò)程各時(shí)刻溫度進(jìn)行求解, 溫度反演結(jié)果如圖6所示。

圖6 三種方法溫度反演結(jié)果Fig.6 Temperature inversion results of three methods

依據(jù)金屬測(cè)溫理論, 實(shí)驗(yàn)直接獲取的輻射對(duì)應(yīng)界面溫度, 并不代表在該壓力下金屬的卸載溫度。 由“理想界面模型”, 金屬卸載溫度、 界面溫度以及窗口沖擊溫度三者間的關(guān)系如式(13)。

(13)

利用式(13)計(jì)算得到本實(shí)驗(yàn)沖擊溫度、 卸載溫度以及界面溫度數(shù)值, 將其與本算法、 最小二乘法以及內(nèi)點(diǎn)罰函數(shù)法反演得到的溫度值進(jìn)行對(duì)比。

表6是本次實(shí)驗(yàn)相關(guān)量的理論值與實(shí)驗(yàn)計(jì)算結(jié)果。 最小二乘法計(jì)算出的界面溫度絕對(duì)誤差約為16.44%, 內(nèi)點(diǎn)罰函數(shù)法計(jì)算出的界面溫度絕對(duì)誤差約為11.37%, 本算法求解出的界面溫度絕對(duì)誤差約為4.41%, 得到的結(jié)果更接近理論計(jì)算值。

表6 三種方法的計(jì)算結(jié)果比較Table 6 Comparison of calculation results of three methods

3 結(jié) 論

要獲得材料真實(shí)的沖擊溫度值, 除了精密的測(cè)試手段與方法, 溫度反演模型與溫度反演算法的選取至關(guān)重要。 針對(duì)以往沖擊溫度測(cè)量方法適用性差, 溫度反演速度慢的問(wèn)題, 將普朗克理論與目標(biāo)優(yōu)化相結(jié)合, 建立基于約束優(yōu)化的溫度反演模型, 并利用平衡優(yōu)化器算法對(duì)模型進(jìn)行求解以獲取材料溫度值, 提高了溫度反演的運(yùn)算效率以及溫度反演方法的適用性; 針對(duì)平衡優(yōu)化器算法求解精度較差且穩(wěn)定性偏低的問(wèn)題, 利用多起點(diǎn)序列二次規(guī)劃法進(jìn)一步完善求解, 避免平衡優(yōu)化器算法陷入局部最優(yōu), 以提高溫度反演的精度和效率。 結(jié)果表明, 本文提出的算法是一種有效的沖擊波壓縮下材料溫度反演算法, 實(shí)現(xiàn)了高效、 高精度以及適用性廣的溫度反演。