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

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

0 引言

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

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

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

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

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

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

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

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

(1)

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

(2)

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

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

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

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

假設(shè)系統(tǒng)某一質(zhì)量特性Y和m個(gè)可控輸入變量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對(duì)模型輸出的影響。因此，因子xi基于方差的一階全局靈敏度指數(shù)計(jì)算如下：

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

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

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

其中：b=2.702 7;b1=-0.762 13;b2=-0.864 13;當(dāng)n≤100時(shí)，b3=b4=0,當(dāng)n>100，且n為奇數(shù)時(shí)，b3=0.299 6，b4=-149.357，當(dāng)n為偶數(shù)時(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}的一個(gè)子集，且1≤i1<…

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

假設(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]中的條件滿足時(shí)，上式的分解是唯一的，對(duì)其兩邊平方并積分，可得：

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

Sobol’法的總方差為

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

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

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

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

由上式可得M個(gè)穩(wěn)健的條件輸出位置統(tǒng)計(jì)量，然后用穩(wěn)健的尺度參數(shù)統(tǒng)計(jì)量Shamos和MAD估計(jì)M個(gè)位置參數(shù)統(tǒng)計(jì)量的穩(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ì)量改進(jìn)的RDLR方法的一階全局靈敏度指數(shù)，即一階Sobol’指數(shù)的計(jì)算公式為：

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

(9)

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

(10)

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

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

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

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

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

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

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

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

3 算例分析

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

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

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

3.1 Ishigami函數(shù)

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

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

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

令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次仿真試驗(yàn)輸出在9種不同的污染參數(shù)組合情形的均值、HL、中位數(shù)、標(biāo)準(zhǔn)差、Shamos和MAD等6個(gè)統(tǒng)計(jì)量的條形圖。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 Sobol G-函數(shù)

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

其中：輸入因子xi(i=1,2,…,d)在[-1,1]d上均勻分布，ai是非負(fù)常數(shù)。Sobol G—函數(shù)可以通過控制非負(fù)常數(shù)來調(diào)節(jié)每一個(gè)輸入對(duì)輸入方差的貢獻(xiàn)。當(dāng)非負(fù)常數(shù)ai小時(shí)，其對(duì)應(yīng)的輸入因子對(duì)輸出的方差貢獻(xiàn)就更大;當(dāng)非負(fù)常數(shù)ai大時(shí)，其對(duì)應(yīng)的輸入因子對(duì)輸出的方差貢獻(xiàn)就小。本文設(shè)置輸入因子數(shù)目為6，每個(gè)輸入因子對(duì)應(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)計(jì)量HL和Shamos的極限BP點(diǎn)為29.3%，在此驗(yàn)證數(shù)據(jù)中污染數(shù)據(jù)占比較高的情形下，基于穩(wěn)健統(tǒng)計(jì)量的全局靈敏度方法的有效性，設(shè)置污染參數(shù)最大比例為25%。如表7所示為不同污染情形下6維Sobol G—函數(shù)基于不同統(tǒng)計(jì)量的1 000次仿真試驗(yàn)所得輸入變量靈敏度指數(shù)均值統(tǒng)計(jì)結(jié)果。

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

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

4 結(jié)束語(yǔ)

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

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

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

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

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

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