改進的全球衛(wèi)星導(dǎo)航系統(tǒng)陣列抗干擾方法

張彥軍 ，栗　蘋

( 1.北京理工大學(xué)機電工程與控制國家重點實驗室，北京 100081；

2.中北大學(xué)電子測試技術(shù)國家重點實驗室，山西 太原 030051)

?

改進的全球衛(wèi)星導(dǎo)航系統(tǒng)陣列抗干擾方法

張彥軍1,2，栗蘋1

( 1.北京理工大學(xué)機電工程與控制國家重點實驗室，北京 100081；

2.中北大學(xué)電子測試技術(shù)國家重點實驗室，山西 太原 030051)

摘要:針對全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)陣列接收機采用多級維納濾波進行降秩處理時，若在期望方向混入連續(xù)波干擾后會導(dǎo)致性能變差甚至是信號相消的現(xiàn)象，提出了基于前置自適應(yīng)陷波器的組合改進方法。該方法在基于廣義旁瓣相消結(jié)構(gòu)的均勻圓陣上進一步采用Householder變換實現(xiàn)了算法的數(shù)值計算。仿真和試驗表明整個系統(tǒng)構(gòu)架具有很好的抗期望方向連續(xù)波干擾和抗寬帶干擾的性能。

關(guān)鍵詞:全球?qū)Ш叫l(wèi)星系統(tǒng)；自適應(yīng)陷波器；陣列信號處理；多級維納濾波

0引言

隨著導(dǎo)航衛(wèi)星技術(shù)的發(fā)展與應(yīng)用，其地面接收機面臨的窄帶或?qū)拵褐剖礁蓴_問題也日益突出。基于陣列信號處理技術(shù)的接收機可以通過調(diào)整各接收陣元的加權(quán)復(fù)系數(shù)來自適應(yīng)地消除這些干擾，但是由于涉及到大規(guī)模矩陣求逆或求解矩陣特征值等問題，在實際應(yīng)用中受到一定的限制。

多級維納濾波器(multistage Wiener filter, MWF)作為近些年一種新的降秩方法[1]，既不需要相關(guān)矩陣的特征結(jié)構(gòu)，也無需矩陣求逆運算，與其他基于特征值的算法(主分量法、交叉譜法)相比，在降秩更多的同時仍然能夠保持很好的性能[2]，因此在空域濾波[3]和空時濾波等陣列信號處理領(lǐng)域得到了廣泛應(yīng)用[4]。文獻[5]將MWF引入到廣義旁瓣相消的算法結(jié)構(gòu)中，形成了線性約束多級維納濾波器(linearly constrained MWF,LC-MWF)。但是，廣義旁瓣相消器要求期望信號只能從非自適應(yīng)的主支路通過，而自適應(yīng)的輔助支路中僅含有干擾和噪聲分量。在高信噪比的情況下，將有一部分期望信號泄露到輔助支路而出現(xiàn)信號相消的現(xiàn)象[6]。

對于衛(wèi)星信號而言，它到達接收機輸入端時信號極其微弱，以GPS衛(wèi)星信號的C/A碼為例，其典型值約為-160 dBw(環(huán)境溫度：290 K);低于熱噪聲約20 dB，因此不會出現(xiàn)高信噪比的情況。但是如果在期望方向出現(xiàn)大能量的干擾就會等效成高的信噪比，所以在進行LC-MWF算法之前必須予以剔出，文獻[7]討論了GPS接收機在期望方向產(chǎn)生非平穩(wěn)的沖激毛刺干擾下處理方法。此外，連續(xù)波干擾(continuous wave interference，CWI)也是常見的一種干擾[8],高頻TV發(fā)射器、無線電導(dǎo)航中使用的高頻全向信標等系統(tǒng)都可以產(chǎn)生CWI諧波[9]。如果在期望方向或靠近期望方向上混有大能量的CWI，這時可以前置自適應(yīng)窄帶IIR陷波器進行預(yù)消除。同時，它也保證了后續(xù)的空域濾波或空時濾波只需處理寬帶干擾，提高了其抗干擾的自由度和性能。

本文對實際受CW干擾的導(dǎo)航數(shù)據(jù)進行上述方法處理，并通過捕獲結(jié)果驗證了該改進方法的可行性。隨后在多級維納濾波器中引入householder變換(Householder MWF,HMWF)，進一步降低了陣列算法的運算量并提高了算法的穩(wěn)健性[10]，形成了衛(wèi)星信號比較有效的陣列抗干擾方法。

