寧小磊, 吳穎霞, 趙新, 趙軍民, 張凱, 郝跳鋒,呂梅柏, 胡亞峰, 陳韻

(1.中國華陰兵器試驗中心, 陜西 華陰 714200; 2.西安現(xiàn)代控制技術(shù)研究所, 陜西 西安 710065；3.西北工業(yè)大學(xué) 航天學(xué)院, 陜西 西安 710072)

關(guān)聯(lián)分析是發(fā)現(xiàn)、查詢數(shù)據(jù)之間關(guān)聯(lián)性或相關(guān)性的一種實用分析技術(shù)，描述了事物中某些屬性同時出現(xiàn)的規(guī)律和模式[1-2]。關(guān)聯(lián)分析在經(jīng)濟學(xué)、軍事學(xué)、互聯(lián)網(wǎng)、航空航天及人工智能等領(lǐng)域有著廣泛應(yīng)用[3-4]，其中時間序列分析是數(shù)據(jù)關(guān)聯(lián)分析領(lǐng)域的熱點研究內(nèi)容。

對不同工程背景及應(yīng)用需求，目前已經(jīng)提出了多種數(shù)據(jù)關(guān)聯(lián)分析方法[5-6]。文獻[1]提出了一種信息不確定條件下的時間序列關(guān)聯(lián)分析法，解決了時間序列中含有不確定信息時的數(shù)據(jù)關(guān)聯(lián)分析問題。文獻[7]提出了一種灰色關(guān)聯(lián)分析方法，解決了時間序列曲線相似性關(guān)聯(lián)分析問題。常見的時間序列關(guān)聯(lián)分析方法還有TIC法、相關(guān)系數(shù)法、誤差分析法、EARTH方法等[7]，這些方法操作簡單，應(yīng)用廣泛，但主要適用于單樣本時間序列之間的分析，不適合工程中經(jīng)常遇到的多樣本(多元)時間序列情況[4]。針對多樣本關(guān)聯(lián)分析，文獻[4]提出了一種概率關(guān)聯(lián)度模型和概率關(guān)聯(lián)分析方法，實現(xiàn)了由序列曲線關(guān)聯(lián)向立體面板關(guān)聯(lián)的方法拓展，但計算時需要由比較序列樣本構(gòu)建累積分布函數(shù)。實際應(yīng)用中由于試驗抽樣的昂貴性、復(fù)雜性和困難性等，比較序列在實際情況中常常是小樣本數(shù)據(jù)容量，構(gòu)建的累積分布函數(shù)與真實分布差別一般較大，若直接應(yīng)用概率關(guān)聯(lián)度模型，關(guān)聯(lián)分析效果較差，嚴(yán)重影響了概率關(guān)聯(lián)度模型的工程應(yīng)用范圍。

針對上述問題，本文提出一種小樣本概率關(guān)聯(lián)度模型。使用樣本函數(shù)或樣本統(tǒng)計量代替樣本構(gòu)造累積經(jīng)驗分布函數(shù)并計算概率關(guān)聯(lián)系數(shù)，發(fā)展了概率關(guān)聯(lián)度模型的理論框架。當(dāng)比較序列樣本容量n<5時，使用樣本函數(shù)或樣本統(tǒng)計量構(gòu)建累積分布函數(shù)計算概率關(guān)聯(lián)系數(shù)；當(dāng)比較序列樣本容量n≥5時，使用Bootstrap方法重抽樣擴展樣本后構(gòu)建積累分布函數(shù)計算概率關(guān)聯(lián)系數(shù)，改進了概率關(guān)聯(lián)度模型中關(guān)聯(lián)系數(shù)的計算方法，從而提高了概率關(guān)聯(lián)分析的效率。通過仿真案例和實際應(yīng)用驗證了本文方法的正確性和有效性。

1 問題描述與分析

假設(shè)有參考序列x為

