杜配冰，劉 鈺，陳志華，關(guān) 奇，夏洪富，蔡利兵

(西北核技術(shù)研究所，西安710024)

激光在大氣中傳輸時會出現(xiàn)光束擴(kuò)展及功率密度下降等一系列湍流效應(yīng)，使激光系統(tǒng)的使用效果發(fā)生變化[1-2]。數(shù)值模擬是研究激光大氣傳輸湍流效應(yīng)的重要手段之一，數(shù)值模擬中常用的方法是采用多層相位屏模擬大氣湍流特性，再通過傅里葉變換求解激光的傳輸過程[3-4]，一般只針對確定的輸入?yún)?shù)進(jìn)行計算。由于現(xiàn)實(shí)條件下激光器與大氣環(huán)境均非絕對穩(wěn)定，激光器及大氣參數(shù)均存在著較大的隨機(jī)不確定性，這嚴(yán)重影響了數(shù)值模擬結(jié)果的可信度。近年來，張建柱等分析了大氣參數(shù)的測量不確定度對激光大氣傳輸?shù)桨衅骄β拭芏鹊挠绊?，結(jié)果表明，當(dāng)湍流越強(qiáng)且熱暈越弱時，大氣相干長度的測量不確定度對計算結(jié)果的影響越大[5]。陳志華等研究了激光大氣傳輸?shù)桨衅骄β拭芏鹊耐馔朴嬎惴椒?，并分析了近場激光功率密度和大氣傳輸修正系?shù)的不確定度，結(jié)果表明，當(dāng)測量路徑上的大氣相干長度大于8 cm且外推路徑上的發(fā)射仰角大于35°時，可以采用幾何外推方法計算到靶平均功率密度[6]。

對激光大氣傳輸湍流效應(yīng)的數(shù)值模擬結(jié)果進(jìn)行不確定度量化，可減少不確定性帶來的風(fēng)險，并有效提高數(shù)值模擬的預(yù)測能力[7-8]。蒙特卡羅(MC)方法是不確定度量化中最常用的方法，具有操作簡單和適用性強(qiáng)等特點(diǎn)[9]。但MC方法的均值收斂率過慢，需要使用大量樣本進(jìn)行數(shù)值模擬才能得到較精確的結(jié)果，因此計算成本過高。近年來，多項(xiàng)式混沌方法在不確定度量化計算中的應(yīng)用越來越廣泛。例如，王瑞利等針對爆轟問題研究了JWL模型的多項(xiàng)式混沌不確定度量化方法[10]；伍月千等基于多項(xiàng)式混沌方法研究了電磁波入射角對孔縫結(jié)構(gòu)機(jī)箱中心點(diǎn)電場強(qiáng)度幅值的影響[11]；萬華平等研究了針對結(jié)構(gòu)動力特性的多項(xiàng)式混沌不確定度量化方法[12]。多項(xiàng)式混沌方法的核心思想是將隨機(jī)解表示為隨機(jī)輸入?yún)?shù)的正交多項(xiàng)式展開形式，從而構(gòu)建與數(shù)值模擬對應(yīng)的多項(xiàng)式混沌代理模型。由于正交多項(xiàng)式序列存在高精度快速計算方法[13]，所以，用代理模型取代數(shù)值模擬可有效降低數(shù)值模擬生成樣本的計算開銷。目前，國內(nèi)外在激光大氣傳輸?shù)牟淮_定度量化方面，尚缺乏多項(xiàng)式混沌方法的研究報道。

本文將多項(xiàng)式混沌方法與拉丁超立方抽樣方法[14]相結(jié)合，建立了激光大氣傳輸湍流效應(yīng)數(shù)值模擬的不確定度量化方法；以地空垂直上行激光傳輸為例，考慮了激光器出口光束質(zhì)量和大氣相干長度的不確定性，計算了遠(yuǎn)場功率密度的不確定度，并與MC方法的計算結(jié)果進(jìn)行了對比，驗(yàn)證了本文不確定度量化方法的精確性與有效性。

1激光大氣傳輸湍流效應(yīng)的數(shù)值模擬策略

激光大氣傳輸?shù)臄?shù)值模擬通常包含真空傳輸和相位屏構(gòu)造2個部分。激光在大氣中的傳輸可用拋物型方程描述為

(1)