1期望方向含有CWI的LC-MWF

假設(shè)有16均勻線陣，陣元間距為導(dǎo)航載波信號的半波長，期望信號為帶寬20 M的導(dǎo)航信號從方向角0°射入，信噪比是-30 dB，噪聲為帶寬20 M的熱噪聲；3個同樣帶寬的寬帶高斯干擾方向角分別是-30°、20°和60°，其信干比分別是-48 dB、-53 dB和-58 dB。

在輸出最大信干噪比的條件下

(1)

求得

(2)

其中Wopt為所求系數(shù)向量W的最優(yōu)權(quán)值，Ri+n為干擾加噪聲的自相關(guān)矩陣，α(θ0)為歸一化的期望信號導(dǎo)向矢量，Rs為信號的自相關(guān)矩陣。

為了避免對協(xié)方差矩陣Ri+n求逆，構(gòu)建如圖1所示的基于廣義旁瓣相消結(jié)構(gòu)的多級維納濾波器，即LC-MWF。

圖1　線性約束多級維納濾波器Fig.1　Linearly constrained multistage Wiener filter

圖1中X(k)為16陣元的輸入矢量，B為期望矢量α(θ0)的阻塞矩陣，T為前向分解的轉(zhuǎn)換矩陣，Wd為后向合成的系數(shù)矢量。

則有

WMWF=THWd

(3)

WLC-WMF=α(θ0)-BWMWF

(4)

如圖2實線所示，假設(shè)LC-MWF為全秩轉(zhuǎn)換，在導(dǎo)航信號期望方向不含干擾的理想條件下，當取Wopt=WLC-WMF時，在3個干擾方向獲得極深的零陷，同時保持期望方向(0°)的功率值最大；當從期望方向混入信干比為-60 dB的強連續(xù)波干擾信號時，如圖2虛線所示，發(fā)現(xiàn)在期望方向出現(xiàn)了信號相消的現(xiàn)象，同時在3個干擾的方向上也出現(xiàn)了零陷變淺的現(xiàn)象，導(dǎo)致了整個抗干擾性能的下降。

圖2　線性約束多級維納濾波器的期望信號相消Fig.2　Desired signal cancellation of LC-MWF

為了防止導(dǎo)航方向的干擾對LC-MWF造成信噪比的嚴重下降，必須對輸入信號X(k)進行濾波預(yù)處理，剔除CW干擾，同時保證處理后的導(dǎo)航信號不能失真嚴重，影響后續(xù)的捕獲性能。顯然，F(xiàn)IR濾波器具有線性相位的優(yōu)點，并且沒有反饋回路，無需考慮穩(wěn)定性的問題。但同時要達到很好的截止特性需要的階次較高，會帶來計算量和存儲量的增加。與FIR濾波器相比，低階的IIR濾波器可以達到同樣的濾波截止特性，能夠滿足導(dǎo)航信號后續(xù)信號處理的要求。

2基于前置自適應(yīng)IIR格型陷波器的改進方法

2.1基于互相關(guān)的自適應(yīng)IIR格型陷波器

利用IIR自適應(yīng)格型陷波器去除CW干擾是直接序列擴頻通信中的一種有效辦法。如圖3 所示，u(n)為輸入信號，經(jīng)過一個陷波器輸出y(n)，經(jīng)過一個點頻帶通濾波器輸出x(n)，其傳遞函數(shù)分別為：

(5)

(6)

并且有：α=(1-tan(B/2))/(1+tan(B/2)),β=cos(ωN),ωN∈[0,π]。其中，α參數(shù)控制-3 dB陷波帶寬B，β參數(shù)控制陷波頻率ωN。

文獻[11]主要是通過FFT將采樣的時域信號轉(zhuǎn)換到頻域然后尋找幅值最大的頻率值，通過該值求解得β，然后再利用陷波帶寬α參數(shù)的求解公式獲得α。由于該方法需要FFT運算，并且為防止頻率泄露需要加窗處理，所以增加了硬件資源和時間的復(fù)雜度。