x=[x(1),x(2), …,x(k), …,x(T)]

(1)

比較序列y為

(2)

式中:k為時刻;T為序列長度。

參考序列x和比較序列y之間的數(shù)據(jù)關(guān)聯(lián)度為C,一般地,C可由如(3)式計算

C(x,y)=consistency(x,y)

(3)

式中:consistency(·)為一致性分析函數(shù),也稱一致性檢驗算子,一般有C(x,y)∈[0,1]。文獻[4]總結(jié)分析了目前常見的一致性分析函數(shù)consistency(·),并提出了一種適合多樣本時間序列的概率關(guān)聯(lián)算子和概率關(guān)聯(lián)分析方法。

參考序列x和比較序列y之間的概率關(guān)聯(lián)度計算步驟如下[4]:

步驟1計算k時刻概率關(guān)聯(lián)系數(shù)poperator(k)

poperator(k)=Fy(k)(x(k))

(4)

步驟2計算概率關(guān)聯(lián)度p

p=c(poperator)

(5)

式中:Fy(k)(·)為由比較序列樣本y1∶n(k)確定的累積分布函數(shù);Fy(k)(x(k))為x(k)在累積分布函數(shù)Fy(k)中的累積分布函數(shù)值;c(·)為概率關(guān)聯(lián)度模型中的綜合概率度計算公式,可取均勻檢驗結(jié)果。

從(4)～(5)式可見,步驟1中需要由比較序列樣本y1∶n(k)構(gòu)建累積分布函數(shù)Fy(k),當(dāng)比較序列y1∶n(k)樣本容量較大時(樣本量n較大時),可以構(gòu)建比較精確的累積分布函數(shù)Fy(k),綜合計算得到的概率關(guān)聯(lián)度p比較符合實際。當(dāng)比較序列y1∶n(k)樣本容量較小時,由比較序列y1∶n(k)構(gòu)建的累積分布函數(shù)Fy(k)因抽樣誤差與實際分布函數(shù)差距較大,此時,若應(yīng)用綜合計算概率關(guān)聯(lián)度p進行分析決策風(fēng)險較大。工程實際中,由于試驗抽樣的昂貴性、復(fù)雜性和困難性,經(jīng)常會遇到小樣本情況的比較序列數(shù)據(jù),限制了概率關(guān)聯(lián)度模型的實際應(yīng)用效果和應(yīng)用范圍。

2 小樣本概率關(guān)聯(lián)度模型

2.1 改進思路

綜上可知,小樣本條件下概率關(guān)聯(lián)度模型應(yīng)用遇到的主要問題是:比較序列y1∶n(k)樣本容量n有限,基于比較序列樣本y1∶n(k)不容易構(gòu)建精確的累積分布函數(shù)Fy(k),導(dǎo)致概率關(guān)聯(lián)度模型使用效果較差。解決該問題的有效方法是使用小樣本比較序列y1∶n(k)構(gòu)建滿足要求的累積分布函數(shù)Fy(k)。

Bootstrap方法是Efron于1979年提出的一種逼近復(fù)雜系統(tǒng)統(tǒng)計量估計值分布的統(tǒng)計方法[8],是目前被廣泛采用的小樣本數(shù)據(jù)處理方法,可以用來解決概率關(guān)聯(lián)度模型面臨的小樣本問題,即使用Bootstrap方法重抽樣擴展樣本后構(gòu)建積累分布函數(shù)計算概率關(guān)聯(lián)系數(shù)。但通常Bootstrap方法樣本容量以n≥5較合適[9-10],當(dāng)n<5時,其計算結(jié)果的任意性很大,構(gòu)造的經(jīng)驗分布函數(shù)和統(tǒng)計效果較差。針對樣本容量n<5的情況,文獻[3,9]使用了一種樣本順序比率統(tǒng)計量K,在樣本容量n<5的條件下應(yīng)用效果較好[3],認(rèn)為相對于樣本總體分布函數(shù)或經(jīng)驗分布函數(shù),使用樣本順序比率統(tǒng)計量K分布函數(shù)相容性檢驗效率更高,可見小子樣樣本構(gòu)建樣本順序比率統(tǒng)計量K經(jīng)驗分布函數(shù)比小子樣樣本構(gòu)建樣本經(jīng)驗分布函數(shù)精度要好。因此,可以使用樣本順序比率統(tǒng)計量K的經(jīng)驗分布函數(shù)代替概率關(guān)聯(lián)度模型中的經(jīng)驗分布函數(shù)。

