非平穩(wěn)時(shí)序數(shù)據(jù)的分段辨識(shí)及其遞推算法①

黃雄波

(佛山職業(yè)技術(shù)學(xué)院 電子信息系,佛山 528137)

非平穩(wěn)時(shí)序數(shù)據(jù)的分段辨識(shí)及其遞推算法①

黃雄波

(佛山職業(yè)技術(shù)學(xué)院 電子信息系,佛山 528137)

在實(shí)際生活中,廣泛地存在著一類(lèi)在整體上屬于非平穩(wěn)但又可轉(zhuǎn)化為數(shù)段局部平穩(wěn)的時(shí)序數(shù)據(jù),對(duì)該類(lèi)非平穩(wěn)時(shí)序數(shù)據(jù)的辨識(shí)問(wèn)題進(jìn)行了研究,并提出了一種具有遞推機(jī)制的分段辨識(shí)算法.該算法從平穩(wěn)時(shí)序數(shù)據(jù)的定義出發(fā),以均值、方差及自相關(guān)函數(shù)等數(shù)字統(tǒng)計(jì)特征為校驗(yàn)統(tǒng)計(jì)量,構(gòu)造了具有遞推機(jī)制的均值突變點(diǎn)、方差突變點(diǎn)及自相關(guān)函數(shù)突變點(diǎn)的析出算法,在此基礎(chǔ)上,從被辨識(shí)的非平穩(wěn)序列中劃分出數(shù)段局部平穩(wěn)的子序列,進(jìn)一步,應(yīng)用Burg算法對(duì)各局部平穩(wěn)子序列進(jìn)行了自回歸的遞推辨識(shí).實(shí)驗(yàn)表明,新設(shè)計(jì)的算法能以較小的位置偏差析出各局部平穩(wěn)子序列的分界點(diǎn),同時(shí),在保證較高精度的辨識(shí)條件下,計(jì)算效能獲得了顯著的提升.

非平穩(wěn)時(shí)序數(shù)據(jù);局部平穩(wěn);分段辨識(shí);遞推;分界點(diǎn);自回歸模型

1 引言

時(shí)序數(shù)據(jù)根據(jù)其統(tǒng)計(jì)結(jié)構(gòu)是否具有時(shí)變的特點(diǎn),可分為平穩(wěn)序列和非平穩(wěn)序列兩大類(lèi),通常,人們所獲取到的時(shí)序數(shù)據(jù)大多為非平穩(wěn)序列,也就是說(shuō),其有關(guān)的統(tǒng)計(jì)特性是隨著時(shí)間t的變化而變化的.對(duì)于平穩(wěn)時(shí)序數(shù)據(jù)而言,其線(xiàn)性辨識(shí)模型主要有三種, AR(自回歸)模型、MA(滑動(dòng)平均)模型和自回歸滑動(dòng)平均(ARMA)模型,這些模型可以相互轉(zhuǎn)換,而辨識(shí)的主要任務(wù)有兩個(gè)方面,辨識(shí)模型階次的確定和辨識(shí)模型參數(shù)的估計(jì).近年來(lái),眾多專(zhuān)家學(xué)者圍繞這三種辨識(shí)模型展開(kāi)了系統(tǒng)而深入的研究[1-9],例如,張賢達(dá)針對(duì)帶有高斯ARMA噪聲的平穩(wěn)時(shí)序數(shù)據(jù),提出了一種高階累量的MA模型辨識(shí)算法;熊淵博研究了一種線(xiàn)性算法,將ARMA模型辨識(shí)問(wèn)題中的非線(xiàn)性求解問(wèn)題化簡(jiǎn)為求解兩組線(xiàn)性方程組;周毅等通過(guò)用AR模型等價(jià)ARMA模型的思想,提出了一種確定ARMA模型參數(shù)估計(jì)方法,并給出了計(jì)算等價(jià) AR模型參數(shù)估計(jì)的依階次遞增遞推算法.

