基于混合蒙特卡洛/多項(xiàng)式混沌展開(kāi)方法的多參數(shù)隨機(jī)等離子體不確定性分析

劉江凡 劉曉妹 李錚 焦子涵 徐聰 席曉莉

（1.西安理工大學(xué), 西安 710100；2.中國(guó)運(yùn)載火箭技術(shù)研究室空間物理重點(diǎn)實(shí)驗(yàn)室, 北京 100076）

0 引 言

由于電磁模型的復(fù)雜性不斷增加，需要對(duì)直接影響輸出響應(yīng)量的不確定性進(jìn)行量化和分析。對(duì)于此類分析，多項(xiàng)式混沌展開(kāi)(polynomial chaos expansion, PCE)方法已經(jīng)展示出了相當(dāng)大的潛力，有望替代耗時(shí)的蒙特卡洛(Monte Carlo, MC)方法。PCE方法的有效性已經(jīng)在各種數(shù)值仿真中得到了驗(yàn)證，在這些仿真試驗(yàn)中，具有少量隨機(jī)變量輸入的隨機(jī)輸出響應(yīng)量可以通過(guò)少量多項(xiàng)式展開(kāi)項(xiàng)充分逼近[1-4]。然而，隨著隨機(jī)變量維數(shù)的增加，在相同的多項(xiàng)式總階數(shù)D下，展開(kāi)項(xiàng)的數(shù)量迅速增加。例如，D=1,2,3時(shí)，有2個(gè)隨機(jī)輸入變量的PCE模擬，分別需要3，6，10個(gè)多項(xiàng)式展開(kāi)項(xiàng)。相比之下，對(duì)于相同的總階數(shù)D，具有10個(gè)隨機(jī)變量的模擬需要11，66和286個(gè)展開(kāi)項(xiàng)。隨著隨機(jī)變量數(shù)量的快速增長(zhǎng)，多項(xiàng)式展開(kāi)項(xiàng)的數(shù)量迅速增加，會(huì)顯著降低PCE方法的效率。目前已有的比較成熟的減少PCE項(xiàng)以及所需樣本個(gè)數(shù)的方法主要有稀疏多項(xiàng)式混沌法[5-8]和稀疏節(jié)點(diǎn)法[9-11]。

為克服多參數(shù)PCE模擬所產(chǎn)生的巨大計(jì)算成本，我們采用一種混合MC/PCE方法[12-14]，該方法是一種使用已知近似隨機(jī)函數(shù)的特殊方法。由于傳統(tǒng)的MC方法通過(guò)簡(jiǎn)單地增加迭代次數(shù)或者降低隨機(jī)變量的標(biāo)準(zhǔn)差σ來(lái)提高計(jì)算結(jié)果估計(jì)的精度，σ越小，計(jì)算結(jié)果的精度越高，但是計(jì)算成本過(guò)高。因此，在本文中，使用非侵入式多項(xiàng)式混沌(nonintrusive polynomial chaos，NIPC)方法得到相對(duì)應(yīng)的近似隨機(jī)函數(shù)，并與少量MC方法聯(lián)合仿真，從而達(dá)到多次MC模擬結(jié)果的精度。同時(shí)基于回歸法的PCE方法不使用越來(lái)越高的多項(xiàng)式階數(shù)來(lái)提高精度，避免了大量的多項(xiàng)式展開(kāi)項(xiàng)數(shù)，從而達(dá)到緩解“維數(shù)災(zāi)難”問(wèn)題的目的。該混合方法不需要MC的大量樣本，同時(shí)提高了低階多項(xiàng)式混沌方法的計(jì)算精度，降低了計(jì)算復(fù)雜度，節(jié)約了時(shí)間和成本。

在本文中，將混合MC/PCE方法應(yīng)用于一維電磁波傳播問(wèn)題中，即平面波入射非均勻分層等離子體平板，每一層的等離子體電子密度相互獨(dú)立且服從高斯分布；并利用該模型來(lái)驗(yàn)證所提出的混合MC/PCE算法的準(zhǔn)確性和計(jì)算效率。