其中，k=2π/λ，λ為激光波長；u為入射光波在自由空間傳輸至(x,y,z)處光波電場的復(fù)振幅；n′為湍流引起折射率的起伏值。求解式(1)時，通常將大氣傳輸路徑分成若干段，每一段均設(shè)置薄相位屏以體現(xiàn)湍流的影響，其余部分設(shè)置為真空，結(jié)合傅里葉變換分段循環(huán)求解。

傅里葉變換譜反演法是構(gòu)造相位屏的常用方法，主要是通過湍流的功率譜對隨機(jī)矩陣進(jìn)行濾波，再通過傅里葉逆變換得到大氣擾動的相位[15]。

(2)

其中，κx,κy分別為x，y方向的相空間的波數(shù)；C為常數(shù)；R(κx,κy)為服從N(0,1)的高斯隨機(jī)數(shù)；Fφ(κx,κy)為相位畸變的功率譜密度。

假設(shè)相鄰的相位屏相距Δz，則第i個相位屏處光波電場的復(fù)振幅為

(3)

其中，f表示傅里葉變換；f-1表示傅里葉逆變換；τ(r)為相鄰的相位屏之間由湍流引起的相位調(diào)制函數(shù)。

2多項(xiàng)式混沌方法

多項(xiàng)式混沌方法最早于1938年由Wiener提出[16]。該方法使用Hermite多項(xiàng)式基函數(shù)逼近帶有Gauss型隨機(jī)參數(shù)問題的隨機(jī)解。隨后，Xiu等通過Askey法則將多項(xiàng)式混沌方法推廣到帶有非Gauss型隨機(jī)參數(shù)的問題，可根據(jù)隨機(jī)參數(shù)的不同分布類型找到對應(yīng)l2范數(shù)最優(yōu)逼近的正交多項(xiàng)式基函數(shù)[17]。表1給出了常見概率分布函數(shù)對應(yīng)的正交多項(xiàng)式基函數(shù)及其分布區(qū)間。

表1常見概率分布函數(shù)對應(yīng)的正交多項(xiàng)式基函數(shù)及其定義域Tab.1Corresponding polynomial chaos basis of generalcontinuous probability distributions and their domain

對于帶有n個獨(dú)立隨機(jī)參數(shù)X的不確定性系統(tǒng)Y(X)，若X=(X1,X2,…,Xn) ，隨機(jī)參數(shù)的聯(lián)合概率分布函數(shù)FX(x)=P(X1≤x1,X2≤x2,…,Xn≤xn)，隨機(jī)參數(shù)Xi的邊緣概率分布函數(shù)FXi(xi)=P(Xi≤xi)，則

(4)

(5)

其中，i=(i1,i2,…,in)，|i|=i1+i2+…+in，且0≤|i|,i1,i2,…,in≤N。構(gòu)建一個關(guān)于目標(biāo)函數(shù)y(X)的l2范數(shù)，逼近N階正交多項(xiàng)式展開式，即

(6)

一元正交多項(xiàng)式基函數(shù)的正交性可表示為

E[bi(Xk)bj(Xk)]=

(7)

E[Bi(X)Bj(X)]=

(8)

E[y(X)]≈E[yN(X)]=

(9)

目標(biāo)函數(shù)的方差為

(10)

其中，ci可根據(jù)基函數(shù)的正交性由式(7)得到，即

(11)

采用多項(xiàng)式混沌方法逼近光滑函數(shù)的收斂速度隨N呈指數(shù)增長，構(gòu)建代理模型的計算開銷主要用于求解多項(xiàng)式展開系數(shù)[18]。由于MC方法的收斂速度低于線性收斂速度，所以采用式(11)求解系數(shù)所需的樣本量過高，嚴(yán)重降低了計算效率。因此，可隨機(jī)抽取少量樣本，采用最小二乘回歸方法求解系數(shù)，即

(12)

3算例

本文對激光大氣傳輸湍流效應(yīng)數(shù)值模擬進(jìn)行不確定度量化時，根據(jù)文獻(xiàn)[5]，主要考慮輸入?yún)?shù)帶來的不確定性，以激光器出口光束質(zhì)量β0和大氣相干長度r0為不確定性來源。

選取激光為準(zhǔn)直Gauss光束，大氣能見度和傳輸距離均為10 km，接收面為5 000 mm×5 000 mm的正方形，傳輸路徑為垂直上行傳輸，大氣折射率結(jié)構(gòu)常數(shù)選用Hufnagel-Valley5/7參數(shù)模型，折射率起伏譜選用Kolmogorov模型。結(jié)合不確定性參數(shù)的統(tǒng)計分布隨機(jī)特性[20]，選β0為3倍衍射極限，r0為6.75 cm。根據(jù)文獻(xiàn)[6]，實(shí)驗(yàn)測量大氣相干長度的不確定度應(yīng)不大于7%，所以，本文假設(shè)在高斯分布場合下β0和r0的不確定度均為5%，即β0和r0均以99.73%的概率分別分布于3(1±5%)與6.75(1±5%)的參數(shù)空間中。