而在空時自適應(yīng)濾波處理中，陷波作為一種前置的輔助方法，不宜耗費過多的資源和時間。因此在本文中主要采用了基于時域互相關(guān)值的迭代求解實時獲得β值，從而避開FFT運算。

設(shè)輸入信號u(n)的數(shù)學(xué)表達式是：

u(n)=s(n)+w(n)+j(n)

(7)

其中s(n)為衛(wèi)星導(dǎo)航信號，w(n)為高斯加性白噪聲，j(n)為連續(xù)波干擾，取單頻正弦波為其表達式：j(n)=Acos(ωjn+φ0),則輸入信號u(n)經(jīng)過傳遞函數(shù)HB(z)，濾除了s(n)+w(n)成分，經(jīng)過傳遞函數(shù)HN(z)，濾除了j(n)成分，因此x(n)和y(n)的互相關(guān)值可以看作干擾和信號與噪聲的互相關(guān)值的即時映射，利用它可以調(diào)整陷波頻率參數(shù)β。

定義φ[n]為x[n]干擾功率瞬時估計值，由于在采集和數(shù)字量化過程中存在高頻擾動，所以對干擾功率進行一階遞歸數(shù)字濾波，即φ[n]=ρφ[n-1]+(1-ρ)x2[n]，其中ρ是遺忘因子。由文獻[12]可知，ρ越接近于1，越有利于抑制高頻區(qū)的噪聲功率，所以通常取0<ρ≈1。為了使β[n]初始時快速收斂而在逼近穩(wěn)定值時慢速收斂，令參數(shù)λ[n]與干擾功率φ[n]成反比，并與調(diào)整收斂率的經(jīng)驗系數(shù)μ相乘以保證穩(wěn)定收斂，則有λ[n]=μ/φ[n]作為時變步長參數(shù)。這樣自適應(yīng)求解的迭代公式如下[13]：

β[n+1]=β[n]-λ[n]×(y[n]x[n])

(8)

2.2圓陣條件下的線性約束Householder多級維納濾波

通過前置的自適應(yīng)陷波器，可以判斷是否有CW干擾并可以將之有效地去除，保證了進入后續(xù)的空域或空時自適應(yīng)多級維納濾波能夠處理相對理想的數(shù)據(jù)。

為了使多級維納濾波器有更好的降秩性能，文獻[14]和文獻[15]分別提出了兩種不同的阻塞矩陣構(gòu)成方法，但都是分兩步進行，即先求每一級的期望矢量，再構(gòu)造相應(yīng)的阻塞矩陣。為了進一步減低陣列計算的復(fù)雜度和提高陣列算法的魯棒性，在圖1的基礎(chǔ)上引入HMWF算法，形成圖4的LC-HMWF。

圖4　線性約束Householder多級維納濾波器Fig.4　Linearly constrained Householder MWF

1)迭代初始化：d0(k)=a(θ0)X(k)，X0(k)=B0X(k)

2)前向分解流程：

fori=1,2,…,D

vi=pi-kiδie1(e1是單位矢量)；

3)后向迭代流程：

εD=dD

fori=D,D-1,…,1

文獻[10]討論了基于均勻線陣的空時Householder多級維納濾波方法。但由于均勻線陣只能提供180°無模糊的方位角能力，而均勻圓陣作為一個平面陣不僅具有360°的全向方位角覆蓋能力，同時還有俯仰角指向能力，有二維空間波束形成能力。因此為了進一步拓寬多級維納濾波器應(yīng)用的廣泛性，本文對之進行了圓陣的推導(dǎo)應(yīng)用。

假設(shè)有M陣元組成均勻圓陣，陣元半徑為導(dǎo)航載波信號的半波長，當每個陣元的輸入信號經(jīng)過前述的IIR自適應(yīng)預(yù)陷波后，進入后續(xù)的基于HMWF的二維空時濾波結(jié)構(gòu)。信號的空時輸入模型為X(k)=[x11,x12,…,x1N,x21,…,xM1,xM2,…,xMN](M為陣元數(shù)，N為延遲環(huán)節(jié)數(shù))，每個陣元后的延遲時間τ<1/B，B為信號的帶寬，此處為20 M。