2.2 變量函數(shù)概率關(guān)聯(lián)度模型改進

文獻[4]提出的概率關(guān)聯(lián)度模型,使用比較序列樣本y1∶n(k)直接構(gòu)建經(jīng)驗分布函數(shù),但樣本順序比率統(tǒng)計量K是比較序列樣本的函數(shù),不能直接使用。為了解決該問題,拓展概率關(guān)聯(lián)度模型的使用范圍,提出了變量函數(shù)概率關(guān)聯(lián)度改進模型,使其能同時應(yīng)用于樣本或樣本函數(shù)構(gòu)建經(jīng)驗分布函數(shù)。

為便于描述,首先引入2個引理。

引理1設(shè)X是一連續(xù)隨機變量,其分布函數(shù)為F(X),則F(X)服從[0,1]上的均勻分布。

引理2設(shè)Y=f(X)是一連續(xù)隨機變量,其分布函數(shù)為F(f(X)),則F(f(X))服從[0,1]上的均勻分布。

從上述2個引理可見,變量或變量函數(shù)(變量統(tǒng)計量)均可以代人概率關(guān)聯(lián)度模型中計算概率關(guān)聯(lián)系數(shù),因此對基本概率關(guān)聯(lián)度模型改進如下:

步驟1計算k時刻概率關(guān)聯(lián)系數(shù)poperator(k)

poperator(k)=Fy(k)(f(x(k)))

(6)

步驟2計算概率關(guān)聯(lián)度p,算子如下

p=c(poperator)

(7)

式中:f(·)為樣本函數(shù)或樣本統(tǒng)計量,一般為線性或非線性函數(shù)關(guān)系式,根據(jù)問題背景靈活選擇。

變量函數(shù)概率關(guān)聯(lián)度改進模型,可以采用樣本或樣本函數(shù)構(gòu)建經(jīng)驗分布函數(shù),來滿足小樣本概率關(guān)聯(lián)度模型改進使用需求。

注1:通過上述改進,概率關(guān)聯(lián)度模型適應(yīng)性更廣,但需要注意的是,并不是每一種樣本/變量函數(shù)f(·)都可以參與計算概率關(guān)聯(lián)系數(shù),這是因為函數(shù)計算會引入誤差,或者說抽樣誤差會通過函數(shù)計算放大,所以選擇合適的變量函數(shù)或樣本統(tǒng)計量很關(guān)鍵。(雖然很難,但仍可以尋找到一些有用統(tǒng)計量,即樣本函數(shù)f(·),在小樣本條件下應(yīng)用效果較好,比如文獻[9]找到的樣本順序比率統(tǒng)計量K。)

2.3 經(jīng)驗分布函數(shù)構(gòu)造

2.3.1 經(jīng)驗分布函數(shù)

設(shè)x(1),x(2),x(3), …,x(n)為來自分布函數(shù)F的隨機樣本,其經(jīng)驗分布函數(shù)Fn(x)定義為[11]

(8)

式中:I[·]為示性函數(shù);#A為集合A中元素的個數(shù)。Fn(x)為x的右函數(shù),共有n個跳躍點,跳躍度為1/n,即Fn(xi)-Fn(xi-1)=1/n,i=1,…,n,且有Fn(-∞)=0,Fn(+∞)=1。

