謝 恩，馬義中，劉麗君，林成龍

(南京理工大學(xué) 經(jīng)濟管理學(xué)院，江蘇 南京 210094)

0 引言

全局靈敏度分析通常假設(shè)試驗數(shù)據(jù)服從正態(tài)分布，但數(shù)據(jù)污染使得實際試驗數(shù)據(jù)偏離假設(shè),伴隨著大數(shù)據(jù)時代海量數(shù)據(jù)的出現(xiàn)，數(shù)據(jù)被污染的風(fēng)險加大[1]。同樣，基于污染數(shù)據(jù)的仿真試驗可能會歪曲事物本來的面目，降低說服力，甚至可能得出謬誤的結(jié)論。因此，如何改進(jìn)統(tǒng)計方法，使得在面對數(shù)據(jù)污染時仍能夠得到穩(wěn)健的結(jié)果，成為近年來研究的熱點問題[2-3]。樣本中存在異常值被認(rèn)為是數(shù)據(jù)污染的一種，傳統(tǒng)的異常值檢驗技術(shù)[4]通過識別異常值并直接舍棄的處理方法，往往導(dǎo)致表示過程質(zhì)量特性的重要信息丟失，影響后續(xù)工作[5]。RIPLEY[6]指出在多元或高維數(shù)據(jù)情形下難以識別數(shù)據(jù)中的異常值，需要采用穩(wěn)健統(tǒng)計方法來抵抗污染值的干擾。

基于穩(wěn)健模型和穩(wěn)健統(tǒng)計技術(shù)建立模型并估計模型參數(shù)，能更好地描述數(shù)據(jù)集中大多數(shù)數(shù)據(jù)的統(tǒng)計特性，能抵御異常值對仿真試驗的影響[7]。SERFLING[8]指出樣本數(shù)據(jù)存在污染的情形時，用中位數(shù)和中位數(shù)絕對偏差(Median Absolute Deviation, MAD)代替均值和方差估計樣本位置和尺度參數(shù)可以獲得穩(wěn)健性較高的結(jié)果。韓云霞等[9]研究了基于穩(wěn)健統(tǒng)計量(Hodges-Lehmann, HL)和Shamos估計數(shù)據(jù)污染情形下最優(yōu)參數(shù)置信區(qū)間，結(jié)果表明，基于穩(wěn)健統(tǒng)計量的方法估計精度更高，同時能夠抵抗污染值的干擾。PARK等[10]研究了數(shù)據(jù)污染情形下不同的穩(wěn)健統(tǒng)計量估計位置參數(shù)和尺度參數(shù)，通過分析蒙特卡洛仿真試驗結(jié)果，給出幾種穩(wěn)健統(tǒng)計量在數(shù)據(jù)污染情形下的有效性，其中分別用HL和Shamos估計位置參數(shù)和尺度參數(shù)進(jìn)行穩(wěn)健設(shè)計的最優(yōu)解優(yōu)于其他穩(wěn)健統(tǒng)計量。劉麗君等[11-12]用中位數(shù)估計位置參數(shù)，用MAD估計尺度參數(shù)，改進(jìn)了序貫分支試驗的顯著性檢驗程序，在不增加額外試驗次數(shù)的前提下，解決了不同類型數(shù)據(jù)污染情形下的因子篩選問題。

靈敏度分析方法主要分為3類[13-14]：①旨在區(qū)分模型輸入因子是否對模型輸出有影響的因子篩選方法(screening method)；②局部靈敏度分析方法(local sensitivity analysis)，考慮圍繞特定設(shè)計點的微小擾動對模型輸出不確定性影響；③全局靈敏度分析方法(global sensitivity analysis)，考慮整個輸入空間內(nèi)的因子變化對模型輸出不確定性的影響，旨在對模型所有輸入因子基于其對模型輸出不確定效應(yīng)的貢獻(xiàn)進(jìn)行重要度排序。Sobol’指數(shù)法作為基于方差的全局靈敏度分析方法[15-16]，獨立于模型輸入和輸出，容易解釋和實現(xiàn)，適用于模型輸入因子排序和重要變量篩選，在工業(yè)工程、環(huán)境和大氣工程及化工工程等領(lǐng)域被廣泛研究和使用[17-19]。但在面對復(fù)雜問題時，Sobol’指數(shù)通過大量的模型估計也很難得到收斂的合理解，無法識別因子重要度。隨著計算機技術(shù)的發(fā)展和進(jìn)步，使得幾乎所有工程領(lǐng)域復(fù)雜系統(tǒng)的行為都能用數(shù)值模型來模擬和預(yù)測[20]，因此，基于代理模型的全局靈敏度分析方法獲得了高度關(guān)注[20-22]。另一方面,基于方差的Sobol’指數(shù)蒙特卡洛仿真計算技術(shù)也得到了改進(jìn)[23-25]，雙循環(huán)重排序方法(Double Loop Reordering approach, DLR)作為Sobol’蒙特卡洛仿真技術(shù)改進(jìn)方法的一種，使得仿真計算效率得到很大的提高[26-29]。