用MC方法進(jìn)行激光大氣傳輸數(shù)值模擬不確定度量化時，對β0和r0隨機(jī)抽取10 000組樣本，通過數(shù)值模擬可以得到10 000組遠(yuǎn)場光強(qiáng)分布的樣本，歸一化后得到的遠(yuǎn)場相對光強(qiáng)分布，如圖1所示。

圖1用MC方法得到的10 000組樣本的遠(yuǎn)場相對光強(qiáng)分布Fig.110 000 samples of far-field relative optical intensity distribution calculated by MC method

采用隨機(jī)抽樣方法和拉丁超立方抽樣(LHS)方法對β0和r0隨機(jī)抽取20組樣本，抽樣結(jié)果的對比如圖2所示。其中，圖2(a)為隨機(jī)抽樣方法抽取的20組樣本，圖2(b)為LHS方法抽取的20組樣本，圖2(c)和圖2(d)為隨機(jī)抽樣方法和LHS方法抽取樣本對應(yīng)的累積分布函數(shù)(CDF)，由于LHS方法將CDF等分成若干區(qū)域后再分別進(jìn)行隨機(jī)抽樣，所以抽取的樣本更均勻。

(a)20 samples by random sampling

(b)20 samples by LHS

(c)CDFs of 20 samples by random sampling

(d)CDFs of 20 samples by LHS

抽取樣本后進(jìn)行數(shù)值模擬，根據(jù)多項(xiàng)式混沌(PC)方法和表1，選取Hermite多項(xiàng)式構(gòu)建遠(yuǎn)場功率密度的代理模型，用代理模型代替數(shù)值模擬，將隨機(jī)抽取的β0和r0樣本代入代理模型，最終得到歸一化后遠(yuǎn)場相對功率密度均值，如表2所列。遠(yuǎn)場相對功率密度的不確定度及其與MC方法的相對偏差，如表3所列。由表2可見，PC方法與MC方法得到的遠(yuǎn)場相對功率密度均值比較接近。由表3可見，PC方法與MC方法的相對偏差基本上隨著PC展開階數(shù)與樣本數(shù)量的增加而減小，當(dāng)階數(shù)為3且樣本量為50組、階數(shù)為4且樣本量不小于30組及階數(shù)為5時，相對偏差均小于5%，說明此時建立的代理模型具有足夠高的計算精度。

表2遠(yuǎn)場相對功率密度均值Tab.2Far-field means of relative powerdensity by MC and PC methods

由于所需生成的樣本量越少，數(shù)值模擬的計算開銷越小，所以采用20組樣本和5階PC方法構(gòu)建的代理模型，將10 000組樣本代入代理模型，計算得到樣本數(shù)與相對功率密度之間的柱狀圖，如圖3所示。其中，圖3(a)為MC方法的結(jié)果，圖3(b)為PC方法的結(jié)果。根據(jù)圖3分布，估計得到遠(yuǎn)場相對功率密度的概率密度函數(shù)(PDF)，如圖4所示。由圖4可得到MC方法與PC方法計算的遠(yuǎn)場相對功率密度的不確定度分別為18.7%與18.6%，二者依然吻合較好。

(a)Relative power density by MC method

(b)Relative power density by 5-order PC method

圖4遠(yuǎn)場相對功率密度的概率密度函數(shù)Fig.4Probability density function of far-field relative power density

表3遠(yuǎn)場相對功率密度的不確定度及其與MC方法的相對偏差Tab.3Uncertainties of far-field relative power density and their relative deviations with MC method

4結(jié)論

本文針對現(xiàn)實(shí)條件下激光器與大氣環(huán)境非絕對穩(wěn)定的問題，建立了一種激光大氣傳輸湍流效應(yīng)數(shù)值模擬的不確定度量化方法。以出口光束質(zhì)量與大氣相干長度為不確定性來源，利用少量樣本構(gòu)建了多項(xiàng)式混沌代理模型，實(shí)現(xiàn)了遠(yuǎn)場相對功率密度的不確定度量化，在不損失精度的前提下大大提升了計算效率。