基于概率關聯分析的仿真模型驗證方法研究

寧小磊, 趙新, 吳穎霞, 趙軍民, 呂梅柏, 陳韻

(1.中國華陰兵器試驗中心， 陜西華陰 714200;2.西北工業大學航天學院， 陜西西安 710075;3.西安現代控制技術研究所， 陜西西安 710065)

實彈飛行試驗是裝備鑒定/定型試驗檢驗、考核的最終手段，但隨著武器裝備單子樣價格的日益昂貴，現場實彈飛行試驗樣本越來越少，若單純依據現場試驗進行裝備鑒定/定型評估風險增大。除實彈飛行試驗外，仿真試驗是裝備鑒定/定型工作最重要的試驗手段[1]。實踐證明，在裝備鑒定/定型工作中綜合應用仿真技術，能夠有效提高裝備鑒定/定型質量、節省試驗經費、縮短試驗周期。

與實彈試驗不同，仿真試驗是基于模型的活動[2-4]，仿真模型可信度是仿真手段在裝備鑒定/定型應用中的生命線。目前要檢驗仿真模型可信度，最基本、最直接的方法是比較相同輸入條件下飛行試驗數據和仿真數據的一致性[4]。動態數據一致性檢驗方法主要有灰色關聯分析、TIC、相關系數法、誤差分析法、譜分析方法等[3]，這些方法原理簡單，計算易操作，但應用在裝備鑒定/定型工作時仍存在一些不足，主要是：這些方法僅適應單樣本(一維)試驗數據之間的一致性檢驗[5-8]，而裝備鑒定/定型模型驗證要檢驗單樣本(一維)飛行試驗數據與多樣本(多維)仿真數據的一致性。現有文獻[3-8]主要研究提高單樣本參考數據和單樣本比較數據的一致性檢驗精度和可靠性，并不適合解決裝備鑒定/定型面臨的單樣本實彈飛行數據和多樣本仿真數據這一類數據一致性檢驗問題。

針對上述問題，本文提出了一種基于概率關聯分析的仿真模型驗證方法。建立的概率關聯度模型和概率關聯分析方法，實現了由序列曲線關聯向立體面板關聯的方法拓展。建立的基于概率關聯分析的仿真模型驗證方法，解決了現有模型驗證方法應用時數據利用不充分、驗證算法不合適的問題。通過數值測試和應用實例驗證了本文方法的正確性和有效性。

1 問題描述與分析

在裝備鑒定/定型試驗中，分別獲取了實彈飛 行試驗數據x(i)(t),i=1,2,…,m和仿真數據y(i)(t),i=1,2,…,n,其中,上標i為第i次試驗;m為實彈飛行試驗次數;n為仿真試驗次數。根據試驗特點和數據來源，分析如下：

受試驗消耗、試驗周期等限制,尤其是價格昂貴的導彈類裝備鑒定/定型試驗中,現場實彈飛行科目設置時,1組試驗條件組合僅設置少量甚至1發導彈進行實彈試驗,實彈飛行試驗數據x是小樣本甚至是單樣本。此外實彈飛行試驗均是在確定性試驗條件下進行,但是確定性試驗條件(如射向、射角、射程、目標特性、使用環境、作戰方式等)僅是決定飛行試驗數據狀態的部分因素,大量隨機因素(如測試誤差、器件誤差、非線性誤差和漂移、安裝誤差、延時誤差、質心誤差等)同時也是決定飛行試驗數據狀態的重要因素,實彈飛行試驗數據x包含隨機因素的影響效應不全面。仿真技術具有經濟性、方便性和可重復性等優勢,使得在相同組合試驗條件下可以進行大量重復試驗,仿真試驗數據y是大樣本,同時仿真試驗數據y還包含有大量隨機因素對試驗數據狀態的影響信息。

根據上述分析,試驗數據的特點是:一般地,在一類試驗條件下,實彈飛行試驗數據x是單樣本數據,極少情況下是多樣本數據;仿真試驗數據y在這一類試驗條件下是多樣本數據。因此,裝備鑒定/定型模型驗證要解決的問題是:通過檢驗單樣本(維數=1)實彈試驗數據x和多樣本(維數>1)仿真數據y的一致性,判斷仿真系統是否可信。