2.3.2 樣本容量n<5經(jīng)驗分布函數(shù)構(gòu)造

設(shè)樣本x=[x(1),x(2),x(3), …,x(n+1)](其中1個樣本模擬參考序列樣本,其他樣本模擬比較序列樣本),對樣本按自小至大順序排列,得到樣本順序統(tǒng)計量x′=[x(1),x(2),x(3),…,x(n+1)],則樣本順序比率統(tǒng)計量Kijk=(x(j)-x(i))/(x(k)-x(i)),1≤i

當(dāng)n=2時,有1個統(tǒng)計量K123=(x(2)-x(1))/(x(3)-x(1)),且0≤K123≤1,對于樣本x來自總體X～N(μ,σ2),K123的累積分布函數(shù)為[3,9]

(9)

2.3.3 樣本容量n≥5經(jīng)驗分布函數(shù)構(gòu)造

2.3.3.1 經(jīng)典Bootstrap方法步驟

(10)

(11)

(12)

(13)

2.3.3.2 Bootstrap方法改進

Bootstrap方法通過大量再生子樣進行統(tǒng)計推斷,緩解了小樣本問題,但經(jīng)典的Bootstrap方法的采樣方式具有一定局限性,主要是:①樣本的累積經(jīng)驗分布函數(shù)將樣本的取值范圍限制在[x(1)x(n)]中,且樣本的取值是離散的,對于連續(xù)取值的變量無法獲取樣本點之外的信息。②從公式(10)可見,當(dāng)i=n或x=x(n)時,有pn=1,但理論上應(yīng)是當(dāng)n→∞,才有pn=1;同理x=x(1),有pn=0,但理論上應(yīng)是當(dāng)x→-∞,有pn=0。圖1給出不同樣本構(gòu)造的經(jīng)驗分布函數(shù)及真實分布函數(shù)對比。為了普適性應(yīng)用概率關(guān)聯(lián)度模型,對基本Bootstrap方法進行修正,改進的主要思路是:

1) 使用樣條函數(shù)代替原經(jīng)驗分布函數(shù)構(gòu)造使用的階躍連接,從而解決了對于連續(xù)取值的變量無法獲取樣本點之外的信息問題。

圖1 經(jīng)驗分布函數(shù)

2.4 小樣本概率關(guān)聯(lián)度模型

按照以下步驟計算參考序列x和比較序列y之間的小樣本概率關(guān)聯(lián)度。

步驟1計算k時刻概率關(guān)聯(lián)系數(shù)poperator(k)

(14)

步驟2計算概率關(guān)聯(lián)度p

p=c(poperator)

(15)

2.5 小樣本概率關(guān)聯(lián)度模型的基本性質(zhì)

性質(zhì)1小樣本概率關(guān)聯(lián)度具有以下基本性質(zhì)。

1) 規(guī)范性,即0≤p(x,y)≤1;

2) 整體性,對于不同的相關(guān)因素序列xi,xj,一般有p(xi,xj)≠p(xj,xi),i≠j;

3) 可比性和唯一性;

4) 干擾因素獨立性。

性質(zhì)2概率關(guān)聯(lián)度不滿足偶對稱性,即χ={x,y},有p(x,y)≠p(y,x)。

性質(zhì)3概率關(guān)聯(lián)度模型不滿足數(shù)乘變換一致性和平移變換一致性。

3 適用于小樣本問題的概率關(guān)聯(lián)分析步驟

步驟1在相同初始條件下,分別得到參考序列x和比較序列y。

(16)

(17)

式中:m為參考序列x1∶m(k)樣本容量,n為比較序列y1∶n(k)樣本容量，此處可以使用矩陣型概率關(guān)聯(lián)度模型[4]。

步驟2對參考序列x和比較序列y進行預(yù)處理,使其滿足等步長、等長度的數(shù)據(jù)序列要求。