去除CW干擾后的衛(wèi)星期望信號是一個寬帶信號，因此其導(dǎo)向矢量是帶內(nèi)各個頻點的空時函數(shù)，其頻點fi所對應(yīng)的導(dǎo)向矢量為ai(θ0,φ0)=Ss?St,其中?為Kronecker product符號，St時間導(dǎo)向矢量和Ss空間導(dǎo)量矢量分別表示為Ss=[1,ejωs,…,ej(M-1)ωs]和St=[1,ejωt,…,ej(M-1)ωt]，其中ωs和ωt分別表示歸一化的空間頻率和時間頻率；θ0和φ0分別是期望方向的俯仰角和方位角。將 (θ0,φ0)方向k個頻點的ai(θ0,φ0)求均值得到該方向的a(θ0,φ0)將之代入圖4中，就可以計算得基于LC-HMWF算法的空時自適應(yīng)濾波的權(quán)系數(shù)。

3數(shù)值仿真與試驗

3.1IIR格型陷波器試驗數(shù)據(jù)處理與性能分析

在陣列前端采用IIR格型陷波器，既要保證能夠?qū)W干擾濾除，同時又不能夠影響處理后數(shù)據(jù)與本地碼的相關(guān)性，從而導(dǎo)致捕獲性能變差。因此控制-3 dB帶寬參數(shù)的α(0<α<1)至關(guān)重要。如圖5所示，α越大，陷波器的陷波帶寬越小;當α=1時，IIR格型陷波器就會變成一個全通濾波器。所以，通常α的設(shè)置值接近于1,以保證在濾除干擾的前提下，不能夠?qū)罄m(xù)的數(shù)據(jù)產(chǎn)生很大的非線性誤差。

當GNSS信號從天線進入前端采樣后，由帶通采樣定理可知其表達式u(n)為一組中頻帶的實數(shù)數(shù)據(jù)序列，此時可以將數(shù)據(jù)正交下變頻后在基帶復(fù)數(shù)域進行自適應(yīng)陷波[16]，也可以直接在中頻域處理。考慮到硬件資源的開銷，在中頻域處理更為合適。

圖6所示的為一組陣元采集的導(dǎo)航數(shù)據(jù)的頻譜圖，可以看到載波中心頻率附近處混雜有一條比較大的譜線幅值，在此連續(xù)波的干擾下采用匹配濾波器相關(guān)捕獲的圖形如圖7，其相關(guān)峰值混雜在雜亂的相關(guān)值中，難有可靠的捕獲。同時更為嚴重的后果是當數(shù)據(jù)進入陣列信號處理算法后會造成如圖2所示的期望信號相消現(xiàn)象，從而導(dǎo)致抗干擾算法的性能降低。

為了減少濾波后信號的失真，取陷波帶寬參數(shù)α=0.98，在此條件下由仿真可知陷波器幅頻曲線中的陷波深度約為-25dB。因此在干噪比為0～-25dB的范圍中設(shè)置不同的連續(xù)波干擾功率來確定經(jīng)驗收斂系數(shù)μ，以確保式(8)能夠在此范圍內(nèi)收斂于干擾頻率，從而獲得μ=2-6，同時取ρ=1-2-6≈1。在此條件下，經(jīng)表達式(8)和(6)進行迭代濾波后獲得的穩(wěn)定數(shù)據(jù)y(n)的頻譜如圖8所示，頻譜中心的連續(xù)波干擾譜線被剔出。對該數(shù)據(jù)進行匹配濾波器相關(guān)捕獲后相關(guān)峰值如圖9所示，其捕獲的矩陣均峰比值可以達到約為4.35，超過了該接收機實際捕獲時預(yù)置的經(jīng)驗門限值3.5，能夠為后續(xù)的陣列信號處理提供可靠的處理數(shù)據(jù)。文獻[13]指出，IIR格型陷波器的-3dB帶寬在小于30kHz時，其固有的非線性相位響應(yīng)不會對后續(xù)的捕獲和跟蹤產(chǎn)生影響。

3.2LC-HMWF算法數(shù)據(jù)處理與性能分析