2 現有模型驗證方法缺陷

動態數據一致性檢驗統一描述為2個步驟:

步驟1計算k時刻關聯系數,算子如下所示

R(k)=r(x(k),y(k))

(1)

式中:x,y為時間序列；r(·)為關聯系數計算函數;R(k)為k時刻關聯系數;k為時刻,k=1,2,…,T,T為時間序列長度。

步驟2計算綜合關聯度,算子如下所示

C(x,y)=c(R)

(2)

式中:c(·)為綜合關聯度計算函數;C(x,y)為綜合關聯度;R為關聯系數。

(1)式、(2)式便完成了動態數據的一致性檢驗,目前常見的動態數據一致性檢驗方法匯總見表1所示。

表1 常見的一致性檢驗方法

由表1可知,現有動態數據一致性檢驗方法能夠直接處理單樣本(一維)參考序列與單樣本比較序列之間的一致性問題,當檢驗數據為多樣本(多維)序列,或參考序列和比較序列其一是多樣本序列時,傳統算法會得到多個關聯度計算結果,如何綜合應用多個關聯度結果是一個亟待解決的新問題。而裝備鑒定/定型模型驗證是檢驗單樣本(一維)實彈飛行試驗數據x和多樣本(多維)仿真數據y的一致性,因此表1中的方法無法直接使用。如果選取多樣本仿真數據中的一條進行分析,勢必造成試驗信息浪費,模型驗證結果不充分。同時，若逐一關聯分析多樣本仿真數據后，如何應用出現的多個關聯度也是一個難題。

3 概率關聯度及基本性質

3.1 概率關聯度

為解決單樣本實彈飛行數據和多樣本仿真數據的一致性檢驗問題,提出了概率關聯度模型,在數據關聯分析時能夠同時處理多樣本序列,避免了現有方法抽取多樣本比較序列中的1條進行分析的不足，拓展了數據關聯分析的應用領域。

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

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

引理2設X、Y是不相關的連續隨機變量,其分布函數分別為F(X)、G(Y),則[F(X)、G(Y)]服從[0,1]上的二維均勻分布。

根據以上3個引理得到如下結論:

假設有2組數據

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

FY(k)為由數據Y(k)=[Y1(k),Y2(k),…,Yn(k)]T確定的概率分布函數,p(X(k))=FY(X(k))為X(k)的累積分布函數值,若2組數據X與Y一致,則有p(X(k)),k=1,2,…,T服從[0,1]上的均勻分布,其中,T為數據序列長度。

根據上述結論,單樣本實彈飛行數據Y和多樣本仿真數據Y的一致性檢驗可以轉換為如下數學問題的求解:判斷p(X(k)),k=1,2,…,T是否服從[0,1]上的均勻分布,若p(X(k))服從[0,1]上的均勻分布,則仿真模型通過驗證;否則未通過驗證。

根據以上分析,提出了概率關聯度,具體歸納為:

假設有參考序列x和比較序列y,按照以下步驟計算參考序列x和比較序列y的概率關聯度。

步驟1計算k時刻概率關聯系數poperator(k),算子如下所示

poperator(k)=FY(x(k))

(3)

步驟2步驟計算概率關聯度p,算子如下所示

p=c(poperator)

(4)

式中:FY(·)為由y(k)確定的概率分布函數;FY(x(k))為x(k)在分布函數FY(·)中的累積概率分布函數值;c(·)為概率關聯度模型中的綜合概率度計算公式,取均勻檢驗結果,工程實現時可以調用chi2gof函數。

3.2 概率關聯度的矩陣型數據應用改進

假設參考數據x為多樣本數據,即

(5)

為了計算矩陣型(多樣本)參考序列x和矩陣型比較序列y的概率關聯度,對(3)～(4)式描述的概率關聯度模型改進如下:

步驟1計算k時刻的概率關聯系數pi,operator(k),i=1,2,…,m,算子如下所示

pi,operator(k)=FY(xi(k))

(6)

步驟2計算概率關聯度p,算子如下所示

p=c(pi,operator)

(7)

(6)～(7)式中符號含義同(3)～(4)式。