然而，在數(shù)據(jù)存在污染的情形下，不能準(zhǔn)確反映樣本的統(tǒng)計特性，導(dǎo)致基于方差的靈敏度分析方法無法識別因子重要度或?qū)σ蜃又匾冗M(jìn)行排序，無法正確度量模型輸出不確定性的來源，從而使得模型更加復(fù)雜、計算成本更高。針對仿真試驗中數(shù)據(jù)污染(本文指數(shù)據(jù)中存在異常值)的問題，引入穩(wěn)健統(tǒng)計量改進(jìn)Sobol’法中的DLR蒙特卡洛仿真方法，提出了穩(wěn)健雙循環(huán)重排序方法(Robust Double Loop Reordering approach, RDLR)，并通過仿真試驗驗證了所提方法的有效性和穩(wěn)健性。

1 穩(wěn)健統(tǒng)計量及其性質(zhì)

1.1 穩(wěn)健統(tǒng)計量及其相對效率

1.2 兩組穩(wěn)健統(tǒng)計量的漸近分布特性

基于穩(wěn)健統(tǒng)計量HL-Shamos組合，可構(gòu)造如下統(tǒng)計量[30]：

(1)

基于穩(wěn)健統(tǒng)計量Med-MAD組合，可構(gòu)造如下統(tǒng)計量[33]：

(2)

為了驗證兩組穩(wěn)健統(tǒng)計量的分布情況，在正態(tài)分布假設(shè)條件下，分別在樣本量(n)為10和100時進(jìn)行1 000次隨機抽樣，基于樣本點繪制統(tǒng)計量TMM，THS的正態(tài)分布QQ圖，以及樣本量為10的經(jīng)驗累計分布圖和概率密度函數(shù)圖如圖1所示。

對比圖1a和圖1b可知，在樣本量較大時，兩組統(tǒng)計量更趨向于正態(tài)分布；在樣本量較小時，基于統(tǒng)計量THS的分布更趨向于正態(tài)分布。如圖1c所示為基于統(tǒng)計量THS、TMM以及均值—方差的經(jīng)驗累計分布圖，如圖1d所示為基于統(tǒng)計量THS、TMM以及均值—方差的密度函數(shù)，分析圖1c和圖1d可得，統(tǒng)計量THS、TMM的分布均近似服從正態(tài)分布，且基于統(tǒng)計量THS的分布更接近正態(tài)分布。因此，在樣本量較小時建議優(yōu)先選擇THS統(tǒng)計量。基于穩(wěn)健統(tǒng)計量的漸近正態(tài)分布特性，本文采用穩(wěn)健統(tǒng)計量代替均值和標(biāo)準(zhǔn)差估計樣本的位置參數(shù)和尺度參數(shù)是合理的。

2 基于穩(wěn)健統(tǒng)計量的全局靈敏度分析

2.1 穩(wěn)健統(tǒng)計量替代均值和方差的全局靈敏度分析

假設(shè)系統(tǒng)某一質(zhì)量特性Y和m個可控輸入變量x=(x1,…,xm)之間的關(guān)系為Y=f(x)，則模型輸出的全方差公式為：

V[Y]=Vxi[Ex-i[Y|xi]]+Exi[Vx-i[Y|xi]]。