假設(shè)有七均勻圓陣，一個陣元在圓心，其余六個均布在陣元半徑為載波半波長的圓周上，每個陣元的延時環(huán)節(jié)數(shù)為8。在CW干擾已經(jīng)被抑制的條件下，只有兩個寬帶干擾的信干比分別為-48dB和-53dB，熱噪聲信噪比是-30dB。期望信號俯仰角和方位角設(shè)為30°和50°，其余2個干擾的方向分別是：(50°，100°)，(70°，300°)。采用LC-HMWF算法，設(shè)置定點采樣精度為16bit，在秩rank=15的條件下，不失一般性，觀測期望信號俯仰角30°時方位角與頻率的聯(lián)合方向圖和干擾信號方位角100°時俯仰角與頻率的聯(lián)合方向圖分別如圖10和圖11所示。可見算法既在期望方向形成了主瓣，同時在干擾方向形成了零陷。圖12顯示了LC-HMWF算法的秩與輸出信噪比的關(guān)系，在秩等于15的時候就能夠取得全秩的性能。圖13則表明了最終取得的信噪比能夠保證信號獲得很好的捕獲性能，其中橫軸對碼相位的取值范圍進行了按比例縮小。

在實際試驗中，分別有兩個干擾源對接收機七均勻圓陣天線陣實施干擾，一個為連續(xù)波干擾(頻點值為衛(wèi)星載波頻率，信干比約-63dB)；另一個為寬帶BPSK干擾(帶寬20MHz，覆蓋衛(wèi)星信號的全部帶寬，信干比約-80dB)。測試過程中，固定寬帶干擾天線的位置不變，然后移動連續(xù)波干擾天線，使其與一顆可接收的衛(wèi)星波達方向基本一致。在經(jīng)過具有前置自適應(yīng)陷波的多級維納濾波后，有7顆星被成功捕獲。

圖5　自適應(yīng)IIR格型陷波器幅頻曲線Fig.5　Amplitude-frequency curve of adaptive IIR lattice notch filter

圖6　CW干擾下的數(shù)據(jù)頻譜圖Fig.6　Spectrogram under CW jamming

圖7　CW干擾下的二維捕獲圖Fig.7　Two-dimensional capture under CW jamming

圖8　自適應(yīng)濾波后的數(shù)據(jù)頻譜圖Fig.8　Spectrogram after adaptive filtering

圖9　自適應(yīng)濾波后的二維捕獲圖Fig.9　Two-dimensional capture after adaptive filtering

圖10　方位角-頻率方向圖Fig.10　The azimuth-frequency pattern

圖11　俯仰角-頻率方向圖Fig.11　The pitch-frequency pattern

圖12　LC-HMWF的秩-輸出信噪比關(guān)系圖Fig.12　LC-HMWF rank-output SJNR

圖13　陣列輸出信號的碼相關(guān)圖Fig.13　Code correlation output of array

4結(jié)論

本文提出了改進的GNSS陣列抗干擾方法。該方法針對多級維納濾波器在衛(wèi)星導(dǎo)航信號期望方向混有大能量的變化連續(xù)波干擾時出現(xiàn)的信號相消現(xiàn)象，在不影響信號的捕獲性能條件下，對采樣數(shù)據(jù)進行了前置自適應(yīng)IIR陷波器的預(yù)處理方法來保證降秩算法的有效性，同時提高了抗寬帶干擾的空間自由度。仿真和試驗表明，預(yù)處理后的數(shù)據(jù)在引入LC-HMWF的圓陣空時抗干擾算法結(jié)構(gòu)中抑制了寬帶干擾，在降秩的條件下實現(xiàn)了對信噪比的提升和信號的捕獲，比較有效地實現(xiàn)了一種衛(wèi)星信號的陣列抗干擾方法。

參考文獻:

[1]GoldsteinJS,ReedIS,ScharfLL.AmultistagerepresentationoftheWienerfilterbasedonorthogonalprojections[J].IEEETransactionsonInformationTheory, 1998, 44(7): 2943-2959.

[2]RosenbergLuke,GrayDouglasA.ConstrainedFast-TimeSTAPforInterferenceSuppressioninMultichannelSAR[J].IEEETransactionsonAerospaceAndElectronicSystems, 2013, 49(3): 1792-1805.

[3]SongNuan,LamareRodrigoCde,HaardtMartin,etal.AdaptiveWidelyLinearReduced-RankInterferenceSuppressionBasedontheMultistageWienerFilter[J].IEEETransactionsonsignalprocessing, 2012, 60(8): 4003-4016.