3.3 概率關聯度的性質

定理1概率關聯度具有以下基本性質。

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

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

3) 可比性和唯一性;

4) 干擾因素獨立性。

證明

1) 根據(4)式、(7)式所示概率關聯度定義,其內涵描述的是T個概率關聯系數服從[0 1]均勻分布概率,其值范圍為[0 1],故有0≤p(x,y)≤1。

2) 若xi,xj∈x={xS|s=0,1,…,m;m≥2},對于任意xS1,xS2,一般由xS1,xS2確定的分布函數FxS1,FxS2不同,有FxS1(x(k))≠FxS2(x(k)),即p(xS1(k),xS2(k))≠p(xS2(k),xS1(k))成立,故模型滿足整體性。

3) 由于p(x,y)計算公式中不含有其他未知參數,因此具有唯一性和可比性。

4) 由于p(x,y)的計算過程,只與數據x和y有關,與其他因素無關,因此具有干擾因素獨立性。

定理2概率關聯度不滿足偶對稱性,即χ={x,y},有p(x,y)≠p(y,x)。

證明分別將x,y作為參考序列,根據整體性性質推導過程,一般有p(x(k),y(k))≠p(y(k),x(k))成立,故概率關聯度不滿足偶對稱性。

定理3概率關聯度模型不滿足數乘變換一致性和平移變換一致性。

證明對y分別作數乘變換和平移變換得y′=cy,y″=y+c。很顯然,由y確定的分布函數Fy與y′,y″確定的分布函數Fy′,Fy″一般并不一致,故x在y,y′,y″確定的分布函數中的累積分布概率值不同,故p(x,y)≠p(x,cy),p(x,y)≠p(x,y+c),即概率關聯度對數乘變換和平移變換均敏感,不滿足數乘變換和平移變換一致性。

4 基于概率關聯分析的仿真模型驗證方法

4.1 基本步驟

步驟1在相同初始條件下,分別得到實彈飛行試驗數據x和仿真試驗數據y。

(8)

(9)

式中:m為飛行試驗樣本量,當m=1時為單樣本實彈數據;n為仿真試驗樣本量;T為數據序列長度。

步驟2對飛行試驗數據x和仿真試驗數據y進行預處理,使其滿足等步長、等長度的數據序列要求。

步驟3計算k時刻由仿真試驗數據y(k)確定的分布函數FY(k)。

步驟4計算k時刻飛行試驗數據x(k)在FY(k)中的累積分布函數值p(k)=FY(x(k))。

步驟5檢驗p在一定置信水平α下是否服從[0,1]上的均勻分布。

步驟6若p服從[0,1]上的均勻分布,則模型通過驗證。否則模型未通過驗證,需進行仿真模型修正后,重新進行仿真模型驗證。

4.2 優勢分析

1) 在關聯系數計算時可以同時處理多樣本仿真數據,有效綜合了系統隨機效應影響,此外對實彈飛行數據也沒有維度限制(單樣本或多樣本均可),這正是本文方法區別于灰色關聯分析、TIC、譜分析、誤差分析法等現有方法的主要優勢,也是本文方法的核心進步點和創新點。

2) 因仿真步長和外測設備測試步長一般在0.02～0.04 s,因此試驗數據序列長度T一般會較大,參與本文方法中均勻分布檢驗的數據量很大,均勻性檢驗結果是穩定的,具有統計意義。

3) 灰色關聯分析關注曲線幾何形狀相似性、TIC、相關系數關注曲線距離、譜分析關注數據頻域組成與分布,而本文方法關注參考數據與比較數據是否來源于同一隨機過程,更符合實際情況。

5 數值測試

5.1 場景描述

為了檢驗本文方法的有效性,選取幾個典型案例進行了試驗測試。測試過程如下:

步驟1假定真實分布為F(x),從真實分布F(x)中直接抽取樣本x(i)(t)模擬實彈飛行數據,其中i=1,2,…,m,m為樣本量,t為飛行時間。

步驟2假定仿真試驗誤差為dF(y),用G(y)=F(x)+dF(y)模擬仿真試驗分布,從G(y)中抽取樣本y(i)(t)模擬仿真試驗數據,其中i=1,2,…,n,n為樣本量。很顯然,當dF(y)=0時,仿真試驗分布與真實分布相同。