步驟3計算k時刻概率關(guān)聯(lián)系數(shù)poperator(k)

1) 當(dāng)y1∶n(k)樣本容量n<5時

①構(gòu)建樣本順序比率統(tǒng)計量Kijk=(x(j)-x(i))/(x(k)-x(i)),1≤i

②將參考樣本x(k)帶入經(jīng)驗分布函數(shù)Fy1∶n(K)計算累積分布函數(shù)值,得到關(guān)聯(lián)關(guān)聯(lián)系數(shù)poperator(k)。

2) 當(dāng)y1∶n(k)樣本容量n≥5時

②將參考樣本x(k)帶入經(jīng)驗分布函數(shù)Fy1∶n(K)計算累積分布函數(shù)值,得到關(guān)聯(lián)系數(shù)poperator(k)。

步驟4決策。檢驗poperator在一定置信水平α下是否服從[0,1]上的均勻分布。若通過檢驗,說明通過關(guān)聯(lián)分析。否則,未通過概率關(guān)聯(lián)分析。

4 仿真測試與分析

通過幾個仿真測試案例驗證本文改進模型的正確性和有效性,仿真案例分別測試樣本容量n=2(Kijk有具體理論分布解析式)、n=3(Kijk沒有具體理論分布解析式)和樣本容量n≥5時(Bootstrap方法改進)的應(yīng)用場景。

仿真1參考時間序列X和比較時間序列Y為:

X=[0.053 0.059-1.912 0.402 1.837-0.618

-0.043 1.125];

Y=[-0.437-2.843 0.846 0.196 0.870-0.524

1.573-0.191;-0.014-0.465-0.950-1.404

-1.257 0.447-1.328-1.241]

檢驗參考時間序列X和比較時間序列Y的一致性。

比較時間序列Y樣本容量n=2,有1個統(tǒng)計量K123=(x(2)-x(1))/(x(3)-x(1)), 且0≤K123≤1,累積分布函數(shù)為

(18)

圖2給出了參考時間序列和比較時間序列樣本;將參考時間序列X帶入(18)式所示累積分布函數(shù),得到各個時刻的概率關(guān)聯(lián)系數(shù),如圖3所示;圖3b)給出了排序后的概率關(guān)聯(lián)系數(shù)。最后計算得到的概率關(guān)聯(lián)度P=0.608 3,結(jié)果:H=0。其中,H表示仿真模型驗證結(jié)果(kstest函數(shù)的計算返回值),當(dāng)H=0時,表示兩組數(shù)據(jù)一致;當(dāng)H=1時,表示兩組數(shù)據(jù)不一致。

圖2 樣本序列圖3 概率關(guān)聯(lián)系數(shù)

仿真2參考時間序列X和比較時間序列Y如下:

X=[0.503 1.231 1.103 0.078 0.447-0.280

-0.419-0.384-1.851-0.406];

Y=[0.607 0.337 0.119 0.311 0.623-0.824

1.801-2.039 0.029 0.745;0.333-0.352

0.356-0.914 1.588 1.765 1.429-1.196

-0.479 1.567;0.985 0.968 1.506 1.341

1.526 0.403-1.318-0.110 0.567 0.213]

檢驗參考時間序列X和比較時間序列Y的一致性。

由于比較時間序列Y樣本容量n=3時,樣本順序比率統(tǒng)計量K沒有顯式表達式,首先通過數(shù)值模擬的方法給出了累積概率密度分布如圖4所示。得到累積分布函數(shù)后,計算參考樣本在其中的函數(shù)值便可得到概率關(guān)聯(lián)系數(shù)。圖5給出了各個時刻的概率關(guān)聯(lián)系數(shù)。最后計算得到的概率關(guān)聯(lián)度P=0.716 1,結(jié)果:H=0。

圖4 n=3時的累積概率密度分布

圖5 概率關(guān)聯(lián)系數(shù)