其中：E[·]表示期望，V[·]表示方差，x-i表示除xi之外的輸入因子，Vxi[Ex-i[Y|xi]]用來定量地度量輸入因子xi對模型輸出的影響。因此，因子xi基于方差的一階全局靈敏度指數(shù)計算如下：

Si=Vxi[Ex-i[Y|xi]]/V[Y]。

當(dāng)模型輸出中含有異常值時，考慮采用穩(wěn)健統(tǒng)計量HL-Shamos或者M(jìn)ed-MAD替代均值—方差計算全局靈敏度指數(shù)，使得基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法在面對數(shù)據(jù)污染時仍能準(zhǔn)確識別模型輸入因子的重要度。

其中：a=1.158 75;a1=0.414 253 297;a2=0.442 396 799;當(dāng)n≤100時，a3=a4=0，當(dāng)n>100時,a3=2.822，a4=12.238。

其中：b=2.702 7;b1=-0.762 13;b2=-0.864 13;當(dāng)n≤100時，b3=b4=0,當(dāng)n>100，且n為奇數(shù)時，b3=0.299 6，b4=-149.357，當(dāng)n為偶數(shù)時,b3=-2.417，b4=-153.01。

基于穩(wěn)健的位置參數(shù)—尺度參數(shù)為HL-Shamos組合的全局靈敏度指數(shù)定義為Si1,i2,…,ik(HS)，

Si1,i2,…,ik(HS)=

(3)

基于穩(wěn)健的位置參數(shù)—尺度參數(shù)為Med-MAD組合的全局靈敏度指數(shù)定義為Si1,i2,…,ik(MM)，

Si1,i2,…,ik(MM)=

(4)

式(3)和式(4)中{i1,i2,…,ik}是{1,2,…,m}的一個子集，且1≤i1<…

2.2 基于穩(wěn)健統(tǒng)計量的RDLR蒙特卡洛仿真計算方法

假設(shè)函數(shù)f在積分空間Im=[0,1]m上平方可積，函數(shù)f通過Sobol’分解可分解為：

fi,j(xi,xj)+…+f1,2,…,m(x1,x2,…,xm)。

當(dāng)文獻(xiàn)[15]中的條件滿足時，上式的分解是唯一的，對其兩邊平方并積分，可得：

假定兩互補子集y和z構(gòu)成輸入變量x=(y,z)，令子集y=(xi1,xi2,…,xik)，子集y的方差

Sobol’法的總方差為

其中{i1,i2,…,ik}是{1,2,…,m}的一個子集，即1≤i1<…

在仿真試驗中，假定x和x′是樣本空間中兩個N×m維相互獨立的抽樣數(shù)組，令x=(y,z)，x′=(y′,z′)，對應(yīng)于輸入因子子集y=(xi1,…,xik)的方差為：

針對樣本均值和方差對異常值敏感的問題，提出穩(wěn)健雙循環(huán)重排序方法(RDLR)，采用穩(wěn)健統(tǒng)計量替代均值—方差改進(jìn)傳統(tǒng)的DLR方法。用穩(wěn)健的位置統(tǒng)計量估計每個分區(qū)中Nm個模型輸出的位置參數(shù)，即對第k(1≤k≤M)個分區(qū)分別計算HL統(tǒng)計量和中位數(shù)Med統(tǒng)計量來估計每個分區(qū)中模型條件輸出的位置參數(shù)，形式如下：

p,q=1,…,Nm,k=1,…,M；

由上式可得M個穩(wěn)健的條件輸出位置統(tǒng)計量，然后用穩(wěn)健的尺度參數(shù)統(tǒng)計量Shamos和MAD估計M個位置參數(shù)統(tǒng)計量的穩(wěn)健尺度參數(shù)：

(5)

(6)

模型非條件輸出的穩(wěn)健尺度參數(shù)為：

f(yq,zq)|))2/A(N),p,q=1,…,N；

(7)

(8)

因此，由式(5)～式(8)可得，輸入變量y=xi的基于穩(wěn)健統(tǒng)計量改進(jìn)的RDLR方法的一階全局靈敏度指數(shù)，即一階Sobol’指數(shù)的計算公式為：

Sy(HS)=Dy(HS)/DHS；

(9)

Sy(MM)=Dy(MM)/DMM。

(10)

2.3 RDLR方法的實現(xiàn)步驟

