基于混合蒙特卡洛/多項式混沌展開方法的多參數隨機等離子體不確定性分析

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

（1.西安理工大學, 西安 710100；2.中國運載火箭技術研究室空間物理重點實驗室, 北京 100076）

0 引 言

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

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

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

1 PCE方法

1.1 多項式混沌基礎

多項式混沌法起源于Wiener各項同性混沌理論中的隨機變量譜展開[15]。近年來，PCE方法被廣泛應用于不確定性分析問題中。多項式混沌方法是一種基于不確定性譜表示的隨機方法。譜表示不確定性的一個重要方面是可以將隨機函數分成確定成分和隨機成分。例如，對于隨機媒質電磁散射問題中的任何隨機函數(即α?)，如雷達散射截面、透射系數等，混沌多項式理論下輸出響應可表示為p階截斷的混沌多項式模型[14]：

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

本文隨機變量被定義為高斯分布，則對應的多項式是厄米特(Hermite)多項式，展開項數P與Hermite多項式最高階d的關系為

式中，n為隨機變量的數量。可以看出，只有1個隨機變量時，展開項數P和多項式最高階d是相等的，即P=d。表1給出了只有1個隨機變量且多項式最高階d=4的PCE基函數。而有兩個及兩個以上數量的隨機變量時，展開項數P會成倍增加。例如，如果有3個隨機變量，多項式最高階d={1,2,3}時，展開項數為{4,10,20}。表2為總階數d=2和3個隨機參數的PCE基函數。

表1 總階數d=4和1個隨機參數的PCE基函數Tab.1 PCE basis function with total order d=4 and one random parameter

表2 總階數d=2和3個隨機參數的PCE基函數Tab.2 PCE basis function with total order d=2 and three random parameters

1.2 NIPC方法

混沌多項式系數求取是采用多項式展開方法進行不確定性分析的關鍵，決定著α?的精度。通常，可通過侵入式法(intrusive method)和非侵入式法(nonintrusive method)兩種方法來確定PCE系數[17]。計算PCE系數的過程中，侵入式方法需要對原始模型進行改進和調整，而非侵入式方法則不需要改變原始模型，因此NIPC方法廣受關注。為了提高計算效率，本文采用點配置NIPC方法。

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

求出混沌多項式系數αi，就可以直接計算出輸出響應的隨機概率特性?；贜IPC展開各項系數，即可計算隨機函數α?的統計特性，可得：

隨機函數的均值：

標準差：

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

2 混合MC/PCE方法

在許多工業問題中經常遇到多參數不確定性分析的情況，一種MC與PCE的混合算法被應用于多參數不確定性分析。

2.1 基于PCE的混合方法理論

對于隨機變量X，其具有精確均值E[X]、精確標準差std[X]、精確方差var[X]、概率密度函數P[ξ]，經過N次迭代，且N足夠大，則MC均值可被定義為

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

將μ[X-λG]定義為

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

對應的最小值為

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

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

2.2 混合MC/PCE方法模型

混合MC/PCE方法模型構建出直接描述多參數不確定性和輸出響應不確定性之間的函數關系。該模型的輸入由兩部分組成：一是由電磁仿真算法模擬為MC提供大量樣本點從而獲得的均值和標準差；二是由電磁仿真算法模擬為NIPC模型提供輸出響應量，從而獲得多項式混沌系數進而重構混沌多項式展開表達式并得到NIPC多項式的均值和標準差。具體構建流程如圖1所示。

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

具體實施步驟如下：

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

2) 根據算例中不確定性參數的個數N，采用隨機抽樣方法對各個參數進行M組抽樣，本次抽樣設置樣本點500組。

3) 采用電磁仿真算法對抽樣的樣本點進行MC仿真試驗，得到MC樣本以及MC的均值和標準差。

4) 根據算例中不確定性參數的個數N，確定多項式階數為d=1以及多項式展開項數N+1，對各個參數進行N+1次隨機抽樣，即輸入參數的N+1組樣本{ξ(1),ξ(2),···,ξ(N+1)}，并利用電磁仿真算法得到輸出響應量α?的N+1組樣本{α?(1),α?(2),···,α?(N+1)}。

5) 對第4步得到的N+1組樣本采用NIPC方法對算例進行模擬，利用最小二乘法回歸分析構建NIPC模型并計算得到一組PCE系數，由此系數得到NIPC的均值和標準差。

6) 根據NIPC得到的多項式展開系數重構PCE表達式，并將第2步得到的隨機抽樣的500組樣本點帶入重構的表達式中得到隨機函數G(ξ)樣本。

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

8) 根據第7步得到的協方差計算λ。

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

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

11)利用式(14)得到輸出響應量的標準差，這樣就得到了混合MC/PCE模型的統計特性。

所構建的混合MC/PCE方法模型代表的是不確定參數與輸出響應均值和標準差之間的函數關系，通過給定不確定性參數均值和標準差，便可通過已構建的混合模型快速計算出輸出響應的均值和標準差。

3 仿真結果

算例1平面波斜入射在非均勻等離子體平板上，平面波頻率為30 GHz，電磁波入射等離子體平板的入射角度為θ0=30°。仿真模型如圖2所示。等離子體區域被平均劃分為3層，非均勻各向同性等離子體厚度為d1=d2=d3=50 mm，每層等離子體板的電子密度相互獨立且呈高斯分布。本算例采用HMM算法[22]來計算等離子體對平面波的透射系數。每層等離子體均值分別為μ[ne1]=1×1018m-3，μ[ne2]=1×1018m-3，μ[ne3]=1×1018m-3，標準差分別為σ[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層非均勻等離子體平板的透射系數統計特性?？梢钥吹?，本文混合MC/PCE方法計算的結果精度和10 000次MC方法可以很好地匹配，尤其是均值計算結果。算例采用一階NIPC方法和500次MC進行了聯合，同時分別將一階NIPC方法、兩階NIPC方法和500次MC的計算結果與本文方法進行了比較。從圖3(b)可以明顯看到，混合方法相比于一階NIPC方法和500次MC有更高的精度，與10 000次MC方法相比，計算量明顯減少。兩階NIPC方法比一階NIPC方法精度更高，這是因為NIPC方法可以通過增加多項式階數來提高計算精度，但是高階NIPC方法需要更多的展開項數，增加了算法的復雜度和計算成本。

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

算例2平面波斜入射在非均勻等離子體板上，仿真模型如圖4所示。模擬區域被劃分為10層，非均勻各向同性等離子體厚度均為d=6 mm，每層等離子體板的電子密度獨立且呈高斯分布。10層等離子電子密度值如表3所示。標準差均為均值的十分之一，碰撞頻率為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分別給出了電磁波通過10層隨機等離子體平板透射系數的幅度和相位的統計特性。算例采用一階NIPC方法和500次MC進行了聯合，同時分別將一階NIPC方法，500次和10 000次MC的計算結果與本文混合MC/NIPC方法進行了比較?？梢钥闯?，4種方法得到的透射系數幅度和相位的均值結果較一致；而對于透射系數幅度和相位的標準差來說，相比于一階NIPC方法和500次的MC，本文中混合MC/PCE方法具有明顯的優勢，計算結果與10 000次MC方法的結果更加匹配。

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

表4給出了該算例中NIPC，MC以及混合方法的模擬次數。可以看出：一階NIPC方法仿真次數最少只有11次，這是因為只需要計算11個多項式展開系數；500次的MC仿真相對于10 000次MC精度較差；本文所采用的混合方法具有跟10 000次MC方法相比擬的計算結果，但是仿真次數只需要500次，避免了大量的運算。

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

4 結 論

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