1 PCE方法

1.1 多項(xiàng)式混沌基礎(chǔ)

多項(xiàng)式混沌法起源于Wiener各項(xiàng)同性混沌理論中的隨機(jī)變量譜展開(kāi)[15]。近年來(lái)，PCE方法被廣泛應(yīng)用于不確定性分析問(wèn)題中。多項(xiàng)式混沌方法是一種基于不確定性譜表示的隨機(jī)方法。譜表示不確定性的一個(gè)重要方面是可以將隨機(jī)函數(shù)分成確定成分和隨機(jī)成分。例如，對(duì)于隨機(jī)媒質(zhì)電磁散射問(wèn)題中的任何隨機(jī)函數(shù)(即α?)，如雷達(dá)散射截面、透射系數(shù)等，混沌多項(xiàng)式理論下輸出響應(yīng)可表示為p階截?cái)嗟幕煦缍囗?xiàng)式模型[14]：

式中：αi(x)為待求的混沌多項(xiàng)式展開(kāi)系數(shù)；ψi(ξ)為由隨機(jī)變量的分布函數(shù)ξ確定的正交多項(xiàng)式基函數(shù)的第i項(xiàng)展開(kāi)項(xiàng)。假設(shè)α?是確定性自變量向量x和n個(gè)獨(dú)立隨機(jī)變量的向量ξ=[ξ1,···,ξi,···,ξn]的函數(shù)，ξ的分布函數(shù)類型由文獻(xiàn)[16]給出。

本文隨機(jī)變量被定義為高斯分布，則對(duì)應(yīng)的多項(xiàng)式是厄米特(Hermite)多項(xiàng)式，展開(kāi)項(xiàng)數(shù)P與Hermite多項(xiàng)式最高階d的關(guān)系為

式中，n為隨機(jī)變量的數(shù)量?？梢钥闯?，只有1個(gè)隨機(jī)變量時(shí)，展開(kāi)項(xiàng)數(shù)P和多項(xiàng)式最高階d是相等的，即P=d。表1給出了只有1個(gè)隨機(jī)變量且多項(xiàng)式最高階d=4的PCE基函數(shù)。而有兩個(gè)及兩個(gè)以上數(shù)量的隨機(jī)變量時(shí)，展開(kāi)項(xiàng)數(shù)P會(huì)成倍增加。例如，如果有3個(gè)隨機(jī)變量，多項(xiàng)式最高階d={1,2,3}時(shí)，展開(kāi)項(xiàng)數(shù)為{4,10,20}。表2為總階數(shù)d=2和3個(gè)隨機(jī)參數(shù)的PCE基函數(shù)。

表1 總階數(shù)d=4和1個(gè)隨機(jī)參數(shù)的PCE基函數(shù)Tab.1 PCE basis function with total order d=4 and one random parameter

表2 總階數(shù)d=2和3個(gè)隨機(jī)參數(shù)的PCE基函數(shù)Tab.2 PCE basis function with total order d=2 and three random parameters

1.2 NIPC方法

混沌多項(xiàng)式系數(shù)求取是采用多項(xiàng)式展開(kāi)方法進(jìn)行不確定性分析的關(guān)鍵，決定著α?的精度。通常，可通過(guò)侵入式法(intrusive method)和非侵入式法(nonintrusive method)兩種方法來(lái)確定PCE系數(shù)[17]。計(jì)算PCE系數(shù)的過(guò)程中，侵入式方法需要對(duì)原始模型進(jìn)行改進(jìn)和調(diào)整，而非侵入式方法則不需要改變?cè)寄Ｐ?，因此NIPC方法廣受關(guān)注。為了提高計(jì)算效率，本文采用點(diǎn)配置NIPC方法。

點(diǎn)配置NIPC方法首先用方程中的PCE替換隨機(jī)響應(yīng)或隨機(jī)函數(shù)，然后在隨機(jī)空間中選擇一組ξ=[ξ1,···,ξn]的樣本，最后建立一個(gè)有N個(gè)方程組成的線性方程組?？赏ㄟ^(guò)試驗(yàn)或者數(shù)值模擬獲取對(duì)應(yīng)的α?(x,ξ)，由式(1)可得式(3)所示的方程組，再通過(guò)最小二乘法求解系數(shù)αi。