基于穩(wěn)健統(tǒng)計量的RDLR方法具有較高的穩(wěn)健性，能夠抵御異常值的影響，其仿真算法實現(xiàn)步驟如下：

步驟1算法初始化(測試函數(shù)選擇)，隨機生成N個m維的樣本點xj，計算對應(yīng)的輸出f(xj)。

步驟2令y=xi(1≤i≤m)獲得排序后的樣本點x(j),j=1,2,…,N，將排序后的輸入變量對應(yīng)的輸出f(xj)分成M(M

步驟4計算所有樣本點xj對應(yīng)的非條件輸出的穩(wěn)健尺度參數(shù)DHS和DMM。

步驟5計算基于穩(wěn)健統(tǒng)計量的輸入因子靈敏度指數(shù)Dy(HS)/DHS和Dy(MM)/DMM。

步驟6重復(fù)步驟2～步驟5，令y取每一個輸入變量，計算所有輸入變量的穩(wěn)健靈敏度指數(shù)。

傳統(tǒng)的DLR和改進(jìn)的RDLR方法能計算單個輸入變量的靈敏度指數(shù)，對變量子組無效。

3 算例分析

為了更好地說明所提方法的有效性，選用3個測試函數(shù)對提出的RDLR方法進(jìn)行驗證。Ishigami函數(shù)和Sobol G—函數(shù)都是全局靈敏度分析研究中廣泛采用的測試函數(shù)。Ishigami函數(shù)具有輸入變量間相互關(guān)系復(fù)雜，能夠代表大多數(shù)模型輸入變量間關(guān)系類型的特點[34]；Sobol G—函數(shù)輸入變量數(shù)目可變，并且可以通過調(diào)節(jié)函數(shù)表達(dá)式中非負(fù)常數(shù)，進(jìn)而達(dá)到改變模型輸入變量的靈敏度指數(shù)的目的，具有較高的靈活性[35]?；诜讲畹姆椒o法識別非線性偏態(tài)分布輸出函數(shù)的輸入變量的重要度[36]，為了驗證本文基于穩(wěn)健統(tǒng)計量的方法針對非正態(tài)問題的有效性，選擇該二維輸入的非線性偏態(tài)分布輸出函數(shù)。

考慮數(shù)據(jù)污染對基于方差的全局靈敏度分析方法的影響，人為增加一個初始輸出的污染函數(shù)，記為fcon(y,λ)，污染后的模型輸出為：fcon(y,λ)=y+m·s·I(y),數(shù)據(jù)污染的相關(guān)函數(shù)λ=(p,m,s)與3個參數(shù)有關(guān)，其中p表示樣本數(shù)據(jù)中變異數(shù)據(jù)的概率，本文取值分別為0(不含異常值)、0.05和0.1，此處給定污染數(shù)據(jù)占比均小于各統(tǒng)計量的失效點；m表示變異量，取值分別為100和300；s表示變異符號，取值為1和-1；其中I(y)為示性函數(shù)，即I(y≥0)=1，I(y<0)=-1?；诟鹘o定變異函數(shù)污染參數(shù)的取值，可得9種污染參數(shù)組合，分別命名為Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ、Ⅶ、Ⅷ和Ⅸ，第Ⅰ種變異概率為0，污染參數(shù)組合為(0,-,-)，模型輸出無污染,不同污染情形的參數(shù)組合如表1所示。

表1 不同數(shù)據(jù)污染情形的參數(shù)組合

3.1 Ishigami函數(shù)

Ishigami函數(shù)[34]包含3個輸入因子，表達(dá)式為：

f(x1,x2,x3)=sinx1+a·sin2x2+

其總方差D和條件方差D1,D2,D3的理論計算值如下所示：

令a=7，b=0.1，Ishigami函數(shù)基于方差的輸入因子x1,x2,x3的全局靈敏度指數(shù)解析解分別為：S1=0.313 8,S2=0.442 4,S3=0，因此S2>S1>S3,即因子x2重要度為第一，因子x1重要度為第二，因子x3重要度為第三。

如圖2所示為Ishigami函數(shù)1 000次仿真試驗輸出在9種不同的污染參數(shù)組合情形的均值、HL、中位數(shù)、標(biāo)準(zhǔn)差、Shamos和MAD等6個統(tǒng)計量的條形圖。

如圖2a所示為不同的位置參數(shù)統(tǒng)計量的條形圖，對比分析可知樣本均值對異常值非常敏感，而HL和中位數(shù)幾乎不受異常值影響；圖2b所示為不同尺度參數(shù)統(tǒng)計量的條形圖，其中標(biāo)準(zhǔn)差受到異常值的影響最大，穩(wěn)健尺度統(tǒng)計量Shamos和MAD受異常值影響較小。因此，面對數(shù)據(jù)污染問題，基于穩(wěn)健統(tǒng)計量的仿真試驗結(jié)果不易受異常值干擾。

如表2所示為Ishigami函數(shù)輸出在9組污染情形下的不同統(tǒng)計量的方差，受異常值影響，9組輸出的均值和標(biāo)準(zhǔn)差的方差遠(yuǎn)遠(yuǎn)大于穩(wěn)健統(tǒng)計量的方差，即均值和標(biāo)準(zhǔn)差對異常值更加敏感。

表2 Ishigami函數(shù)輸出在9種污染參數(shù)情形下不同統(tǒng)計量的方差

為驗證數(shù)據(jù)污染情形下基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法的有效性和穩(wěn)健性，本試驗分別在9種不同污染參數(shù)組合情形下各進(jìn)行1 000次獨立重復(fù)試驗，得到因子靈敏度指數(shù)。在每種污染參數(shù)組合條件下、每次獨立試驗中，在模型輸入空間中隨機抽取1 000個樣本點，計算得到1 000個輸出，執(zhí)行一次靈敏度指數(shù)計算，得到一組靈敏度指數(shù)。Ishigami函數(shù)輸入因子靈敏度指數(shù)基于方差、HL/Shamos和Med/MAD等方法在不同污染情形下，試驗所得靈敏度指數(shù)的均值統(tǒng)計結(jié)果如表3所示。

由表3可知，基于方差的全局靈敏度分析方法，僅在第I種污染參數(shù)組合，即輸出無污染的情形下能識別因子的重要度；在輸出受到污染的情形下，因子x1的靈敏度指數(shù)S1和因子x2的靈敏度指數(shù)S2幾乎相等，無法識別輸入因子重要度。基于HL/Shamos和Med/MAD的全局靈敏度指數(shù)在9種污染參數(shù)組合情形下，都能識別輸入因子的重要度，可得輸入因子靈敏度指數(shù)關(guān)系為S2>S1>S3，與解析解的結(jié)果一致；對比基于兩組不同穩(wěn)健統(tǒng)計量的因子靈敏度指數(shù)仿真試驗結(jié)果，基于HL/Shamos方法的性能更優(yōu)，各因子靈敏度指數(shù)更接近解析解。

表3 基于不同統(tǒng)計量的Ishigami函數(shù)在不同污染情形下因子靈敏度指數(shù)

如圖3所示為9種污染情形下，Ishigami函數(shù)基于不同統(tǒng)計量各進(jìn)行1 000次獨立重復(fù)試驗的因子靈敏度指數(shù)箱線圖。由圖3a可知，基于方差的方法在數(shù)據(jù)污染的情形下，不能有效識別因子重要度，基于方差方法對異常值敏感，仿真試驗結(jié)果受到異常值影響。圖3b和圖3c分別為基于穩(wěn)健統(tǒng)計量HL/Shamos和Med/MAD方法的靈敏度指數(shù)箱線圖，結(jié)果表明：在9種污染情形下基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法都能正確識別因子重要度，因子靈敏度指數(shù)關(guān)系為S2>S1>S3，與解析解的結(jié)果一致。因此，基于穩(wěn)健統(tǒng)計量的靈敏度指數(shù)計算方法能夠?qū)崿F(xiàn)數(shù)據(jù)污染情形下因子重要度排序，相比基于方差的方法更穩(wěn)健。比較圖3b與圖3c中第Ⅰ種污染參數(shù)組合情況可知在輸出無污染時，基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法同樣能夠正確識別輸入因子的重要度。

