趙錄峰,呂震宙,張磊剛,王新維
(1.西北工業(yè)大學 航空學院, 陜西 西安 710072; 2.中國運載火箭技術研究院, 北京 100076;
3.中國人民解放軍93363部隊, 遼寧 沈陽 110141)
?
多輸出模型確認中的混合矩指標*
趙錄峰1,呂震宙1,張磊剛2,王新維3
(1.西北工業(yè)大學 航空學院, 陜西 西安710072; 2.中國運載火箭技術研究院, 北京100076;
3.中國人民解放軍93363部隊, 遼寧 沈陽110141)
摘要:在不確定性條件下,同時考慮到多維輸出之間的相關關系和單輸出的均值,構建由多輸出數(shù)學期望列陣和協(xié)方差矩陣組成的多輸出模型確認局部混合矩指標和全局混合矩指標。其中局部混合矩指標包括絕對指標(LA-3M)和相對指標(LR-3M),它們適合單點位置的多輸出局部模型確認;全局混合矩指標也包括絕對指標(GA-3M)和相對指標(GR-3M),它們適合多點位置的多輸出全局模型確認。通過數(shù)字算例和工程算例,所提指標可行有效,能夠方便地度量計算模型和物理實驗之間的差異程度。
關鍵詞:模型確認;混合矩;多輸出模型;相關性;不確定性
隨著計算機仿真技術的飛速發(fā)展和復雜產品工程實驗費用的不斷增加,很多復雜產品的物理實驗逐漸被計算模型所替代。計算模型是對真實物理過程的一個抽象描述,但建模過程受各種不確定性因素的影響,模型預測與物理實驗結果之間可能存在一定的差異。如何科學地定量度量它們之間的差異程度,是當前“模型確認”領域的重點研究內容,也是國內外學術界和工業(yè)界普遍關注和討論的熱點。模型確認被定義為:“從目標用途角度出發(fā),客觀地評估計算模型在多大程度上準確描述真實物理過程”[1-3]。在模型確認過程中,首要的核心問題是構建一個“模型確認指標”,該指標要求能夠科學有效地定量度量計算模型和實際物理過程之間的差異程度。模型確認指標被定義為“用來度量模型響應量與實驗結果之間差異程度的指標”[4-5]。它是模型確認工作的基礎和依據(jù),在工程設計與分析領域,廣泛用于多模型方案決策、模型預測能力評估和模型校準效果評價等工作之中。因此,開展模型確認指標研究,對推進模型確認工作的科學發(fā)展,提高計算機數(shù)字仿真技術的工程應用水平具有重要的理論價值和實踐意義。
近幾年國內外學術界對模型確認指標進行了大量研究,文獻[6]對現(xiàn)有指標進行了系統(tǒng)分類和分析。根據(jù)模型的關鍵特性,模型確認指標可分為不同的類型。按不確定性類型,模型確認指標可分為確定型和不確定型指標[7-8];按響應量的數(shù)目,可分為單輸出和多輸出指標;按輸入可控變量的位置,可分為單點位置和多點位置指標。本文主要研究不確定性條件下的模型確認指標構建問題。概括起來,現(xiàn)有的不確定性條件下的模型確認指標,主要分為四類:假設檢驗類、貝葉斯因子類、頻率類和面積類。通過對比分析,這些指標都有各自的優(yōu)缺點及適用范圍。文獻 [6,9]分別運用假設檢驗和貝葉斯方法開展模型確認工作,這兩種方法都是在一定置信水平下對模型預測與實驗結果的一致性進行判斷,但沒有對它們之間差異程度進行定量度量,實質上是一種定性與定量相結合的方法,并不是嚴格意義上的模型確認過程。文獻[1]提出的頻率指標是對一定置信水平下模型預測與實驗結果均值之間距離的度量,這一指標雖然是一種嚴格意義上的模型確認指標,且便于設計人員理解,但由于沒有考慮響應量的方差和相關性,因此僅適合響應量分布中心發(fā)生偏離情況下的模型的確認。文獻[1,8,10]運用模型輸出響應量與實驗結果經驗分布函數(shù)之間的面積差異,提出了直接面積和u-pooling兩種模型確認指標,適合于單維和多維獨立響應量的模型確認,但由于該方法建立在求解響應量邊緣分布面積差異的基礎之上,沒有考慮到多維響應量之間的相關性,因此不適合多維相關響應量的模型確認問題。為了解決這一問題,文獻[11]在直接面積和u-pooling指標的基礎上,應用多維概率積分轉換方法,提出了概率積分轉換(Probability Integral Transformation,PIT)和t-pooling模型確認指標,分別適合多維相關響應量在單點位置和多點位置的模型確認。雖然這一指標在多維相關響應量模型確認方面具有一定的優(yōu)越性,但也存在一定的缺陷。一是該指標采用多維概率積分轉換方法,壓縮了大量信息,存在評估過度和評估不足的風險;二是該指標屬于無量綱指標,物理意義不明確,不便于工程人員理解和使用;三是當實驗數(shù)據(jù)不足時,無法正確地求解多輸出響應量的聯(lián)合分布函數(shù),從而使得該指標較難應用。為此,本文針對多維相關響應量的模型確認問題,依據(jù)隨機變量各階矩的優(yōu)良特性,提出了一種物理意義相對明確且便于使用的基于混合矩的多維相關模型確認指標。
1混合矩模型確認指標構建
不確定性條件下,如果將多輸出計算模型的輸入變量或參數(shù)(如:載荷、幾何尺寸或邊界條件等)看作隨機變量,那么它們的輸出響應量(如:應力、應變、撓度或加速度等)也是多維隨機變量;同理,物理實驗數(shù)據(jù)也屬于多維隨機變量。因此,依據(jù)概率論和數(shù)理統(tǒng)計相關理論,不確定性條件下多輸出模型確認指標實質上就是對多維計算模型輸出響應量的聯(lián)合概率分布與實驗結果所服從的聯(lián)合概率分布之間的差異程度的度量[9]。
然而,在物理實驗和數(shù)學建模過程中,受隨機因素和人的認知能力影響,多輸出模型確認指標的構建比較復雜,面臨多種情況的挑戰(zhàn)[1]。一是多輸出響應量的物理意義可能不同,各自的量綱存在差異。如有的響應量用空間指標(應變、撓度)度量,有的用時空指標(速度、加速度)度量。二是多維響應量之間可能具有相關性。當多輸出模型具有相同的輸入變量,或輸入變量之間存在很強的相關性時,這些輸入變量之間的相關性就會傳遞到輸出響應量,使得輸出響應量之間存在很強的相關關系。三是多維輸出響應量的輸入變量可能是可控變量,也可能是隨機變量??煽剌斎胱兞吭谝欢〞r空范圍變化[12],當它在單點位置時,多維輸出響應量服從某一聯(lián)合概率分布;當它在多個位置變動時,多維輸出響應量所服從的聯(lián)合概率分布隨之發(fā)生變化。
隨機變量的數(shù)字特征是由隨機變量的分布確定的,能夠描述隨機變量某一方面的隨機取值特征的常數(shù)。對多維隨機變量來說,最重要的數(shù)字特征是數(shù)學期望(一階中心矩)、方差(二階中心矩)和協(xié)方差(二階混合矩)。數(shù)學期望描述隨機變量取值的平均大小,方差描述隨機變量與自身的數(shù)學期望的偏離程度,協(xié)方差描述隨機變量之間的相互關系。數(shù)學期望、方差和協(xié)方差雖然不像分布函數(shù)、分布律和概率密度函數(shù)那樣完整地描述隨機變量,但它能夠描述隨機變量的重要方面或人們最關心的重要特征。由于不確定性條件下的多輸出模型確認指標的本質是對模型輸出響應量與實驗觀察值所服從的兩個聯(lián)合分布函數(shù)之間差異程度的度量,然而在實際工程應用中,物理實驗觀察值的分布函數(shù)一般很難準確求得,因此,多維模型確認指標可以運用這兩個聯(lián)合分布函數(shù)的數(shù)字特征之間的差異程度來表示。
(1)
(2)
(3)
(4)
(5)
(6)
依據(jù)多維隨機變量各階矩的性質,結合文獻[1]和文獻[10]提出的模型確認指標構建原則,分別建立了如式(7)和式(8)所示的局部混合矩模型確認絕對指標(Local Absolute Metric based on Mixed Moment for model validation,LA-3M)和全局混合矩模型確認絕對指標(Global Absolute Metric based on Mixed Moment for model validation,GA-3M)。當多輸出模型的可控變量在單點位置時,可運用式(7)度量模型響應量與實驗結果之間差異程度;當可控變量在多個位置變動時,則每一個單點位置對應的dLA(Ym,Ye)屬于隨機變量,因此,可運用dLA(Ym,Ye)的數(shù)學期望dGA(Ym,Ye)度量模型響應量與實驗結果之間差異程度。
(7)
dGA(Ym,Ye)=E[dLA(Ym,Ye)]
(8)
(9)
dGR(Ym,Ye)=E(dLR(Ym,Ye))
(10)
上述確認指標不僅滿足文獻[1]和文獻[10]所提的客觀性、物理性、無界限、非負性、對稱性、三角不等性和收斂性等基本要求,而且同時考慮了輸出響應量的均值及二階混合矩,合理地度量了計算模型與實驗數(shù)據(jù)在均值、變異性及相關性等方面的差異程度,因此,該指標具有靈敏度高和便于工程應用等優(yōu)良特性。
2混合矩模型確認指標的求解
基于混合矩方法的多輸出模型確認指標求解的核心,是求解模型和實驗輸出響應量的一階中心矩和二階混合中心矩形成的協(xié)方差矩陣。下面依據(jù)大數(shù)定理,運用MonteCarlo數(shù)字模擬法求解。
基于混合矩方法的單點位置多輸出模型確認指標的求解流程如圖1所示,求解過程可概括為4個步驟。
圖1 單點位置多維響應量模型確認指標計算流程Fig.1 Flowchart of validation metric for models with multivariate output at a single validation site
步驟3:計算模型響應量和實驗觀測值的各階矩。采用式(11)~(14) 計算模型響應量和實驗觀測值的各階矩。
(11)
(12)
(13)
(14)
步驟4:求解模型確認指標。將求解所得的各階矩帶入式(7)和式(9),可得到單點位置多維響應量模型確認指標dLA(Ym,Ye)和dLR(Ym,Ye)。
多點位置指標求解建立在單點位置指標求解的基礎之上,主要包括3個步驟。
步驟1:根據(jù)實驗計劃,確定k個不同位置的可控輸入矢量xh,h=1,…,k。
步驟3:根據(jù)式(8)和式(10)求解多點位置的模型確認指標。其中dGA(Ym,Ye)和dGR(Ym,Ye) 分別由式(15)和式(16)求解。
(15)
(16)
3算例分析
下面針對基于混合矩的多輸出模型確認指標,通過兩個算例驗證它們的可行性和有效性。
本算例的物理實驗數(shù)據(jù)通過下式獲?。?/p>
(17)
式中:x(0≤x≤6)為可控變量,θ(θ=1.5)為模型參數(shù),ε1~N(0,0.22)和ε2~N(0,0.22)為兩個實驗響應量的測量誤差,ε1與ε2之間的相關系數(shù)ρε1,ε2=0.5。構建的預測模型共5個,它們分為兩組(見表1)。設模型響應量的樣本量用n表示,物理實驗觀測值的樣本量用m表示。
表1 兩個測試組預測模型
3.1.1測試1
包括3個預測模型。模型1與實驗過程一致,是一個正確的模型;模型2的參數(shù)θ與實驗模型有差異;模型3與實驗模型相比,不僅模型參數(shù)θ發(fā)生了變化,而且相關系數(shù)也發(fā)生了變化。顯然模型1的準確性高于模型2,模型2高于模型3。這組測試的目的:在模型參數(shù)或相關性存在差異的情況下,討論所建指標的適應性、可行性和有效性。
3.1.1.1單點位置測試
m固定,n變化。在x=2.0處,按照2.1部分提出的單點位置指標求解方法和步驟,分別由式(17)生成m=1000組實驗觀測數(shù)據(jù),由表1第一測試組模型生成n組模型響應量,其中n在區(qū)間[100,12 000]變化,計算模型確認指標值。不確定性條件下,指標值具有一定的隨機性,指標均值的置信度為95%的置信區(qū)間隨n的變化情況如圖2、圖3所示。
圖2 dLR均值的置信區(qū)間在x=2.0處隨n的變化曲線Fig.2 Confidence interval of the mean of dLRversus the n at validation site x=2.0
圖3 dLA均值的置信區(qū)間在x=2.0處隨n的變化曲線Fig.3 Confidence interval of the mean of dLAversus the n at validation site x=2.0
由如圖2、圖3可以看出,在單點位置,隨著n的增大,模型確認指標dLA(Ym,Ye)和dLR(Ym,Ye)的均值及其置信區(qū)間都迅速減小,當n>2000時,它們趨于穩(wěn)定狀態(tài),收斂到一個固定值或范圍,這表明n>2000時,模型響應量的計算結果趨于它們的解析解,且3個計算模型與物理實驗之間的差異程度也處于穩(wěn)定狀態(tài)。為了盡可能地提高模型響應量的計算精度,本文在模型確認指標的驗證過程中,選取n=10 000作為模型響應量樣本量。
n固定,m變化。在x=2.0處,按2.1部分提出的單點位置指標求解方法和步驟,由表1第一測試組模型生成n=10 000組模型響應量,由式(17)生成m組實驗觀測數(shù)據(jù),其中m在區(qū)間[5,1000]變化,計算模型確認指標值。指標均值的置信度為95%的置信區(qū)間隨m的變化情況如圖4和圖5所示。
(a) m在區(qū)間[0,50]變化圖(a) m varies on the interval [0,50]
(b) m在區(qū)間[50,1000]變化(b) m varies on the interval [50,1000]圖4 dLR均值的置信區(qū)間在x=2.0處隨m的變化曲線Fig.4 Confidence interval of the mean of dLRversus the m at validation site x=2.0
由圖4,圖5可以看出: 在單點位置,隨著實驗觀測數(shù)據(jù)m的增大,模型確認指標dLA(Ym,Ye)和dLR(Ym,Ye)的均值及其置信區(qū)間都迅速減??;當m>200時,它們開始趨于收斂;當m=1000時,它們收斂到一個固定值或范圍。
由圖4(a)和圖5(a)可以看出,當m在[5,50] 區(qū)間取值時,模型確認指標(尤其相對指標dLR)均值的置信區(qū)間較大,且相互重疊。這一現(xiàn)象表明當實驗觀測數(shù)據(jù)量過少時,運用該指標難以準確區(qū)分計算模型之間的優(yōu)劣,存在判斷錯誤的風險。因此,在實驗數(shù)據(jù)過少,且實驗數(shù)據(jù)誤差過大的情況下,提出的模型確認指標具有一定的局限性。
(a) m在區(qū)間[0,50]變化(a) m varies on the interval [0,50]
(b) m在區(qū)間[50,1000]變化(b) m varies on the interval [50,1000]圖5 dLA均值的置信區(qū)間在x=2.0處隨m的變化曲線Fig.5 Confidence interval of the mean of dLA versus the m at validation site x=2.0
由圖4(b)和圖5(b)可以看出,當實驗觀測數(shù)據(jù)m>50時,局部混合矩模型確認指標都能夠很好地度量計算模型與物理實驗之間的差異程度,客觀地評估計算模型之間的優(yōu)劣。為了從理論上更好地驗證所提指標的可行性和有效性,本文在模型確認指標的驗證過程中,選取實驗觀測數(shù)據(jù)的樣本量m=1000。
n固定,m固定。在x=2.0處,按2.1部分提出的單點位置指標求解方法和步驟,由式(17)生成1000組實驗觀測數(shù)據(jù),由表1第一測試組計算模型生成10 000組模型響應量,計算模型確認指標,并與文獻[11] 提出的PIT面積法指標進行對比,結果見表2。
表2 單點位置第一測試組模型確認指標計算結果
3.1.1.2多點位置測試
可控制變量x在(0≤x≤6)的區(qū)間范圍內,從0開始每間隔0.01取一個x值,即選取600個固定點。同理,對每個x值,由式(17)生成1000組實驗觀測數(shù)據(jù),由表1第一測試組中的3個模型分別仿真生成10 000組模型響應量,按照2.2部分提出的多點位置指標求解方法和步驟,得到的多點位置的模型確認指標值,并與文獻[11] 提出的t-pooling面積指標進行對比,結果見表3。
表3 多點位置第一測試組模型確認指標計算結果
本算例測試1中,通過表2和表3可以清晰地看出,無論在固單點位置,還是在多點位置,在模型參數(shù)和相關系數(shù)存在差異的情況下,本文提出的模型確認指標同PIT和t-pooling面積指標具有同樣的效果,都能夠客觀地判斷出模型1優(yōu)于模型2,模型2優(yōu)于模型3,這一結論與定性分析相符,達到第一組測試的目的,驗證了所提指標的可行性和有效性。
3.1.2測試2
包括2個預測模型。模型4和模型5將確定性參數(shù)θ誤作為不確定性變量,服從均值同為θ=1.5,方差不同的正態(tài)分布。其中模型4中θ的方差較小,而模型5中方差較大,顯然模型4與實驗數(shù)據(jù)的一致性高于模型5。第二組測試的目的:在模型參數(shù)的離散程度存在很小差異時,驗證所建指標是否能夠客觀準確地度量計算模型與物理實驗之間的差異程度。
3.1.2.1單點位置測試
在x=3.0處,按照2.1部分提出的單點位置指標求解方法和步驟,同理,由式(17)生成1000組實驗觀測數(shù)據(jù),由表1第二測試組計算模型生成10 000組模型響應量,得到模型確認指標,并與文獻[11] 提出的PIT面積指標進行對比,結果見表4。
3.1.2.2多點位置測試
在可控制變量x(0≤x≤6)的區(qū)間范圍,從0開始每間隔0.01取一個x值,即選取600個固定點。對應于每個x值,由式(17)生成1000組實驗觀測數(shù)據(jù),由表1第二測試組中的2個模型生成10 000組模型響應量,按照2.2部分提出的多點位置指標求解方法和步驟,得到模型確認指標值,并與文獻[11] 提出的t-pooling面積指標進行對比,結果見表5。
表4 單點位置第二測試組模型確認指標計算結果
表5 多點位置第二測試組模型確認指標計算結果
本算例測試2中,由表4和表5可以清晰地看出,無論在單點位置,還是在多點位置,在確定性參數(shù)θ被錯誤地作為不確定性變量的情況下,本文提出的指標同PIT和t-pooling面積指標具有同樣的效果,都能夠明確地判斷出模型4優(yōu)于模型5,這一結論與定性分析吻合,達到第二組測試的目的,驗證了本文所提指標的可行性和有效性。
對表2~5的計算結果進一步分析可以看出:基于混合矩的多輸出模型確認相對指標(LR-3M和GR-3M)的靈敏度最高,而多輸出模型確認絕對指標(LA-3M和GA-3M)與PIT和t-pooling指標的靈敏度處于同一水平,相對較低。按照模型確認指標構建原則,當模型1與實驗模型完全一致時,其指標值應該為零。從表2和表3的結果來看,模型1的六個指標值都不為零,尤其GR-3M的偏差較大(0.206),這一問題的原因不在于指標本身,而在于實驗室數(shù)據(jù)和模型計算存在不確定性。如果進一步提高模型計算精度,增加實驗數(shù)據(jù)量,減少實驗誤差,模型1的指標值最終會收斂于零。
如圖6所示的矩形截面懸臂梁,自由端承受集中可控力,梁的長度L=2m,截面寬度b=50mm,截面厚度h=37mm,彈性模量E~N(206.8,10.52),輸出響應量為梁固定端A點的正應力σA(GPa)、自由端B點處截面的轉角θB(rad)和撓度yB(mm),實驗數(shù)據(jù)由解析方程式(18)產生,它們的測量誤差分別為εθ~N(0,0.0052),εy~N(0,0.0052),εσ~N(0,0.012)。
圖6 懸臂梁結構Fig.6 Diagram of the cantilever beam structure
(18)
為節(jié)約實驗費用,提高效率,工業(yè)部門建立了3個不確定條件下的懸臂梁模型,這些模型的形式與解析方程相同,但各隨機變量的分布函數(shù)或分布參數(shù)存在差異,具體情況見表6。
表6 懸臂梁輸入隨機變量的分布參數(shù)
表7 懸臂梁模型確認指標計算結果
由表7可以看出,無論在單點位置,還是在多點位置,運用本文提出的模型確認指標都能夠明確判斷出模型1優(yōu)于模型2,模型2優(yōu)于模型3,這一結論與定性分析結果吻合,從而再次驗證了所提指標的可行性和有效性。
4結論
基于混合矩的多輸出模型確認絕對指標和相對指標及求解方法,能夠有效度量模型響應量與實驗結果之間差異程度。
1)由于采用了工程人員比較熟悉且便于應用的各階矩作為指標的組成要素,使得該指標具有工程應用方便、物理意義明顯和分辨率高等優(yōu)點。
2)指標既適合單輸出模型確認,也適合復雜的具有相關關系的多輸出模型確認,可以直接用于工程設計與分析、多模型方案決策、模型預測能力評估和模型校準效果評價等工作之中,對推進工程建模的建設和發(fā)展具有重要意義。
3)指標及其求解方法也存在一定的局限性:一是當實驗數(shù)據(jù)過少且誤差過大,或模型計算結果不準確時,指標計算結果會產生一定的偏差。二是指標求解方法是一種比較簡單的方法,在提高求解效率和準確性方面,指標的求解方法還有待于進一步深化研究。
參考文獻(References)
[1]Oberkampf W L, Barone M F. Measures of agreement between computation and experiment: validation metrics[J]. Journal of Computational Physics, 2006, 217(1): 5-36.
[2]Oberkampf W L, Sindir M,Conlisk A. Guide for the verification and validation of computational fluid dynamics simulations[J]. AIAA,1998.
[3]Oberkampf W L, Roy C J. Verification and validation in scientific computing[M].New York,USA:Cambridge University Press, 2010.
[4]Xiong Y, Chen W, Tsui K L,et al. A better understanding of model updating strategies in validating engineering models[J]. Computer Methods in Applied Mechanics and Engineering, 2009,198(15-16): 1327-1337.
[5]Messer M, Panchal J H, Krishnamurthy V,et al. Model selection under limited information using a value of information based indicator[J]. Journal of Mechanical Design, 2010,132(12): 1-13.
[6]Liu Y, Chen W, Arendt P,et al. Toward a better understanding of model validation metrics[J]. Journal of Mechanical Design, 2011,133(7): 1-13.
[7]Schwer L E. Validation metrics for response histories: perspectives and case studies[J]. Engineering with Computers, 2007,23(4): 295-309.
[8]Ferson S, Oberkampf W L. Validation of imprecise probability models[J]. International Journal of Reliability and Safety, 2009, 3: 3-22.
[9]Rebba R, Mahadevan S. Validation of models with multivariate output[J]. Reliability Engineering and System Safety, 2006, 91(8): 861-871.
[10]Ferson S, Oberkampf W L, Ginzburg L. Model validation and predictive capability for the thermal challenge problem[J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(29-32): 2408-2430.
[11]Li W, Chen W, Jiang Z, et al. New validation metrics for models with multiple correlated responses[J]. Reliability Engineering and System Safety, 2014, 127:1-11.
[12]Jiang X M, Mahadevan S. Bayesian validation assessment of multivariate computational models[J]. Journal of Applied Statistics, 2008, 35(1):49-65.
http://journal.nudt.edu.cn
Mixed moment validation metric for models with multivariate output
ZHAOLufeng1,LYUZhenzhou1,ZHANGLeigang2,WANGXinwei3
(1.School of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China;
2. China Academy of Launch Vehicle Technology, Beijing 100076, China; 3. The PLA Unit 93363, Shenyang 110141, China)
Abstract:Considering the relations among multi-outputs and the mean of single output, the mathematical expectation of single dimensional variable and covariance metric of multi-dimensional variables were introduced into the validation metrics for models. The new metrics of LA-3M and LR-3M were proposed for validating multi-responses at a single validation site, while the metrics of GA-3M and GR-3M were proposed to collect data of multiple responses observed at multiple validation sites. These metrics were examined through a numerical test case and an engineering example to illustrate their feasibility and effectiveness. Results show that the proposed metrics are efficient and they can easily measure the differential degree of multiple responses between calculation model and physical experiment.
Key words:model validation; mixed moment; models with multivariate output; correlation; uncertainty
中圖分類號:TB114.3
文獻標志碼:A
文章編號:1001-2486(2015)06-061-08
作者簡介:趙錄峰(1973—),男,陜西富平人,博士研究生,E-mail:zlf315611@126.com;呂震宙(通信作者),女,教授,博士,博士生導師,E-mail:zhenzhoulu@nwpu.edu.cn
基金項目:國家自然科學基金資助項目(51475370);高等學校博士學科點專項科研基金資助項目(20116102110003)
收稿日期:*2015-01-09
doi:10.11887/j.cn.201506013