相對(duì)地,非平穩(wěn)時(shí)序數(shù)據(jù)還沒(méi)有形成統(tǒng)一而完整的分析方法[10,11].據(jù)此,在實(shí)際應(yīng)用中,通常是運(yùn)用適當(dāng)?shù)臄?shù)學(xué)變換把非平穩(wěn)時(shí)序數(shù)據(jù)轉(zhuǎn)化為平穩(wěn)時(shí)序數(shù)據(jù),并以此作為實(shí)際物理過(guò)程的近似.例如,經(jīng)典的博克斯-詹金斯辨識(shí)方法就是通過(guò)對(duì)非平穩(wěn)時(shí)序數(shù)據(jù)進(jìn)行有限次的差分處理,把非平穩(wěn)序列化為平穩(wěn)序列;朱學(xué)鋒等基于Mallat快速小波算法對(duì)非平穩(wěn)時(shí)序數(shù)據(jù)進(jìn)行了分解和重構(gòu),進(jìn)而獲得了趨勢(shì)項(xiàng)序列和剔除趨勢(shì)項(xiàng)后的零均值平穩(wěn)時(shí)序數(shù)據(jù)[12];林樹(shù)寬等針對(duì)時(shí)序數(shù)據(jù)的非平穩(wěn)性,通過(guò)經(jīng)驗(yàn)?zāi)Ｊ椒纸獾玫搅巳舾蓚€(gè)平穩(wěn)序列和趨勢(shì)項(xiàng),在此基礎(chǔ)上,對(duì)每個(gè)平穩(wěn)序列建立相應(yīng)的辨識(shí)模型[13].

在實(shí)際生活中,廣泛地存在著一類(lèi)在整體上屬于非平穩(wěn)但又可轉(zhuǎn)化為數(shù)段局部平穩(wěn)的時(shí)序數(shù)據(jù), P.M.Djuric等最早應(yīng)用貝葉斯法對(duì)該類(lèi)辨識(shí)問(wèn)題進(jìn)行了研究,并導(dǎo)出了一個(gè)關(guān)于分段數(shù)、各段自回歸模型階數(shù)和各段之間分界點(diǎn)的優(yōu)化方程[14];在此基礎(chǔ)上,王文華等推導(dǎo)出一些具有遞歸關(guān)系的求解表達(dá)式,從而有效地簡(jiǎn)化了優(yōu)化方程的求解[15];為進(jìn)一步提升計(jì)算效能,陳穎等結(jié)合自回歸AR(Auto regressive)模型本身固有的特性,提出了直接遞推多維聯(lián)合分布概率的優(yōu)化方程求解方法[16].總體而言,基于貝葉斯法框架的辨識(shí)算法仍然存在著計(jì)算復(fù)雜的效率問(wèn)題,究其原因是由于在估計(jì)最優(yōu)劃分參數(shù)值時(shí)需要計(jì)算多維條件的分布概率.據(jù)此,本文擬從平穩(wěn)時(shí)序數(shù)據(jù)的定義出發(fā),設(shè)計(jì)實(shí)現(xiàn)一種新的具有遞推機(jī)制的分段辨識(shí)算法,該算法的主要思想是:以均值、方差及自相關(guān)函數(shù)等數(shù)字統(tǒng)計(jì)特征為校驗(yàn)統(tǒng)計(jì)量,構(gòu)造具有遞推機(jī)制的均值突變點(diǎn)、方差突變點(diǎn)及自相關(guān)函數(shù)突變點(diǎn)的析出算法,進(jìn)而能快速有效地析出了各局部平穩(wěn)子序列的分界點(diǎn),最后基于自回歸AR模型對(duì)這些局部平穩(wěn)子序列進(jìn)行遞推辨識(shí).實(shí)驗(yàn)結(jié)果表明,該辨識(shí)算法在保證辨識(shí)精度的基礎(chǔ)上獲得了顯著的計(jì)算效能的提升.

2 問(wèn)題描述

2.1 非平穩(wěn)時(shí)序數(shù)據(jù)的數(shù)字統(tǒng)計(jì)特征描述

從定義1易知,若某一時(shí)序數(shù)據(jù)為非平穩(wěn)時(shí)序數(shù)據(jù),則它的均值、方差及自相關(guān)函數(shù)等數(shù)字統(tǒng)計(jì)特征將部分或全部地隨時(shí)間t的變化而變化.據(jù)此,可用如下的時(shí)變函數(shù)對(duì)上述非平穩(wěn)時(shí)序數(shù)據(jù)的統(tǒng)計(jì)特征進(jìn)行描述:

均值為:

方差為:

自相關(guān)函數(shù)為:

2.2 分段局部平穩(wěn)時(shí)序數(shù)據(jù)的辨識(shí)問(wèn)題

1)在0￡t￡n范圍內(nèi)具有非平穩(wěn)的特征;

從定義2可知,分段局部平穩(wěn)時(shí)序數(shù)據(jù)Yt的部分或全部統(tǒng)計(jì)特性在一些時(shí)刻(如定義中的e,f,…,g)里發(fā)生突變,但對(duì)各獨(dú)立的子序列而言,它們又為平穩(wěn)子序列.通常,也把這些突變時(shí)刻稱(chēng)為分界點(diǎn).

在現(xiàn)實(shí)生活中,嚴(yán)格意義上的平穩(wěn)序列是很難找到的,但很多真實(shí)序列在某一局部的時(shí)間范圍內(nèi)仍可用平穩(wěn)序列來(lái)獲得較為精確的描述,據(jù)此,研究分段局部平穩(wěn)時(shí)序數(shù)據(jù)的辨識(shí)算法有著較為重要的現(xiàn)實(shí)意義.對(duì)分段局部平穩(wěn)時(shí)序數(shù)據(jù)進(jìn)行辨識(shí),其首要的問(wèn)題是找出各局部平穩(wěn)子序列之間的分界點(diǎn),其次是確定這些子序列的模型階次和辨識(shí)參數(shù).

3 非平穩(wěn)時(shí)序數(shù)據(jù)的分段辨識(shí)算法

3.1 算法的設(shè)計(jì)原理

1)分界點(diǎn)的析出方法

對(duì)于某一平穩(wěn)時(shí)序數(shù)據(jù)而言,將它拆分為m段樣本長(zhǎng)度為k的子序列(各段子序列的樣本長(zhǎng)度也可以不相等),則這些子序列的均值、方差和自相關(guān)函數(shù)的樣本參數(shù)估計(jì)值可分別用如下的算式進(jìn)行計(jì)算:

如前所述,這些從平穩(wěn)序列拆分出來(lái)的子序列都應(yīng)具有一致的數(shù)字統(tǒng)計(jì)特征,即任兩個(gè)子序列之間的均值、方差和自相關(guān)函數(shù)等數(shù)字特征不應(yīng)有顯著性的差異.據(jù)此,以上述數(shù)字特征的理論方差為校驗(yàn)統(tǒng)計(jì)量,取顯著水平α=0.05,若式(7)～(9)同時(shí)成立,則可判定u,v兩個(gè)子序列是屬于同一平穩(wěn)序列[18-19].

根據(jù)上述的分析,我們可以得到一種如下的局部平穩(wěn)時(shí)序數(shù)據(jù)分界點(diǎn)的析出算法:按照自左到右的次序往原序列添加Step(Step≥1)個(gè)樣本數(shù)據(jù),并利用式(4)～(6)分別計(jì)算原序列的及添加數(shù)據(jù)后的新序列的同時(shí),利用式(10)～(12)分別計(jì)算出然后根據(jù)式(7)～(9)對(duì)兩個(gè)子序列的均值、方差和自相關(guān)函數(shù)等數(shù)字特征的顯著性進(jìn)行判別,從而決定當(dāng)前所添加的第Step點(diǎn)是否為突變點(diǎn).重復(fù)遍歷處理整個(gè)序列,直至從原序列中析出所有的突變點(diǎn)為止.

從式(10)～(12)易知,對(duì)添加樣本后的序列進(jìn)行相關(guān)的顯著性校驗(yàn)時(shí),均需要顯式調(diào)用均值、方差和自相關(guān)函數(shù)的樣本參數(shù)估計(jì)值,為了更好地提高算法的計(jì)算效能,有必要對(duì)式(4)～(5)進(jìn)行遞推計(jì)算的改進(jìn).

由于:

聯(lián)合式(14),則從式(15)可得到如下的方差遞推計(jì)算表達(dá)式

又由于:

類(lèi)似地,聯(lián)合式(17),則可從式(18)得到如下的自 相關(guān)函數(shù)遞推計(jì)算表達(dá)式:

2)基于AR模型的各分段局部平穩(wěn)時(shí)序數(shù)據(jù)的辨識(shí)

而模型階次 p則可用式(22)所示的 FPE(Final Prediction Error)最小最終預(yù)報(bào)誤差準(zhǔn)則來(lái)確定[21].