樣本量較小時，由于統(tǒng)計量THS相對于統(tǒng)計量TMM的分布更趨向正態(tài)分布，且統(tǒng)計量HL/Shamos相對統(tǒng)計量Med/MAD的效率更高，因此基于統(tǒng)計量HL/Shamos的全局靈敏度分析方法性能更優(yōu)。

3.2 非線性偏態(tài)分布輸出函數(shù)

輸出高度偏態(tài)分布的函數(shù)[36]：

y=x1/x2,x1～χ2(10),x2～χ2(13.978)。

輸入變量x1對函數(shù)輸出的貢獻(xiàn)比x2對函數(shù)輸出的貢獻(xiàn)大，然而基于方差的全局靈敏度分析方法無法識別出兩因子的重要度。

為了驗證輸出偏態(tài)分布且受到污染情形下基于方差、HL/Shamos和Med/MAD等方法的有效性和穩(wěn)健性，采用與3.1節(jié)相同的試驗程序，計算9種污染情形下基于不同統(tǒng)計量的輸入因子靈敏度指數(shù)。在顯著水平為0.05時，對9種污染條件下基于方差方法的兩因子1 000次仿真試驗的靈敏度指數(shù)進(jìn)行檢驗。結(jié)果顯示，除第I種污染參數(shù)組合外，其它污染參數(shù)組合下的P值都大于顯著水平，因此兩因子1 000次仿真試驗的靈敏度指數(shù)無顯著差異，不能對輸入因子進(jìn)行重要度排序，t檢驗P值如表4所示。