求出混沌多項(xiàng)式系數(shù)αi，就可以直接計(jì)算出輸出響應(yīng)的隨機(jī)概率特性。基于NIPC展開(kāi)各項(xiàng)系數(shù)，即可計(jì)算隨機(jī)函數(shù)α?的統(tǒng)計(jì)特性，可得：

隨機(jī)函數(shù)的均值：

標(biāo)準(zhǔn)差：

出于計(jì)算量的考慮，目前仍局限于選取其中部分樣本進(jìn)行回歸，通常選取概率空間出現(xiàn)頻率較大的樣本點(diǎn)。由于樣本點(diǎn)數(shù)目少，可以全部用來(lái)估算混沌多項(xiàng)式系數(shù)，同時(shí)保證不確定性的精度。但是，由于回歸法中涉及矩陣求逆等運(yùn)算，針對(duì)高維問(wèn)題可能較為繁瑣耗時(shí)。所以要將回歸法應(yīng)用于高維問(wèn)題，需要預(yù)先篩選多項(xiàng)式。

2 混合MC/PCE方法

在許多工業(yè)問(wèn)題中經(jīng)常遇到多參數(shù)不確定性分析的情況，一種MC與PCE的混合算法被應(yīng)用于多參數(shù)不確定性分析。

2.1 基于PCE的混合方法理論

對(duì)于隨機(jī)變量X，其具有精確均值E[X]、精確標(biāo)準(zhǔn)差std[X]、精確方差var[X]、概率密度函數(shù)P[ξ]，經(jīng)過(guò)N次迭代，且N足夠大，則MC均值可被定義為

將隨機(jī)變量X的低階PCE近似為隨機(jī)函數(shù)G(ξ)：

將μ[X-λG]定義為

我們可以通過(guò)設(shè)置λ的值使得var[X-λG]最小，即：

對(duì)應(yīng)的最小值為

協(xié)方差cov(X,G)定義為

混合MC/PCE方法有效地將均值估計(jì)量μ[X]轉(zhuǎn)化為兩個(gè)積分，一個(gè)用低階NIPC求解，另一個(gè)用MC模擬求解，具有比傳統(tǒng)MC算法更快的收斂速度。方差估計(jì)量σ2[X]用相同的方法得到：

2.2 混合MC/PCE方法模型

混合MC/PCE方法模型構(gòu)建出直接描述多參數(shù)不確定性和輸出響應(yīng)不確定性之間的函數(shù)關(guān)系。該模型的輸入由兩部分組成：一是由電磁仿真算法模擬為MC提供大量樣本點(diǎn)從而獲得的均值和標(biāo)準(zhǔn)差；二是由電磁仿真算法模擬為NIPC模型提供輸出響應(yīng)量，從而獲得多項(xiàng)式混沌系數(shù)進(jìn)而重構(gòu)混沌多項(xiàng)式展開(kāi)表達(dá)式并得到NIPC多項(xiàng)式的均值和標(biāo)準(zhǔn)差。具體構(gòu)建流程如圖1所示。

圖1 混合MC/PCE方法模型構(gòu)建流程Fig.1 Hybrid MC/PCE method of model building process

具體實(shí)施步驟如下：

1) 確定算例中不確定性參數(shù)的分布類型，本次算例中假設(shè)隨機(jī)參數(shù)均服從高斯分布，故MC方法的抽樣方式采取隨機(jī)抽樣，PCE方法的多項(xiàng)式基函數(shù)選擇Hermite多項(xiàng)式。

2) 根據(jù)算例中不確定性參數(shù)的個(gè)數(shù)N，采用隨機(jī)抽樣方法對(duì)各個(gè)參數(shù)進(jìn)行M組抽樣，本次抽樣設(shè)置樣本點(diǎn)500組。