步驟3給出如下假設:當dF(y)∈ε時,仿真模型滿足使用要求,其中ε為根據試驗要求確定的分布誤差閾值限,一般根據具體工程背景確定。

有效性測試判據為:當dF(y)∈ε條件滿足時,根據實彈數據x和仿真數據y給出仿真模型有效,即分布F(x)與G(y)一致;否則,當dF(x)?ε時,給出分布F(x)與G(y)不一致。

5.2 數值測試及結果分析

試驗1平穩隨機過程測試試驗

用F(t)=N(0,1)模擬實彈飛行抽樣過程,用S(t)=N(0+dμ,1+dσ)模擬仿真抽樣過程,取t=15 s,采樣周期設置為dt=0.02 s,同時取m=1,n=100。

圖1給出了dμ=0∶0.1∶1,dσ=0∶0.1∶1時各種條件組合下121次仿真模型驗證結果示意圖,數值測試中,調用chi2gof函數(取α=0.05)檢驗概率關聯系數是否服從[0 1]上的均勻分布。結果描述中,H表示仿真模型驗證結果(chi2gof函數的計算返回值),當H=0時,表示模型通過驗證;當H=1時,表示模型未通過驗證。從圖1可見,隨著dμ, dσ逐漸增大(仿真系統可信度降低)時,本文方法給出了模型未通過驗證的結論,說明本文方法的有效性。

圖1 H值(α=0.05) 圖2 100次Monte Carlo試驗結果(α=0.01)

為了測試本文方法的穩定性,分別選取dμ=0,dσ=0和dμ=1,dσ=1 2種場景,進行100次Monte Carlo試驗。圖2給出了100次Monte Carlo試驗后的概率關聯系數p(k)。從圖2可見,本文方法能夠正確給出模型驗證通過和未通過的結論。在dμ=0,dσ=0試驗中,概率關聯度p均非常小,H=0。在dμ=1,dσ=1的條件下,給出了96次H=1的結果,說明本文方法檢驗的穩定性。

試驗2非平穩隨機過程測試試驗

用以下非平穩隨機過程場景進行測試。

X(t)=N(μ,σ2)

(10)

式中:μ=20-0.001(0.2t-80)2;σ2=1。隨機過程X(t)均值μ(t)與t相關。

從X(t)直接抽樣模擬實彈數據,從以下隨機過程抽樣模擬仿真數據。

X(t)=N(μ+dμ,σ2)

(11)

圖3給出了ζ=0時的實彈數據和仿真數據。圖4給出了ζ=0時的概率關聯系數p(t)點圖。圖5給出了ζ=0時的概率關聯系數p(t)的頻率圖。圖6給出了ζ=1時的實彈數據和仿真數據。圖7給出了ζ=1時的概率關聯系數p點圖。圖8給出了ζ=1時的概率關聯系數p(t)的頻率圖。從圖4、圖5可見,實彈數據和仿真數據來源于同一非平穩隨機過程X(t)時,概率關聯系數p(t)均勻性很好。從圖7、圖8可見,當實彈數據和仿真數據來源分布不一致時,概率關聯系數p(t)均勻性較差,說明本文方法能夠處理這種情況,給出正確的判斷。

圖3 測試數據 圖4 ζ=0時概率關聯系數圖5 ζ=0時的概率關聯系數頻譜圖

圖6 測試數據 圖7 ζ=1時概率關聯系數圖8 ζ=1時的概率關聯系數頻譜圖

為了測試本文方法對非平穩隨機過程檢驗的穩定性,對2種場景分別進行100次Monte Carlo試驗。圖9給出了100次Monte Carlo試驗后的概率關聯度p。從圖9可見,本文方法能夠正確給出非平穩隨機過程條件下模型驗證通過和未通過的結論。

圖9 100次Monte Carlo試驗結果

試驗3非高斯隨機過程測試試驗

用以下非高斯隨機過程場景進行測試。

X(t)=N(μ,σ2)+β(1,i)

(12)

式中:μ=20-0.001(0.2t-80)2;σ2=1;i=t/0.02。