表4 不同污染情形下基于方差方法的因子靈敏度指數(shù)檢驗

如表5所示為輸出偏態(tài)分布的函數(shù)在9種污染情形下基于不同統(tǒng)計量的全局靈敏度分析方法，各進(jìn)行1 000次獨立重復(fù)試驗，所得因子的靈敏度指數(shù)均值統(tǒng)計結(jié)果。由表5可知，在函數(shù)輸出偏態(tài)分布或被污染情形時，基于方差全局靈敏度分析方法無法識別因子重要度，即因子x1的靈敏度指數(shù)S1和因子x2的靈敏度指數(shù)S2的關(guān)系為S1≈S2。對比基于穩(wěn)健統(tǒng)計量HL/Shamos和Med/MAD方法的靈敏度指數(shù)可知，因子x1的靈敏度指數(shù)S1顯著大于因子x2的靈敏度指數(shù)S2，因此基于穩(wěn)健統(tǒng)計量的靈敏度分析方法不受數(shù)據(jù)污染或偏態(tài)分布的影響，能正確識別因子重要度，具有較高的穩(wěn)健性。

表5 不同污染情形下基于不同統(tǒng)計量的因子靈敏度指數(shù)

如圖4所示為輸出偏態(tài)分布的函數(shù)在9種不同污染情形下，基于不同的統(tǒng)計量的全局靈敏度分析方法，各進(jìn)行1 000次獨立重復(fù)試驗，所得因子靈敏度指數(shù)的箱線圖。

由圖4a可知，在函數(shù)輸出呈偏態(tài)分布條件下，不論輸出是否存在污染的情形，基于方差的方法均無法識別因子的重要度，仿真試驗結(jié)果無法對因子重要度排序；圖4b和圖4c是基于穩(wěn)健統(tǒng)計量的靈敏度指數(shù)箱線圖，結(jié)果顯示函數(shù)輸出在9種污染情形下，輸入因子x1的全局靈敏度指數(shù)S1顯著大于輸入因子x2的全局靈敏度指數(shù)S2，因此基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法在數(shù)據(jù)偏態(tài)分布情形下，不論數(shù)據(jù)是否被污染都能準(zhǔn)確識別輸入因子的重要度；基于穩(wěn)健統(tǒng)計量的全局靈敏度分析具有較高的穩(wěn)健性和更廣泛的適用性。

3.3 Sobol G-函數(shù)

在較高維輸入和較高污染數(shù)據(jù)比例情形下，為了驗證本文基于穩(wěn)健統(tǒng)計量的全局靈敏度分析方法的有效性，選用全局靈敏度分析中常用的測試函數(shù)Sobol G—函數(shù)[35]，其形式為：

