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

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

0 引 言

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

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

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

1 PCE方法

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

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

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

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

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

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

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

1.2 NIPC方法

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

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

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

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

標準差：

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

2 混合MC/PCE方法

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

2.1 基于PCE的混合方法理論

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

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

將μ[X-λG]定義為

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

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

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

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

2.2 混合MC/PCE方法模型

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

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

具體實施步驟如下：

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

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

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

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

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

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

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

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

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

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

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

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

3 仿真結(jié)果

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

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

算例2平面波斜入射在非均勻等離子體板上，仿真模型如圖4所示。模擬區(qū)域被劃分為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層隨機等離子體平板透射系數(shù)的幅度和相位的統(tǒng)計特性。算例采用一階NIPC方法和500次MC進行了聯(lián)合，同時分別將一階NIPC方法，500次和10 000次MC的計算結(jié)果與本文混合MC/NIPC方法進行了比較?？梢钥闯?，4種方法得到的透射系數(shù)幅度和相位的均值結(jié)果較一致；而對于透射系數(shù)幅度和相位的標準差來說，相比于一階NIPC方法和500次的MC，本文中混合MC/PCE方法具有明顯的優(yōu)勢，計算結(jié)果與10 000次MC方法的結(jié)果更加匹配。

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

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

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

4 結(jié) 論

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