從X(t)直接抽樣模擬實彈數據,從以下隨機過程抽樣模擬仿真數據。

X(t)=N(μ+dμ,σ2)+β(1,i)

(13)

選取2種場景分別進行100次Monte Carlo試驗。圖10給出了部分實彈數據和仿真數據。圖11給出了100次Monte Carlo試驗后的概率關聯度p以及仿真模型驗證結果的H值圖譜。從圖11可見,本文方法能夠正確的給出非高斯隨機過程條件下模型驗證通過和未通過的結論,說明本文方法對非高斯隨機過程數據一致性檢驗適應。

圖10 測試數據

圖11 100次Monte Carlo試驗結果

試驗4 對比測試試驗

通過對比測試試驗分析本文方法與現有檢驗方法的性能,試驗場景同試驗1。進行200次試驗,模擬仿真抽樣樣本量為100,試驗時參數設置為:

圖12給出了各種檢驗方法的檢驗結果。從圖12可見,灰色關聯分析方法、相關系數方法、TIC方法、MSE方法給出了由100條曲線組成的曲面,這是因為根據這些方法的定義,分別計算了1組參考模擬數據和100組仿真數據的每一條數據的關聯分析結果,而本文方法給出的是一條關聯曲線結果,可見本文方法可以同時處理全部多組仿真數據。

圖12 不同方法的一致性檢驗結果

上述4個試驗從不同視角測試了本文所提概率關聯分析方法的有效性和正確性,分析上述試驗結果可得:本文方法可以處理平穩隨機過程數據、非平穩隨機過程數據、非高斯過程數據等問題情景,說明本文方法適應性廣泛,對數據特征基本不做要求,使用時不必對數據提前做平穩性、高斯分布等假設;本文方法可以直接處理單樣本數據與多樣本數據的一致性檢驗問題;隨著仿真模型可信度的降低,本文方法能夠正確給出定性的驗證結果。

6 應用實例

應用飛行試驗數據驗證某型反坦克導彈發動機的數字仿真模型。由于發動機推力在導彈飛行試驗中無法直接測試,我們通過導彈飛行試驗V-T數據和仿真試驗V-T數據的一致性檢驗驗證發動機模型。

試驗獲取了1次飛行試驗V-T數據和13次仿真試驗V-T數據,見圖13a)所示,其中紅色曲線為1次飛行試驗V-T數據,黑色曲線為13次仿真試驗V-T數據。圖13b)至13e)給出了灰色關聯分析(序號含義同文獻[3])、相關系數法、TIC、MSE方法的檢驗結果。圖13f)給出了概率關聯分析檢驗時的概率關聯系數圖,顯著水平取α=0.05,H的返回值為1,表示拒絕原假設,模型未通過驗證。分析原因,從圖13可見,仿真試驗V-T數據呈現一定的平行性,這并不符合實際情況。

從圖13可見,由于有13條仿真試驗V-T數據,灰色關聯分析、相關系數法、TIC、MSE方法等現有方法均給出了13個關聯計算結果,但如何綜合應用分析卻未給出策略,這也是這些方法在面向裝備鑒定/定型模型驗證應用時存在的不足。假設給定了判斷閾值,由于13個計算結果中有最大值和最小值,若出現判斷閾值介于最大值和最小值之間,將出現無法判斷的情況。本文方法由于關聯系數計算時,能夠同時處理13條仿真試驗數據,直接給出了驗證結論,此時H=1,模型未通過驗證。

圖13 不同方法的檢驗結果

7 結 論

本文面向裝備鑒定/定型應用的實際需求，針對裝備鑒定/定型試驗數據呈現的多維特點，提出了概率關聯度模型和基于概率關聯分析的仿真模型驗證方法。該方法能夠同時處理多樣本仿真數據，解決了裝備鑒定/定型工作中面臨的單樣本實彈飛行數據和多樣本仿真數據這一類一致性檢驗問題。在計算時融合了試驗過程的隨機因素，在小樣本實彈飛行試驗條件下更能充分利用試驗信息，克服了現有方法驗模時數據利用不充分、驗證算法不合適的不足。最后通過不同場景的數值測試和應用實例驗證了本文方法的合理性和有效性。