概率盒全局靈敏度和活躍子空間跨層降維

胡政文，張保強，鄧振鴻

廈門大學 航空航天學院，廈門 361000

由于制造、測量、計算及模型本身中存在的各種誤差，航空航天系統參數化建模和多學科優化需要考慮不確定性的影響[1-4]。在不確定性分析時，不確定性參數按照類別大致可分為隨機不確定性參數和認知不確定性參數[5]。不同類型不確定性參數混合作用與單種不確定性參數作用導致的模型輸出不確定性存在明顯的差異，因此模型分析過程中須同時考慮隨機和認知不確定性的影響[6-8]。Williamson和Downs[9]在累積分布函數的基礎上引入概率盒來處理隨機和認知不確定性混合量化問題。概率盒可以讓隨機和認知不確定性在模型中分別傳播[10]，并保留2種不確定性各自的特性，因此在混合不確定性量化分析中得到了廣泛應用[11-13]。

同時由于涉及多層次和多學科，航空航天系統參數的維數通常較多[14]。學科間耦合作用下系統通常呈現高度非線性，難以直接判斷輸入和輸出之間的關系，高維參數給系統研究分析帶來了挑戰[15]。因此系統建模前通常須考慮減少輸入參數的維數，以往研究中主要有兩種思路處理維數簡化問題：一種是參數選擇，這種思路認為在所有的輸入參數中只有少數真正與模型輸出相關，其余參數對模型輸出不確定性影響較小可以忽略，參數選擇通常采用靈敏度分析的方法[16]；第二種是降維方法，與參數選擇相比，降維方法假設所有輸入參數都影響模型的輸出結果，但只在參數的某些特定組合下才能最大程度體現出來，而降維方法的目標就是構造這類組合[17]。

靈敏度分析主要研究輸入不確定性對輸出不確定性的貢獻大小。隨著模型復雜性的增加，輸入參數的靈敏度信息逐漸成為建模時需要考慮的因素之一[18-20]。靈敏度分析方法主要分為局部靈敏度分析和全局靈敏度分析，全局靈敏度分析方法主要有Sobol’指標、矩靈敏度指標等[21-22]。傳統的靈敏度分析方法大都基于概率框架，然而實際中模型通常包含認知不確定，須同時考慮區間的靈敏度分析[21, 23]。Oberguggenberger等[24]闡述了概率和非概率框架下的靈敏度分析方法。近年來，隨著概率盒在混合不確定性量化中的應用，概率盒下的全局靈敏度研究越來越多。Ferson和Troy Tucker[25]研究了在混合不確定性概率邊界分析時通過縮減輸入參數的不確定性進行靈敏度分析。Song等[26]提出一種擴展蒙特卡羅模擬方法，計算概率盒下各輸入的基于方差的全局靈敏度指標。Bi等[27]使用Bhattacharyya距離定量描述概率盒下的不確定性參數靈敏度。不確定性縮減法分析概率盒全局靈敏度具有易于理解、工程中實現較為容易的優勢，逐漸得到研究人員的關注[28]。

另一方面，降維方法在航空發動機、衛星、汽車、芯片、圖像處理等許多領域中得到了廣泛應用[14]。Fodor[29]回顧了常用的降維方法，包括主成分分析、因子分析、隨機投影法、連續拓撲映射、回歸等。主成分分析是普遍使用的線性降維方法之一[30-31]，該方法基于參數的協方差矩陣，可以有效地降低均方誤差，也被稱為本征正交分解、Karhunen-Loève (KL) 變換、霍特林變換和經驗正交函數法等[32-33]。Constantine等[34]通過構造輸出的梯度協方差矩陣來判斷輸入參數空間中的主要方向，提出了活躍子空間降維方法。活躍子空間法推廣了傳統的主成分分析思想，對處理高維復雜系統具有一定的優勢[35]。Lukaczyk等[36]將活躍子空間法應用于ONERA-M6跨聲速機翼的形狀設計優化中，在設計空間區域內找到了最佳設計變量。Jefferson等[37]將活躍子空間法應用在綜合水文模型的靈敏度分析。Hu等[38]使用活躍子空間法對衛星系統進行多學科優化設計。