3) 采用電磁仿真算法對(duì)抽樣的樣本點(diǎn)進(jìn)行MC仿真試驗(yàn)，得到MC樣本以及MC的均值和標(biāo)準(zhǔn)差。

4) 根據(jù)算例中不確定性參數(shù)的個(gè)數(shù)N，確定多項(xiàng)式階數(shù)為d=1以及多項(xiàng)式展開(kāi)項(xiàng)數(shù)N+1，對(duì)各個(gè)參數(shù)進(jìn)行N+1次隨機(jī)抽樣，即輸入?yún)?shù)的N+1組樣本{ξ(1),ξ(2),···,ξ(N+1)}，并利用電磁仿真算法得到輸出響應(yīng)量α?的N+1組樣本{α?(1),α?(2),···,α?(N+1)}。

5) 對(duì)第4步得到的N+1組樣本采用NIPC方法對(duì)算例進(jìn)行模擬，利用最小二乘法回歸分析構(gòu)建NIPC模型并計(jì)算得到一組PCE系數(shù)，由此系數(shù)得到NIPC的均值和標(biāo)準(zhǔn)差。

6) 根據(jù)NIPC得到的多項(xiàng)式展開(kāi)系數(shù)重構(gòu)PCE表達(dá)式，并將第2步得到的隨機(jī)抽樣的500組樣本點(diǎn)帶入重構(gòu)的表達(dá)式中得到隨機(jī)函數(shù)G(ξ)樣本。

7) 將第3步得到的MC的樣本和均值以及第6步得到的NIPC的樣本和均值帶入式(13)得到協(xié)方差cov(X,G)。

8) 根據(jù)第7步得到的協(xié)方差計(jì)算λ。

9) 將MC仿真結(jié)果，NIPC樣本以及λ的值帶入式(9)可以得到μ[X-λG]。

10)利用式(8)可以得到輸出響應(yīng)量的均值。

11)利用式(14)得到輸出響應(yīng)量的標(biāo)準(zhǔn)差，這樣就得到了混合MC/PCE模型的統(tǒng)計(jì)特性。

所構(gòu)建的混合MC/PCE方法模型代表的是不確定參數(shù)與輸出響應(yīng)均值和標(biāo)準(zhǔn)差之間的函數(shù)關(guān)系，通過(guò)給定不確定性參數(shù)均值和標(biāo)準(zhǔn)差，便可通過(guò)已構(gòu)建的混合模型快速計(jì)算出輸出響應(yīng)的均值和標(biāo)準(zhǔn)差。

3 仿真結(jié)果

算例1平面波斜入射在非均勻等離子體平板上，平面波頻率為30 GHz，電磁波入射等離子體平板的入射角度為θ0=30°。仿真模型如圖2所示。等離子體區(qū)域被平均劃分為3層，非均勻各向同性等離子體厚度為d1=d2=d3=50 mm，每層等離子體板的電子密度相互獨(dú)立且呈高斯分布。本算例采用HMM算法[22]來(lái)計(jì)算等離子體對(duì)平面波的透射系數(shù)。每層等離子體均值分別為μ[ne1]=1×1018m-3，μ[ne2]=1×1018m-3，μ[ne3]=1×1018m-3，標(biāo)準(zhǔn)差分別為σ[ne1]=0.05×μ[ne1]， σ[ne2]=0.05×μ[ne2]，σ[ne3]=0.05×μ[ne3]，碰撞頻率均為v=1 GHz。

圖2 平面波入射3層非均勻等離子體平板仿真模型Fig.2 The simulation model for plane wave propagation in non-uniformly 3-layered plasma