仿真3參考時間序列X和比較時間序列Y如下:

X=[-1.769 2.155 2.174-0.493-0.268-0.070

-0.803-0.992-1.147-1.917];

Y=[-1.241 1.583 2.739 0.508 2.029 1.004

-0.348-0.546 2.404 0.709;-1.874 2.420

2.878 0.228 0.493 0.434-1.546-0.184

-0.560 0.118;-1.233 2.319 2.485 0.025

-0.559-0.757-0.210-0.699-1.800-1.646;

-1.027 2.307 2.529-1.065 1.811 0.859

1.098 0.295-0.840-0.931;-1.042 1.667

2.707 0.148-0.364 0.194-0.355-1.012

1.490-0.653;-1.619 2.085 2.432

1.055-0.943-0.559-1.386 0.343-0.475

0.568;-1.616 2.210 2.019 1.827 1.334

1.828-0.350 0.440 0.479-0.776]

檢驗參考時間序列X和比較時間序列Y的一致性。

由于k時刻比較序列Y樣本容量n≥5,首先使用Bootstrap方法重抽樣,然后基于重抽樣樣本構(gòu)建累積分布函數(shù),部分結(jié)果如圖6所示。由圖6可見,改進的自助樣本構(gòu)建的累積概率分布函數(shù)明顯比由原始樣本構(gòu)建的累積概率分布函數(shù)光滑。各個時刻的概率關(guān)聯(lián)系數(shù)計算結(jié)果為[0.162 0.455 0.164 0.152 0.911 0.295 0.286 0.115 0.116 0.015],最后計算得到的概率關(guān)聯(lián)度P=0.073,結(jié)果:H=1。

圖6 部分時刻累積分布函數(shù)

5 彈道一致性檢驗應(yīng)用

炮射導(dǎo)彈是由坦克炮發(fā)射的一種精確制導(dǎo)武器，提高了坦克炮的遠距離精確打擊能力。研究不同狀態(tài)下炮射導(dǎo)彈的彈道一致性有利于部隊訓(xùn)練使用。由于炮射導(dǎo)彈價格的昂貴性，現(xiàn)場試驗組織的復(fù)雜性，現(xiàn)場飛行試驗樣本容量一般為小子樣。假設(shè)高原、平原2種狀態(tài)下試驗數(shù)據(jù)如表3所示，根據(jù)炮射導(dǎo)彈的彈道特征[13]，選擇了3個典型彈道特征點(第一波谷、第一波峰、平穩(wěn)點)進行一致性檢驗。

表1 試驗數(shù)據(jù)

表2 檢驗結(jié)果

檢驗結(jié)果見表2，有理由認(rèn)為這2種狀態(tài)滿足彈道一致性，這與實際現(xiàn)場試驗結(jié)果相符。

6 結(jié) 論

針對小樣本條件下的數(shù)據(jù)關(guān)聯(lián)分析需求，研究提出了一種小樣本概率關(guān)聯(lián)度模型。使用樣本函數(shù)或樣本統(tǒng)計量的累積經(jīng)驗分布函數(shù)代替由樣本構(gòu)造的累積經(jīng)驗分布函數(shù)，發(fā)展拓展了概率關(guān)聯(lián)度模型。當(dāng)比較序列樣本容量n<5，使用樣本樣本順序統(tǒng)計量的經(jīng)驗分布函數(shù)計算概率關(guān)聯(lián)系數(shù)；當(dāng)比較序列樣本容量n≥5，使用Bootstrap方法重抽樣構(gòu)造經(jīng)驗分布函數(shù)計算概率關(guān)聯(lián)系數(shù)，解決了小樣本比較序列構(gòu)造的累積概率分布函數(shù)誤差較大導(dǎo)致概率關(guān)聯(lián)系數(shù)計算不準(zhǔn)確的問題。仿真案例和實際應(yīng)用驗證了本文方法的合理性和有效性。