其中：輸入因子xi(i=1,2,…,d)在[-1,1]d上均勻分布，ai是非負(fù)常數(shù)。Sobol G—函數(shù)可以通過控制非負(fù)常數(shù)來調(diào)節(jié)每一個輸入對輸入方差的貢獻(xiàn)。當(dāng)非負(fù)常數(shù)ai小時，其對應(yīng)的輸入因子對輸出的方差貢獻(xiàn)就更大;當(dāng)非負(fù)常數(shù)ai大時，其對應(yīng)的輸入因子對輸出的方差貢獻(xiàn)就小。本文設(shè)置輸入因子數(shù)目為6，每個輸入因子對應(yīng)的非負(fù)常數(shù)是(a1,a2,a3,a4,a5,a6)=(15,3,2,1,0.1,0)，模型輸出的污染參數(shù)組合如表6所示。

表6 Sobol G—函數(shù)的輸出污染參數(shù)組合

考慮穩(wěn)健統(tǒng)計量HL和Shamos的極限BP點為29.3%，在此驗證數(shù)據(jù)中污染數(shù)據(jù)占比較高的情形下，基于穩(wěn)健統(tǒng)計量的全局靈敏度方法的有效性，設(shè)置污染參數(shù)最大比例為25%。如表7所示為不同污染情形下6維Sobol G—函數(shù)基于不同統(tǒng)計量的1 000次仿真試驗所得輸入變量靈敏度指數(shù)均值統(tǒng)計結(jié)果。

表7 不同污染情形下Sobol G—函數(shù)基于不同統(tǒng)計量的因子靈敏度指數(shù)

由表7可知，基于方差的方法僅在數(shù)據(jù)不存在污染的情形下能夠正確地對輸入因子重要度進(jìn)行排序，再次驗證了基于方差的全局靈敏度分析方法不適用于數(shù)據(jù)污染的情形?；诜€(wěn)健統(tǒng)計量的方法既能在數(shù)據(jù)無污染情形下識別因子重要度，又能在數(shù)據(jù)存在污染的情形下正確地識別因子的重要度。當(dāng)數(shù)據(jù)中污染值的比例達(dá)到25%時，基于穩(wěn)健統(tǒng)計量Med/MAD的方法優(yōu)于基于穩(wěn)健統(tǒng)計量HL/Shamos的方法，可能的原因是穩(wěn)健統(tǒng)計量Med/Mad的BP點較高，而25%的污染值的占比使得穩(wěn)健統(tǒng)計量HL/Shamos的性能較低。

4 結(jié)束語

針對數(shù)據(jù)污染情形下基于方差的全局靈敏度分析方法不能正確識別模型輸入因子重要度的問題。本文提出一種新的Sobol’指數(shù)蒙特卡洛仿真計算方法RDLR法。通過蒙特卡洛仿真計算兩個測試函數(shù)在不同的污染情形下基于不同方法的因子靈敏度指數(shù)，對比傳統(tǒng)的DLR方法和本文所提方法的穩(wěn)健性和有效性。分析試驗結(jié)果可得出以下結(jié)論：

(1)本文所提方法具有良好的穩(wěn)健性，在模型輸出存在異常值時，能正確識別因子重要度，有效抵御異常值的影響；

(2)RDLR法具有廣泛的適用性，在模型輸出理想無污染或高度非正態(tài)分布的情形下，能正確識別因子重要度；

(3)RDLR方法不必增加仿真試驗次數(shù)，保證仿真試驗的經(jīng)濟性；

(4)在模型輸入因子數(shù)目較小時，基于穩(wěn)健統(tǒng)計量HL/Shamos的方法性能更優(yōu)，當(dāng)數(shù)據(jù)中污染數(shù)據(jù)占比較高的情形下，基于穩(wěn)健統(tǒng)計量HL/Shamos的方法性能更優(yōu)。

本文主要貢獻(xiàn)是通過引入穩(wěn)健統(tǒng)計量改進(jìn)了Sobol’蒙特卡洛仿真計算方法，進(jìn)而提高了其抗異常值特性，拓展了全局靈敏度分析方法的適用范圍。需要指出的是，本文假設(shè)模型為單個輸出并且輸入因子相互獨立，當(dāng)模型有多個輸出或輸入因子具有相關(guān)關(guān)系時，穩(wěn)健的全局靈敏度分析方法將是一個值得關(guān)注的方向。