圖3給出了電磁波斜入射3層非均勻等離子體平板的透射系數(shù)統(tǒng)計(jì)特性?？梢钥吹剑疚幕旌螹C/PCE方法計(jì)算的結(jié)果精度和10 000次MC方法可以很好地匹配，尤其是均值計(jì)算結(jié)果。算例采用一階NIPC方法和500次MC進(jìn)行了聯(lián)合，同時(shí)分別將一階NIPC方法、兩階NIPC方法和500次MC的計(jì)算結(jié)果與本文方法進(jìn)行了比較。從圖3(b)可以明顯看到，混合方法相比于一階NIPC方法和500次MC有更高的精度，與10 000次MC方法相比，計(jì)算量明顯減少。兩階NIPC方法比一階NIPC方法精度更高，這是因?yàn)镹IPC方法可以通過(guò)增加多項(xiàng)式階數(shù)來(lái)提高計(jì)算精度，但是高階NIPC方法需要更多的展開(kāi)項(xiàng)數(shù)，增加了算法的復(fù)雜度和計(jì)算成本。

圖3 透射系數(shù)幅度的均值和標(biāo)準(zhǔn)差（算例1）Fig.3 The mean and standard deviation of transmission coefficient amplitude (sample 1)

算例2平面波斜入射在非均勻等離子體板上，仿真模型如圖4所示。模擬區(qū)域被劃分為10層，非均勻各向同性等離子體厚度均為d=6 mm，每層等離子體板的電子密度獨(dú)立且呈高斯分布。10層等離子電子密度值如表3所示。標(biāo)準(zhǔn)差均為均值的十分之一，碰撞頻率為v=5 GHz。

圖4 平面波入射10層非均勻等離子體平板仿真模型Fig.4 The simulation model for plane wave propagation in non-uniformly 10-layered plasma

表3 等離子體的電子密度值Tab.3 The electron density value of the plasma

圖5分別給出了電磁波通過(guò)10層隨機(jī)等離子體平板透射系數(shù)的幅度和相位的統(tǒng)計(jì)特性。算例采用一階NIPC方法和500次MC進(jìn)行了聯(lián)合，同時(shí)分別將一階NIPC方法，500次和10 000次MC的計(jì)算結(jié)果與本文混合MC/NIPC方法進(jìn)行了比較?？梢钥闯?，4種方法得到的透射系數(shù)幅度和相位的均值結(jié)果較一致；而對(duì)于透射系數(shù)幅度和相位的標(biāo)準(zhǔn)差來(lái)說(shuō)，相比于一階NIPC方法和500次的MC，本文中混合MC/PCE方法具有明顯的優(yōu)勢(shì)，計(jì)算結(jié)果與10 000次MC方法的結(jié)果更加匹配。

圖5 透射系數(shù)幅度和相位的均值和標(biāo)準(zhǔn)差（算例2）Fig.5 The mean and standard deviation of the transmission coefficient amplitude and phase (sample 2)

表4給出了該算例中NIPC，MC以及混合方法的模擬次數(shù)?？梢钥闯觯阂浑ANIPC方法仿真次數(shù)最少只有11次，這是因?yàn)橹恍枰?jì)算11個(gè)多項(xiàng)式展開(kāi)系數(shù)；500次的MC仿真相對(duì)于10 000次MC精度較差；本文所采用的混合方法具有跟10 000次MC方法相比擬的計(jì)算結(jié)果，但是仿真次數(shù)只需要500次，避免了大量的運(yùn)算。

表4 不同算法仿真次數(shù)的比較Tab.4 Comparison of the number of simulations for different algorithms

4 結(jié) 論

本文采用一種混合MC/PCE方法，進(jìn)行了多參數(shù)隨機(jī)等離子體不確定性分析。采用該混合方法進(jìn)行分層等離子體板多參數(shù)電磁散射特性不確定性分析，并將計(jì)算結(jié)果與MC結(jié)果作對(duì)比，驗(yàn)證了該算法的正確性，避免了隨機(jī)變量增多時(shí)展開(kāi)項(xiàng)數(shù)隨之增多的問(wèn)題，同時(shí)也避免了回歸法在高維問(wèn)題時(shí)由于矩陣求逆運(yùn)算使得多參數(shù)不確定性分析出現(xiàn)繁瑣耗時(shí)的問(wèn)題。對(duì)于本文提出的混合MC/PCE方法，后續(xù)還可以將其應(yīng)用于更復(fù)雜的隨機(jī)等離子體電磁散射特性研究中，如等離子鞘套電磁散射特性。