[4]HuChiaChang,LiaoYungCheng.EfcientMIMODS-UWBDownlinkCh-LevelMMSEEqualizationUsingSubbandAdaptiveInterferenceMitigationTechniques[J].IEEETransactionsonVehicularTechnology, 2012, 61(8): 3783-3790.

[5]丁前軍, 王永良, 張永順,等.自適應(yīng)陣列處理典型算法的內(nèi)在關(guān)系與統(tǒng)一框架[J]. 系統(tǒng)工程與電子技術(shù), 2006, 28(2): 163-174.

[6]張小飛, 徐錕, 徐大專. 一種改進的廣義旁瓣相消的波束形成算法[J]. 中國空間科學(xué)技術(shù)，2006, 26(5):6-9.

[7]黃慶東，張林讓，盧光躍，申東一種改進HMSWF導(dǎo)航接收機空時抗干擾方法[J]. 電子與信息學(xué)報，2011,33(10): 2526-2530.

[8]BorioD,O’DriscollC,FortunyJ.GNSSjammers:Effectsandcountermeasures[C]// 6thESAWorkshopSatelliteNavigationTechnologiesandEuropeanWorkshoponGNSSSignalsandSignalProcessing.NAVITEC,2012:1-7.

[9]VarshneyNeeraj,JainRC.AnAdaptiveNotchFilterforNarrowBandInterferenceRemoval[C]//IEEE, 2013NationalConferenceonCommunications(NCC), 2013.

[10]WernerS,WithM,KoivunenV.Householdermultistagewienerfilterforspace-timenavigationreceivers[J].IEEETransactionsonAerospaceandElectronicSystems, 2007, 43(3):975-988.

[11]龔文飛, 吳嗣亮, 李加琪. 直擴系統(tǒng)中IIR格型濾波器抑制窄帶干擾新方法與性能分析[J]. 電子與信息學(xué)報，2010,32(10):2473-2478.

[12]盛壽麟. 運動系統(tǒng)中的一階遞歸數(shù)字濾波器[J].西安交通大學(xué)學(xué)報,1989,23(1):115-120.

[13]ChienYingRen.DesignofGPSAnti-JammingSystemsUsingAdaptiveNotchFilters[J].IEEESystemsJournal, 2013: 1-10.

[14]RicksDC,GoldsteinJS.Efficientarchitecturesforimplementingadaptivealgorithms[C]//Proceedingsofthe2000AntennaApplicationsSymposium.Illinois:IEEEPress,2000,29-41.

[15]丁前軍, 王永良, 張永順. 自適應(yīng)陣列中多級維納濾波器的有效實現(xiàn)算法[J]. 電子與信息學(xué)報, 2006, 28(5): 936-940.

[16]XiongHailiang,ZhangWensheng,DuZhengfeng,etal.Front-EndNarrowbandInterferenceMitigationforDS-UWBReceiver[J].IEEETransactionsonWirelessCommunications, 2013, 12(9):4328-4337.

An Improved Anti-jamming Method for Global Navigation Satellite System Arrays

ZHANG Yanjun1,2, LI Ping1

(1.National Key Laboratory of Mechatronic Engineering and Control, Beijing Institute of Technology, Beijing 100081,China；2.National Key Laboratory for Electronic Measurement Technology, North University of China, Taiyuan 030051,China)

Abstract:When the global navigation satellite system array receiver uses the multistage wiener filter technology to reduce the rank of covariance matrix , it may resulted in poor performance and even the signal cancellation phenomenon on the condition that the continuous wave interference arrives from the expected direction. To resolve the issue, an improved joint method is proposed based on the front-end adaptive notch filter was proposed. The method further adopted the householder transformation to achieve the numerical calculation based on the generalize side lobe canceller for the uniform circular array. The simulation data and experimental data indicated the proposed architecture had good rejection capabilities for the continuo continuous wave interference from the expected direction and wideband interference.

Key words:global navigation satellite system; multistage wiener filter; adaptive notch filter; array signal processing

中圖分類號:TN911.4

文獻標志碼:A

文章編號：1008-1194(2016)01-0007-06

作者簡介:張彥軍(1975—)，男，山西五臺人，博士研究生，副教授，研究方向：微弱信號抗干擾技術(shù)，陳列信號處理。E-mail: zhangyanjun@nuc.edu.cn。

基金項目:國家自然科學(xué)基金項目資助(51407010)

*收稿日期:2015-09-12