3.2 算法的設(shè)計(jì)實(shí)現(xiàn)

綜上所述,可設(shè)計(jì)如下的分段局部平穩(wěn)時(shí)序數(shù)據(jù)的遞推辨識(shí)算法.

步驟1.在Yt中析出各局部平穩(wěn)子時(shí)序數(shù)據(jù)的分界點(diǎn).

1)從Yt的最左端選取L個(gè)右鄰樣本數(shù)據(jù)構(gòu)成Y1t子序列;

2)利 用 式 (4)～(6)計(jì) 算 Y1t對(duì) 應(yīng) 的

3)往Y1t中添加Step個(gè)右鄰樣本數(shù)據(jù)并合并成子序列,利用式(13)、(16)、(19)遞推計(jì)算Y2t對(duì)應(yīng)的y

5)用式(7)～(9)對(duì)第Step點(diǎn)進(jìn)行突變分界點(diǎn)判別,若第Step點(diǎn)為分界點(diǎn),則保存分界點(diǎn)同時(shí)設(shè)置第Step點(diǎn)為Y1t子序列的左起點(diǎn)并跳轉(zhuǎn)(1.2);否則,并跳轉(zhuǎn)(1.3);若Yt遍歷處理完畢后,跳轉(zhuǎn)步驟2.

步驟2.依據(jù)步驟1中所得的各分界點(diǎn),將Yt劃分為數(shù)段局部平穩(wěn)的子序列在此基礎(chǔ)上,基于AR模型對(duì)各段局部平穩(wěn)子序列進(jìn)行辨識(shí),辨識(shí)參數(shù)和模型階次可由式(21)和式(22)得到.

步驟3.輸出有關(guān)計(jì)算結(jié)果并結(jié)束算法.

4 實(shí)驗(yàn)及結(jié)果分析

為了驗(yàn)證上述算法的合理性及有效性,這里將對(duì)具有分段局部平穩(wěn)的時(shí)序數(shù)據(jù)進(jìn)行相關(guān)的辨識(shí)實(shí)驗(yàn),實(shí)驗(yàn)的主要目的是對(duì)比本文算法與現(xiàn)有算法在辨識(shí)精度及計(jì)算效能上的差異.實(shí)驗(yàn)的硬件環(huán)境為惠普ProDesk 490 G2 MT商 用 臺(tái) 式 機(jī) (CPU:i5-45704*3.2GHz;內(nèi)存:4GB DDR3 1600),軟件環(huán)境及開(kāi)發(fā)工具為Windows 8.1+Microsoft Visual C++2010.

4.1 實(shí)驗(yàn)設(shè)計(jì)

實(shí)驗(yàn)所用的模型如式(23)所示,如圖1所示,該模型由3段樣本長(zhǎng)度均為100的局部平穩(wěn)子序列組成.

圖1 由3段局部平穩(wěn)子序列組成的實(shí)驗(yàn)?zāi)Ｐ?/p>

分別用文獻(xiàn)[15]、文獻(xiàn)[16]及本文的算法對(duì)式(23)進(jìn)行相關(guān)辨識(shí),辨識(shí)的任務(wù)為析出各局部平穩(wěn)子序列的分界點(diǎn),并確定這些子序列的模型階次和辨識(shí)參數(shù).

4.2 實(shí)驗(yàn)結(jié)果與討論

如表1所示,文獻(xiàn)[15]、文獻(xiàn)[16]及本文的算法均能正確地從實(shí)驗(yàn)?zāi)Ｐ椭形龀?段局部平穩(wěn)子序列;在分界點(diǎn)的劃分精度問(wèn)題上,文獻(xiàn)[15]的算法在第1、2個(gè)分界點(diǎn)上分別偏離了實(shí)驗(yàn)?zāi)Ｐ?個(gè)和2個(gè)樣本點(diǎn),文獻(xiàn)[16]的算法則分別偏離了4個(gè)和2個(gè)樣本點(diǎn),而本文的算法則分別偏離了8個(gè)和7個(gè)樣本點(diǎn),上述算法所析出分界點(diǎn)的誤差范圍均在10%以?xún)?nèi).據(jù)此,應(yīng)用均值、方差及自相關(guān)函數(shù)等數(shù)字統(tǒng)計(jì)特征對(duì)時(shí)序數(shù)據(jù)進(jìn)行平穩(wěn)性校驗(yàn)是有效的.對(duì)于各局部平穩(wěn)子序列的模型階次而言,上述算法均能正確地識(shí)別出與實(shí)驗(yàn)?zāi)Ｐ拖嘁恢碌碾A次.