以往處理高維參數的簡化問題的研究中，一般使用靈敏度分析或降維方法兩者之一，對于高維混合不確定性系統單種方法降維效果不明顯。本文將靈敏度分析和降維方法結合在一個框架中，提出一種結合概率盒全局靈敏度和活躍子空間的跨層降維方法。案例分析中將跨層降維方法應用到NASA多學科不確定性量化挑戰問題中，驗證提出的跨層降維方法的有效性。

1 跨層降維方法基本理論

1.1 概率盒全局靈敏度分析

概率盒定義為一組包含變量所有可能累積分布函數(Cumulative Distribution Function, CDF)曲線的邊界，可以有效表示混合不確定性影響下的參數特征，圖1為典型概率盒示意圖。概率盒不確定性分析提供了一種較為直觀的方法來描述受到隨機和認知不確定性綜合作用的輸出特性。

圖1 典型概率盒示意圖Fig.1 Schematic of typical probability box

不確定性縮減法可以用于概率盒框架下不確定性參數的全局靈敏度評估。輸入參數的靈敏度取決于參數自身不確定性的大小以及對模型輸出不確定性的影響，概率盒全局靈敏度可以表示為[25]

(1)

式中：B為原始基準輸入參數；T為縮減不確定性后的輸入參數；unc(·)表示輸出不確定性的度量;S為全局靈敏度指標。若unc(·)為標量，則靈敏度指標S也是標量，由此可以對輸入參數進行靈敏度排序。輸入參數不確定性縮減后，式(1)中分子unc(T)的值會下降。與基于方差的Sobol’指標不同，由于不確定性的減少，對于所有輸入參數，unc(T)均小于unc(B)，因此靈敏度指標S范圍是[0,1]。對各輸入參數依次應用不確定性縮減法，由式(1)可以得到各輸入參數的S值，將所有的S按照大小進行排序即得到對應參數的靈敏度排序。

輸入參數不確定性的縮減有多種方法，不同縮減方式得到的不確定性度量unc(·)的結果不同。通常有以下2種方式：① 使用特定的概率分布替代不確定性輸入參數，特定的概率分布法只消除參數的認知不確定性，不影響其隨機不確定性；② 使用固定值替代不確定性輸入參數，固定值法同時消除參數的隨機不確定性和認知不確定性。

不確定性度量unc(·)也有多種不同的定義方法(如方差、矩等)，不同的定義方式可以用于解決不同的靈敏度分析問題[25]。概率盒是參數所有可能累計分布函數的包絡線，沿著概率盒的上界和下界進行積分可以得到概率盒的面積。本文結合概率盒的特征，采用概率盒的面積作為輸出的不確定性度量unc(·)，此時不確定性縮減法的流程見圖2，式(1)轉化為

圖2 概率盒靈敏度分析的不確定性縮減法流程圖Fig.2 Flow chart of pinching method for probability box sensitivity analysis

(2)

式中：A(·)表示概率盒面積。

不確定性縮減法主要比較縮減輸入參數不確定性后輸出概率盒的面積占比下降值，進而判斷參數的靈敏度信息。

1.2 活躍子空間法