表1 三種算法的辨識(shí)結(jié)果

分界點(diǎn)的析出偏差所引起的辨識(shí)精度問(wèn)題,可用式(24)的平均絕對(duì)百分誤差(MAPE,mean absolute percentage)來(lái)進(jìn)行評(píng)價(jià).如表2所示,文獻(xiàn)[15]的算法的辨識(shí)精度最高,其次是文獻(xiàn)[16]的算法,而本文的算法在辨識(shí)精度上略差于上述2種辨識(shí)算法,究其原因是因?yàn)樵谶f推計(jì)算過(guò)程中本文算法存在著固有的積累誤差,從而使得分界點(diǎn)的析出范圍過(guò)大,并最終導(dǎo)致了Burg算法在模型參數(shù)估計(jì)時(shí)出現(xiàn)了一定的偏差.本文的算法的MAPE數(shù)值雖然比文獻(xiàn)[15]和[16]的算法有所增加,但由于各辨識(shí)曲線(xiàn)的MAPE數(shù)值均沒(méi)有超出10%,故本文的算法仍屬于高精度辨識(shí).

表2 三種算法的辨識(shí)精度評(píng)價(jià)

如表3所示,本文算法在辨識(shí)耗時(shí)上有了顯著的提升,其計(jì)算效能較文獻(xiàn)[15]和文獻(xiàn)[16]的算法分別提升了約80%和70%.事實(shí)上,由于在分界點(diǎn)的析出過(guò)程中引入了遞推機(jī)制,故本文算法的計(jì)算耗時(shí)主要是花費(fèi)在各局部平穩(wěn)子序列的模型定階及模型參數(shù)估計(jì)的過(guò)程中;相對(duì)地,文獻(xiàn)[15]和文獻(xiàn)[16]的算法由于需要計(jì)算多維條件的分布概率來(lái)估計(jì)最優(yōu)分界點(diǎn),故它們的計(jì)算耗時(shí)主要是花費(fèi)在分界點(diǎn)的析出過(guò)程中.

表3 三種算法的辨識(shí)耗時(shí)

從上述實(shí)驗(yàn)結(jié)果及分析中易知,本文提出的分段局部平穩(wěn)時(shí)序數(shù)據(jù)的遞推辨識(shí)算法是有效可行的,在保證高精度的辨識(shí)條件下,本文算法在計(jì)算效能取得了顯著的提升.

5 結(jié)語(yǔ)

本文設(shè)計(jì)實(shí)現(xiàn)了一種具有遞推機(jī)制的分段局部平穩(wěn)時(shí)序數(shù)據(jù)的辨識(shí)算法,并就算法的辨識(shí)精度和計(jì)算效能進(jìn)行了深入的分析.下一步的主要工作有,設(shè)法減少遞推過(guò)程中的積累誤差,以便更精確地析出各局部平穩(wěn)時(shí)序數(shù)據(jù)的分界點(diǎn);同時(shí),研究各分段局部平穩(wěn)時(shí)序數(shù)據(jù)的并行辨識(shí)問(wèn)題,從而更好地提升算法的計(jì)算效能.

1張賢達(dá).用高階累量辨識(shí)MA系統(tǒng)的新方法.電子學(xué)報(bào), 1994,22(10):27–33.

2王正明,易東云.含ARMA噪聲系統(tǒng)模型的參數(shù)辨識(shí)方法,控制理論與應(yīng)用,1996,13(8):471–475.

3熊淵博.ARMA模型參數(shù)的分步估計(jì)方法.湖南大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,30(2):12–15.

4周毅,丁鋒.依等價(jià)AR模型階次遞增的自回歸滑動(dòng)平均模型辨識(shí).華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,34(3): 425–431.

5張新廣,李志農(nóng),王心怡,等.一種基于階比域的AR模型盲辨識(shí)算法.振動(dòng)與沖擊,2009,28(5):41–43.

6博克斯,詹金斯,萊因澤爾.時(shí)間序列分析:預(yù)測(cè)與控制.北京:機(jī)械工業(yè)出版社,2011.