假設原始輸入參數的維數為n，表示為x=[x1,x2,…,xn]T，降維方法的目的是在滿足損失數據信息最少的條件下，通過某種方法將原始的n維參數x用m維的s=[s1,s2,…,sm]T(m

(3)

對于參數x，其標準參數空間可以視為邊界為1的超立方體。非標準的參數空間可以通過正則化轉化為標準參數空間。假設在x的參數空間中隨機抽樣，樣本數為N，則輸入參數x可以表示為矩陣X={xi,j:1≤i≤n,1≤j≤N}。對于任意的xi,j，可正則化為

(4)

目標函數f梯度的協方差矩陣定義為

(5)

(6)

C=WΛWT

(7)

式中：W為n×n的特征向量矩陣；Λ=diag(λ1,λ2, …,λn)是非負特征值組成的對角矩陣。特征值大小反映了目標函數f在x參數空間中沿著對應特征向量微小擾動的平均變化。將特征值按照降序進行排列λ1≥λ2≥…≥λn，特征值越小表示對應的特征向量在梯度協方差矩陣C中的貢獻度越低，降維時可以適當舍棄這些向量。

活躍子空間確定基向量維數有多種方法：通過設定固定閾值λ0，只保留大于λ0的特征值及對應的特征向量；定義κ為前m項特征值占所有特征值的累加比例：

(8)

(9)

(10)

式中：U為活躍子空間的基向量，包含W的前m列；Λ1為對應的特征值對角矩陣。確定活躍子空間的基向量后，原輸入的參數空間可以投影到該活躍子空間：

s=UTx

(11)

此時目標函數f(x)可以近似表示為活躍子空間內降維變量s的函數：

f(x)≈g(UTx)=g(s)

(12)

g(s)的定義域為γ={s=UTx,x∈Rn}∈Rm。在降維變量s的基礎上構建低維代理模型：

g(s)≈g*(s)≡R(s;g1,g2,…,gM)

(13)

式中：M為代理模型訓練樣本數；g*(s)≡R為活躍子空間中定義的代理模型，可以降低計算模型的復雜性和減少計算時間和成本。原函數和代理模型的誤差可以表示為[34]

(14)

系數C1取決于x和ρ(x)。由式(14)可以看出，舍棄的特征值之和越小，即累加比例κ的值越大，原函數和新函數誤差也越小，即活躍子空間降維法所產生的誤差也就越小。若前m項特征值累加比例κ>90%，則可以選擇前m項特征向量組成基向量用于降低輸入參數空間的維度。

1.3 跨層降維方法

在不確定性縮減法和活躍子空間的基礎上，提出結合概率盒全局靈敏度和活躍子空間的跨層降維方法。

圖3是跨層降維方法的流程圖，基本步驟如下：

圖3 跨層降維方法流程圖Fig.3 Flow chart of cross-layer method for dimension reduction

步驟1整理系統的不確定性參數，將其按照隨機、認知和混合不確定性進行分類。

步驟2使用概率盒法對初始模型進行不確定性分析，得到模型輸出的初始概率盒，計算概率盒的面積。

步驟3縮減待分析的輸入參數的不確定性，對模型進行分析得到縮減不確定性后輸出概率盒的面積。

步驟4比較縮減參數不確定性前后輸出概率盒的面積改變量，計算參數的靈敏度指標S值。

步驟5選取另一輸入參數，重復步驟3和步驟4，直至完成所有不確定性參數的靈敏度分析。

步驟6對各輸入參數的S值按照降序排列，將S值較小的輸入參數固定。

步驟7對保留不確定性的輸入參數進行抽樣，取樣本數N，計算模型輸入并得到梯度協方差矩陣。

步驟8對梯度協方差矩陣進行特征分解，按照特征值大小對特征值進行排序，確定保留的活躍子空間的維數和基向量。

步驟9對輸入參數使用活躍子空間的基向量進行降維。

2 NASA多學科不確定性量化挑戰問題

為了驗證提出的概率盒全局靈敏度和活躍子空間的跨層降維方法，本文以NASA Langley研究中心提出的多學科不確定性量化挑戰問題為案例[40]。

2.1 挑戰問題描述

挑戰問題的數學模型描述了地面基站無線電波操控遙控飛機的通用傳輸動力學特性，如圖4所示[40]。要求對不確定性系統模型進行綜合評估，解決不確定性量化、全局靈敏度分析、極限狀態分析和穩健優化設計中的一系列關鍵問題。

圖4 NASA多學科不確定性量化挑戰問題模型示意圖[40]Fig.4 Model schematic of NASA multidisciplinary uncertainty quantification challenge problem[40]

模型中的不確定性參數按照類型劃分為隨機、認知以及混合不確定性參數3類。Ψ表示所研究的系統模型，p表示Ψ中不確定性輸入參數(p∈R21)，d表示設計參數(d∈R14)，g表示Ψ性能約束(g∈R8)，g取決于p和d。不確定性輸入參數p、設計參數d和約束g之間的關系可以表示為

x=h(p)

(15)

g=f(x,d)

(16)

式中：x是參數為p的中間函數(x∈R5)；g和x假定為連續可微函數。式(15)展開為

(17)

x的每個分量都表示系統一個學科的輸出，x1～x5均為標量。若系統Ψ滿足對于g的所有分量均滿足gi≤0,i=1,2,…,8，則可以認為該系統是可靠的，符合設計要求。對于確定設計參數d，滿足約束向量g≤0的輸入參數p的集合稱為安全域，反之稱為失效域。g的分量有一個不滿足約束條件，對應的參數p就落在失效域內。圖5為挑戰問題系統模型Ψ中各參數變量的層次關系示意圖。

圖5 NASA挑戰問題模型參數層次關系Fig.5 Hierarchical relationship of model parameters in NASA challenge problem

2.2 全局靈敏度分析

由于系統模型Ψ的輸入參數中包含隨機、認知以及混合不確定性，首先對Ψ進行概率盒不確定性分析，研究系統的不確定性傳遞和量化的相關問題，確定不確定性輸入參數p對學科輸出x的全局靈敏度，以便后續活躍子空間降維過程中固定靈敏度低的參數。

輸入參數p的不確定性通過函數h傳播，得到x1、x2、x3、x4和x5的概率盒，縮減p中的認知和混合不確定性參數時，對應x概率盒的面積也會改變。

由式(17)可知，子學科Ψ1的輸出x1是p1、p2、p3、p4和p5的函數，表1給出了參數p1～p5的具體分布信息，其中p1為單峰Beta分布族組成的混合不確定性參數，p2是區間不確定性參數，p3服從[0,1]均勻分布，p4和p5為相關的二元正態分布族組成的混合不確定性參數。表中:Δ表示參數的區間范圍;E(·)、V(·)和ρ分別代表均值、方差和相關系數。隨機不確定性是系統模型的固有屬性，對模型輸出影響較大，為減少計算成本，只對認知和混合不確定性參數按照S值大小排序。

表1 x1對應不確定性輸入參數p1～p5

對于子學科Ψ1，對初始條件下的x1進行概率盒不確定性分析。取1 000條x1的樣本累積概率分布曲線，得到初始條件下x1的概率盒，如圖6 所示，此時x1概率盒的面積為0.229。

圖6 初始條件下的x1概率盒Fig.6 Probability box of x1 under initial conditions

使用不確定性縮減法分別縮減參數p1、p2、p4、p5的不確定性。由于p1為混合不確定性參數，為計算p1對x1的概率盒全局靈敏度，不確定性縮減法的固定值取為p1的期望值0.7。p2為認知不確定性參數，對應的縮減值取為區間中值0.5。p4～p5服從二維正態分布，相關系數ρ的絕對值在0～1范圍內，在計算p4的靈敏度時，取p4為期望值0，p5仍然服從正態分布，得到p4的靈敏度，同理計算p5的概率盒全局靈敏度。綜合考慮計算精度以及效率成本，結合x1概率盒邊界情況，每個參數均取1 000條累積分布函數曲線，概率盒邊界趨于穩定，進而得到縮減不確定性后各參數對應的x1概率盒，如圖7所示。運用式(2)計算得到各參數的靈敏度指標S值，計算結果列于表2之中。由圖7和表2可以看出，子學科Ψ1的參數靈敏度排序依次是p1>p5>p4>p2。由靈敏度指標S值的計算結果，p2和p4對x1概率盒的靈敏度均小于0.1，可以將p2和p4視為固定的參數。

表2 x1概率盒的面積及對應輸入參數靈敏度排序

圖7 縮減各參數不確定性對應的x1概率盒Fig.7 Probability boxes of x1 after pinching corresponding uncertain parameters

其余子學科x2、x3、x4和x5對應的16個不確定性輸入參數p6～p21的具體分布信息如表3所示。

表3 x2～x5對應不確定性輸入參數p6～p21Table 3 Description of uncertain input parametersp6-p21 corresponding to x2-x5

表4 x2概率盒的面積及對應輸入參數靈敏度排序

表5 x3概率盒的面積及對應輸入參數靈敏度排序

表6 x4概率盒的面積及對應輸入參數靈敏度排序

圖8 縮減參數p21不確定性對應的x5概率盒Fig.8 Probability boxes of x5 after pinching p21

綜合概率盒靈敏度分析的計算結果，篩選出p1、p5、p6、p7、p12、p16、p17、p18、p21共9個影響程度較大的認知或混合不確定性參數以及p3、p9、p11、p19共4個隨機不確定性參數。

2.3 活躍子空間降維

根據全局靈敏度分析的結果，共有13個影響較大的不確定性參數，應用概率盒全局靈敏度和活躍子空間的跨層降維方法，采取活躍子空間進一步降低不確定性參數的維度。

為描述系統模型的可靠性，定義如下系統性能度量函數：

(18)

性能度量函數ζ描述了系統模型最容易超出約束的g分量，間接反映了系統失效的最大可能性。為構建合適的活躍子空間基向量，將ζ作為目標函數。由拉丁超立方抽樣選取300個輸入參數樣本，分別計算對應目標函數ζ的值，并根據式(6) 構造ζ的近似梯度協方差矩陣C。

對梯度協方差矩陣進行特征分解，將特征值按照從大到小排序。表7列出了C的13個特征值，圖9為特征值隨著維度的衰減曲線圖。由表7可知，矩陣C的一維對應的特征值為0.063 84，遠遠大于二維的特征值1.400×10-17，且一維特征值占所有特征值和的比例κ超過99.99%。由式(8) 和式(14)可知，對于目標函數ζ，活躍子空間維數為一，選取一維對應的特征向量作為活躍子空間的基向量。此時式(11)中基向量，U=[0.256，-0.087 4，-0.010 1，0.015 9，0.025 5，-0.082 6，-0.033 6，0.026 3，-0.019 0，-0.038 2，0.088 8，0.103，0.94 7]。

圖9 梯度協方差矩陣的特征值衰減曲線Fig.9 Curve of eigenvalues of gradient covariance matrix

表7 梯度協方差矩陣的特征值Table 7 Eigenvalues of gradient covariance matrix

為進一步驗證活躍子空間選取的基向量的準確性，可以使用統計學中Bootstrap自助法重復抽取樣本組成新的梯度協方差矩陣進行活躍子空間的基向量計算。取Bootstrap自助法重復次數為1 000，圖10和圖11分別為活躍子空間基向量各分量分布柱狀圖和區間，由圖可知活躍子空間基向量的各分量分布較為集中，可以認為求得的U是準確的。

圖10 Bootstrap自助法(1 000次)活躍子空間基向量分量分布柱狀圖Fig.10 Histograms of components of active subspace vector using Bootstrap (1 000 times)

圖11 Bootstrap自助法(1 000次)活躍子空間基向量分量分布區間Fig.11 Ranges of components of active subspace vector using Bootstrap (1 000 times)

確定活躍子空間的基向量后，由式(11)，原13維輸入的參數空間可以投影到一維的活躍子空間，圖12是活躍子空間降維變量和目標函數ζ在樣本數據的對應關系，由圖可以清晰看出目標函數隨著降維變量的變化趨勢，有利于進一步不確定性代理模型的構建及模型參數優化工作。

圖12 樣本下降維變量和目標函數的關系Fig.12 Relationship between dimension-reduction variables and objective function based on samples

3 結 論

針對航空航天系統不確定性建模過程中存在不確定性參數類型混雜并且維數眾多問題，本文以概率盒為不確定性分析方法，提出了一種結合概率盒全局靈敏度和活躍子空間的跨層降維法。主要研究如下：

1) 針對NASA多學科不確定性量化挑戰問題，使用概率盒不確定性分析方法處理x1～x5混合不確定性量化問題，使得隨機、認知和混合不確定性在模型中有效傳播。

2) 基于概率盒的計算結果，使用不確定性縮減法對p1～p21認知和混合不確定性參數進行全局靈敏度分析，得到對應參數的靈敏度指標S值，將輸入參數按照靈敏度進行排序。

3) 在全局靈敏度分析的基礎上，固定8個S值小的不確定性參數，對其余13個輸入變量使用活躍子空間降維法。根據輸出的梯度協方差矩陣特征值的大小確定活躍子空間的維數和降維基向量，將原始的輸入參數空間映射到一維的活躍子空間內，降低了系統模型輸入參數的維度。

4) 結果顯示提出的跨層降維方法能夠將21維不確定性參數維數降低為一維，證明了該方法降低混合不確定性高維復雜系統輸入參數維度的有效性，顯示出對于處理高維不確定復雜系統、構建精度符合要求且計算效率較高的模型方面的潛力和優勢，為進一步的模型參數優化奠定基礎。