7 Matilainen M,Nordhausen K,Oja H.New independent component analysis tools for time series.Statistics& Probability Letters,2015:80–87.

8 Yin Y,Shang PJ,Xia JN.Compositional segmentation of time series in the financial markets.Applied Mathematics and Computation,2015:399–412.

9 Bernas M,P?aczek B.Period-aware local modelling and data selection for time series prediction.Mechanical Systems and Signal Processing,2016:60–77.

10王宏禹,邱天爽,陳喆.非平穩(wěn)隨機(jī)信號(hào)分析與處理.北京:國(guó)防工業(yè)出版社,2008.

11 Spiridonakos MD,Fassois SD.Adaptable functional series TARMA models for non-stationary signal representation and their application to mechanical random vibration modeling.Signal Processing,2014:63–79.

12朱學(xué)鋒,韓寧.基于小波變換的非平穩(wěn)信號(hào)趨勢(shì)項(xiàng)剔除方法.飛行器測(cè)控學(xué)報(bào),2006,25(5):81–85.

13林樹(shù)寬,楊玫,喬建忠,等.一種非線(xiàn)性非平穩(wěn)時(shí)間序列預(yù)測(cè)建模方法.東北大學(xué)學(xué)報(bào)(自然科學(xué)版),2007,28(3): 325–328.

14 DjuricPM,Kay SM.FayeBoundreaux-BartelsG. Segmentation of nonstationary signals.Proc.of the IEEE ICASSP.1992.161–164.

15王文華,王宏禹.分段平穩(wěn)隨機(jī)過(guò)程的參數(shù)估計(jì)方法.電子科學(xué)學(xué)刊,1997,19(5):311–317.

16陳穎,李在銘.一種改進(jìn)的分段平穩(wěn)隨機(jī)過(guò)程的參數(shù)估計(jì)方法.電子與信息學(xué)報(bào),2003,25(6):735–740.

17金連文,韋崗.現(xiàn)代數(shù)字信號(hào)處理簡(jiǎn)明教程.北京:清華大學(xué)出版社,2003.

18楊叔子,吳雅,軒建平,等.時(shí)間序列分析的工程應(yīng)用(上冊(cè)).武漢:華中科技大學(xué)出版社,2007.

19項(xiàng)靜恬.動(dòng)態(tài)數(shù)據(jù)處理-時(shí)間序列分析.北京:氣象出版社,1986.

20黃雄波.時(shí)序數(shù)據(jù)的周期模式發(fā)現(xiàn)算法的遞推改進(jìn).計(jì)算機(jī)技術(shù)與發(fā)展,2016,26(2):47–51.

21克西蓋斯納,沃特斯,哈斯勒.現(xiàn)代時(shí)間序列分析導(dǎo)論.北京:中國(guó)人民大學(xué)出版社,2015.

SegmentedIdentification and RecursiveAlgorithm forNon-StationaryTime Series Data

HUANG Xiong-Bo

(Department of Electronic and Information Engineering,Foshan Professional Technical College,Foshan 528137,China)

In fact,there widely exists a kind of time series data that is non-stationary but can be transformed into several local stationary time series data,the identification problem of the non-stationary time series data is studied,and then this paper proposes a piecewise recursive identification algorithm with mechanism.Based on the definition of stationary time series data,the precipitation algorithm which has the mean variance point,the mutation point and the mutation point of the autocorrelation function,is constructed based on the statistical characteristics such as the mean and variance and autocorrelation function.On this basis,a series of locally stationary sub sequences are identified from the identified non-stationary sequences,and then,the Burg algorithm is applied to the recursive identification of local stationary subsequences.The experimental results show that the new algorithm can divide the boundary points of the local stationary sub sequences with smaller position deviation.At the same time,the calculation efficiency is improved significantly under the condition of high accuracy.

non-stationary time series data;local stationary;segmented identification;recursion;demarcation point; auto-regressive model

廣東省科技計(jì)劃工業(yè)攻關(guān)項(xiàng)目(2011B010200031);佛山職業(yè)技術(shù)學(xué)院校級(jí)重點(diǎn)科研項(xiàng)目(2015KY006)

2016-08-14;收到修改稿時(shí)間:2016-10-19

10.15888/j.cnki.csa.005772