田　野，練秋生，徐　鶴

(燕山大學(xué)信息科學(xué)與工程學(xué)院，河北秦皇島 066004)

基于稀疏信號(hào)重構(gòu)的DOA和極化角度估計(jì)算法

田野，練秋生，徐鶴

(燕山大學(xué)信息科學(xué)與工程學(xué)院，河北秦皇島 066004)

現(xiàn)有的波達(dá)方向(Direction Of Arrival,DOA)和極化參數(shù)估計(jì)方法大多基于子空間理論.本文從稀疏信號(hào)重構(gòu)角度出發(fā)，提出了一種新的DOA和極化角度估計(jì)算法.該算法首先構(gòu)建一個(gè)只包含DOA信息的累積量矩陣模型，然后基于加權(quán)l(xiāng)1范數(shù)最小化獲得DOA估計(jì).在DOA估計(jì)的基礎(chǔ)上，進(jìn)一步通過求和平均運(yùn)算構(gòu)建三個(gè)包含不同極化信息的累積量向量模型，利用Zhang懲罰進(jìn)行稀疏性約束，獲得近似無(wú)偏的極化角度估計(jì).闡述了如何利用極化信息來(lái)區(qū)分兩個(gè)入射角度一樣的信源信號(hào).計(jì)算機(jī)仿真結(jié)果驗(yàn)證了所提算法的有效性.

DOA和極化角度估計(jì)；稀疏信號(hào)重構(gòu)；加權(quán)l(xiāng)1范數(shù)；Zhang懲罰

1　引言

波達(dá)方向(Direction Of Arrival,DOA)估計(jì)是陣列信號(hào)處理領(lǐng)域的核心研究?jī)?nèi)容之一，它在智能天線、無(wú)線通信、雷達(dá)、聲吶、電子偵查及地震勘探等領(lǐng)域有著廣泛的應(yīng)用前景[1]，受到了國(guó)內(nèi)外學(xué)者的廣泛關(guān)注，各種適用于DOA估計(jì)的算法被相繼提出[2～5].這些算法均基于標(biāo)量傳感器陣列，它們不能充分利用EM信源信號(hào)的極化信息.

在無(wú)線通信和雷達(dá)領(lǐng)域，極化敏感陣列已被證明可以獲得比標(biāo)量傳感器更多的優(yōu)勢(shì)[6,7].鑒于這一事實(shí)，一些適用于DOA和極化參數(shù)估計(jì)的算法也被提出，包括矢量叉積算法[8]，極化類ESPRIT算法[9～11]，極化四元數(shù)MUSIC算法[12,13]，極化雙模MUSIC算法[14]，極化MODE算法[15]等.上述算法大都基于子空間理論，由于理論框架的限制其分辨率性能有限，且不能充分利用信號(hào)的極化信息.

近年來(lái)，稀疏信號(hào)重構(gòu)作為一種新的理論框架逐漸被引入到DOA估計(jì)領(lǐng)域.目前已出現(xiàn)的稀疏信號(hào)重構(gòu)類DOA估計(jì)算法主要包括l1范數(shù)懲罰算法[16]、加權(quán)l(xiāng)1范數(shù)懲罰算法[17,18]、l0范數(shù)逼近算法[19,20]以及迭代重加權(quán)最小范數(shù)算法[21]等.和傳統(tǒng)子空間類DOA估計(jì)算法相比，稀疏信號(hào)重構(gòu)類算法具有分辨率高、噪聲魯棒性好、無(wú)需精確的初始條件等優(yōu)勢(shì).然而值得注意的是，現(xiàn)有的稀疏重構(gòu)類DOA估計(jì)算法大都存在估計(jì)偏的問題[22]，這在一定程度上影響了參數(shù)估計(jì)性能，并嚴(yán)重制約其向極化敏感陣列的拓展.

為了充分利用EM信號(hào)的極化信息并克服現(xiàn)有稀疏信號(hào)重構(gòu)類DOA估計(jì)算法存在的共性問題，本文借助COLD(Concentered Orthogonal Loop and Dipole)陣列，提出了一種新的基于稀疏信號(hào)重構(gòu)的DOA和極化角度估計(jì)算法.該算法步驟如下：(1)構(gòu)建累積量矩陣并基于加權(quán)l(xiāng)1范數(shù)最小化獲得DOA估計(jì)；(2)通過求和平均運(yùn)算構(gòu)建三個(gè)特殊的累積量向量，利用Zhang懲罰進(jìn)行稀疏性約束獲得極化角度估計(jì)；(3)利用估計(jì)的極化信息判斷是否有兩個(gè)信源從相同的方向入射到陣列.最后通過仿真實(shí)驗(yàn)驗(yàn)證了所提算法的有效性.

2　極化遠(yuǎn)場(chǎng)源模型

考慮K個(gè)完全極化遠(yuǎn)場(chǎng)窄帶信號(hào)入射到一個(gè)由L=2M+1個(gè)COLD陣元組成的均勻線陣上.如圖1所示，其第m(m=-M,…,M)個(gè)陣元以間距d均勻分布在y軸上.偶極子平行于z軸，電磁環(huán)平行于x-y平面.橢圓極化方式下的電場(chǎng)輸出可表示為

E=Eθθ+Eφφ

(1)

式(1)中，Eθ代表水平分量，Eφ代表垂直分量，θ和φ分別代表沿方位角θ和俯仰角φ的球面單位向量.對(duì)于給定的信號(hào)極化，電場(chǎng)分量可進(jìn)一步表示為

Eθ=E0cos(γ),

Eφ=E0sin(γ)ejη

(2)

式(2)中，E0代表非零復(fù)信號(hào)幅度，γ∈[0,π/2)和η∈[-π,π]分別代表信號(hào)的極化角度和極化相位差.

基于式(1)和式(2)，得到入射信號(hào)的電場(chǎng)和磁場(chǎng)分別為

E=Eθθ+Eφφ=E0[cos(γ)θ+sin(γ)ejηφ]

(3)

H=Eθφ-Eφθ=E0[cos(γ)φ-sin(γ)ejηθ]

(4)

為方便且不失一般性，假設(shè)天線和入射信號(hào)是共面的，即φ=90°，則式(3)和式(4)可重新表示為

E=E0[cos(γ)(-sin(θ)x+cos(θ)y)-sin(γ)ejηz]

(5)

H=E0[sin(γ)ejη(sin(θ)x-cos(θ)y)-cos(γ)z]

(6)

其中，x，y和z分別代表沿x，y，z方向的單位向量.

(7)

(8)

以第0個(gè)陣元為相位參考點(diǎn)，陣列輸出可表示為

u[l](t)=As[l](t)+n[l](t)

(9)

u[d](t)=As[d](t)+n[d](t)

(10)

其中

(11)

(12)

s[l](t)=-[s1(t)cos(γ1),…,sK(t)cos(γK)]T

(13)

s[d](t)=-[s1(t)sin(γ1)ejη1,…,sK(t)sin(γK)ejηK]T

(14)

(15)

(16)

其中，上標(biāo)T代表轉(zhuǎn)置操作，A=[a(θ1),…,a(θK)]代表L×K的陣列流型矩陣，其第k列代表第k個(gè)信號(hào)的導(dǎo)向矢量

a(θk)=[e-jMωk,…,1,…,ejMωk]T

(17)

為保證參數(shù)估計(jì)的唯一性，本文做如下假設(shè)：

[A1]信源信號(hào){s1(t),…,sK(t)}為窄帶、統(tǒng)計(jì)獨(dú)立過程且其四階累積量非零；

[A2]噪聲分量n[l](t)和n[d](t)為加性統(tǒng)計(jì)獨(dú)立高斯過程，且與信源信號(hào)不相關(guān)；

[A3]為避免相位模糊，陣元間距和信源數(shù)分別滿足d≤λ/2，K

3　DOA和極化角度估計(jì)算法

3.1DOA估計(jì)

為抑制高斯噪聲并保證估計(jì)精度，本文基于陣列輸出數(shù)據(jù)分別構(gòu)建L×L3的自極化和互極化累積量矩陣，表示為

(18)

(19)

(20)

(21)

為獲得更好的估計(jì)結(jié)果，定義

(22)

式(22)說明信源信號(hào)在累積量域已被完全接收，這直接保證了算法的DOA估計(jì)精度.同時(shí)由于矩陣C只包含信源的DOA信息，因此DOA估計(jì)性能將不受極化參數(shù)的影響.

矩陣C的列向量cp,q,n可表示為

(23)

式(23)中，sp,q,n代表K×1的向量

(24)

(25)

(26)

式(26)中，X是包含所有列向量xp,q,n的矩陣.假定信源數(shù)K已知或已通過AIC或MDL準(zhǔn)則準(zhǔn)確估計(jì)[23]，則可通過對(duì)矩陣C進(jìn)行奇異值分解(Singular Value Decomposition,SVD)來(lái)降低計(jì)算復(fù)雜度，得到

(27)

其中，C=ULVH，CSV=CVDK，XSV=XVDK，DK=[IK,0]T，IK代表K×K的單位矩陣，0代表K×(L3-K)的零值矩陣.上標(biāo)H代表共軛轉(zhuǎn)置操作.

(28)

(29)

(30)

(31)

(32)

(33)

和獲得式(22)類似，定義

(34)

式(34)的矩陣形式為

Cw=AC4sAH

(35)

(36)

(37)

優(yōu)化問題(37)可直接通過內(nèi)點(diǎn)法二階錐規(guī)劃(SOCP)[24]進(jìn)行求解，其標(biāo)準(zhǔn)形式為

s.t.

(38)

3.2極化角度估計(jì)

為獲得極化角度估計(jì)，本文構(gòu)建三個(gè)特殊的累積量向量，分別表示為式(39)～(41).

(39)

(40)

(41)

向量c1到c3的稀疏表示為

(42)

(43)

(44)

(45)

(46)

(47)

(48)

證明當(dāng)只有一個(gè)信源信號(hào)由方向θ入射到陣列時(shí)，有

ρk=cos4(γ)sin4(γ)-(cos2(γ)sin2(γ))2=0

(49)

當(dāng)有兩個(gè)信源信號(hào)由相同的方向θ入射到陣列時(shí)，有式(50).

ρk=[(cos4(γ1)+cos4(γ2))(sin4(γ1)+sin4(γ2))

-(cos2(γ1)sin2(γ1)+cos2(γ2)sin2(γ2))2]

=cos4(γ1)sin4(γ2)+cos4(γ2)sin4(γ1)

-2cos2(γ1)sin2(γ1)cos2(γ2)sin2(γ2)

=(cos2(γ1)sin2(γ2)-cos2(γ2)sin2(γ1))2

=sin2(γ2+γ1)sin2(γ2-γ1)

(50)

由于γ∈[0,π/2)，并且γ1和γ2不同，則γ1+γ2∈(0,π)，γ1-γ2∈(-π/2,0)∪(0,π/2)，因此1>ρk>0，并且ρk隨著γ1→π/2,γ2→0或γ2→π/2,γ1→0而趨于最大值1.證畢.

注1：定理1說明當(dāng)極化角度不同時(shí)，本文所提算法可有效地區(qū)分兩個(gè)入射角度相同的信源信號(hào)，并且區(qū)分性能隨著γ1→π/2,γ2→0或γ2→π/2,γ1→0而變好.

注2：本文極化參數(shù)估計(jì)的目的是為了在極化域區(qū)分兩個(gè)入射角度相同的信源信號(hào)，由于在實(shí)施過程中只依靠極化角度即可完成區(qū)分(見定理1)，故本文未進(jìn)行極化相位差的估計(jì).

3.3算法步驟

本文所提算法步驟可總結(jié)如下：

步驟1構(gòu)建累積量矩陣C，Cw及累積量向量ck；

步驟3對(duì)矩陣Cw進(jìn)行奇異值分解，并利用信號(hào)子空間與噪聲子空間的正交性構(gòu)建權(quán)值矩陣w；

步驟4基于加權(quán)l(xiāng)1范數(shù)最小化方法和SOCP獲得DOA估計(jì)；

步驟5借鑒Zhang懲罰思想，在LASSO估計(jì)初值的基礎(chǔ)上，利用式(45)進(jìn)行稀疏信號(hào)重構(gòu)獲得極化角度估計(jì)；

步驟6通過ρk的值判斷是否有兩個(gè)信源信號(hào)由相同的方向入射到陣列.

(51)

4　仿真與分析

本章通過仿真實(shí)驗(yàn)驗(yàn)證所提算法的性能，并與交叉偶極子陣下最具代表性的ESPRIT算法和COLD陣下的MUSIC算法以及MODE算法進(jìn)行比較.仿真中采用由7個(gè)陣元組成的均勻線陣，陣元間距d=λ/2，入射信號(hào)功率相等，建模為ejζt，其中相位ζt均勻分布在[0,2π].除實(shí)驗(yàn)3外，噪聲均假定為高斯白噪聲.以1°間隔先對(duì)-90°到90°的空間進(jìn)行粗網(wǎng)格劃分，進(jìn)而在估計(jì)出的角度周圍逐步細(xì)化網(wǎng)格.DOA估計(jì)的均方根誤差(RMSE)通過500次獨(dú)立的Monte Carlo仿真實(shí)驗(yàn)獲得，定義為

(52)

實(shí)驗(yàn)1角度分辨率性能比較

為展示本文所提算法的角度分辨率性能，考慮兩個(gè)空間間隔很近的信源信號(hào)入射情況，入射角度分別為θ1=25°，θ2=30°，相對(duì)應(yīng)的極化參數(shù)(γ,η)分別為(20°,0°)和(60°,0°).SNR=5dB，快拍數(shù)為500條件下的空間譜輸出如圖2(a)所示，可以看到本文所提算法和MUSIC算法均能分辨兩個(gè)信源信號(hào).然而當(dāng)把SNR降為-5dB時(shí)，如圖2(b)所示，只有本文所提算法仍能很好地區(qū)分兩個(gè)信源信號(hào)，而MUSIC算法則模糊了兩個(gè)信號(hào)的譜峰.這說明本文所提算法具有更高的分辨率特性.

實(shí)驗(yàn)2高斯白噪聲下角度估計(jì)性能比較

考慮兩個(gè)窄帶不相關(guān)信源DOA估計(jì)情況，入射角度分別為θ1=-10°，θ2=30°，相對(duì)應(yīng)的極化參數(shù)與實(shí)驗(yàn)1相同.DOA估計(jì)均方根誤差隨SNR和快拍數(shù)的變化關(guān)系分別如圖3(a)和圖3(b)所示.可以看出，本文所提算法的DOA估計(jì)性能雖然略低于MODE算法，但在整個(gè)SNR和采用快拍數(shù)域全面優(yōu)于ESPRIT算法.同時(shí)需要說明的是，MODE算法需要一個(gè)較好的初始估計(jì)才能保證全局最優(yōu)性，而本文所提算法由于基于凸優(yōu)化理論框架，無(wú)需較好的初始條件.

實(shí)驗(yàn)3高斯色噪聲下角度估計(jì)性能分析

除背景噪聲假定為高斯色噪聲外，其它實(shí)驗(yàn)條件與實(shí)驗(yàn)2相同.噪聲協(xié)方差矩陣N的第(i1,i2)個(gè)元素為

(53)

實(shí)驗(yàn)4兩種特殊入射情況下的估計(jì)結(jié)果

考慮兩種特殊的信源信號(hào)入射情況，即(1)入射角度不同，極化參數(shù)相同；(2)入射角度相同，極化角度不同.SNR=10dB，快拍數(shù)為1000下10次獨(dú)立的仿真結(jié)果分別如圖5(a)和圖5(b)所示.圖5(a)中三個(gè)信源的入射角度分別為θ1=-10°，θ2=30°，θ3=45°，極化參數(shù)(γ,η)均為(30°,0°)，可以看出信號(hào)入射角度不同時(shí)本文算法不僅可以很容易地區(qū)分三個(gè)信源信號(hào)，而且可以獲得很好的極化角度估計(jì).圖5(b)中三個(gè)信源的入射角度分別為θ1=-10°，θ2=30°，θ3=30°，相對(duì)應(yīng)的極化參數(shù)(γ,η)為(30°,0°)，(45°,0°)和(0°,0°).可以看出，雖然圖5(b)中只出現(xiàn)了兩個(gè)信源位置信息，但我們發(fā)現(xiàn)ρ1?ρ2=0.2480，這說明在極化角度不同時(shí)，即便兩個(gè)信源的入射角度一致，所提算法仍能有效地區(qū)分它們.換句話說，本文算法可以更好的利用極化信息，這是其它算法無(wú)法比擬的.

5　結(jié)論

本文在稀疏信號(hào)重構(gòu)理論框架下提出了一種新的DOA和極化角度估計(jì)算法.該算法在累積量域構(gòu)建稀疏表示模型，并分別基于加權(quán)l(xiāng)1范數(shù)最小化和Zhang懲罰策略實(shí)現(xiàn)高斯白/色噪聲下的DOA和極化角度估計(jì).借助估計(jì)的極化角度信息，該算法還可以有效區(qū)分兩個(gè)入射角度一樣的信源信號(hào).計(jì)算機(jī)仿真實(shí)驗(yàn)驗(yàn)證了所提算法的有效性，并展示了其在分辨率、色噪聲下估計(jì)精度等方面的優(yōu)勢(shì).

[1]Krim H,Viberg M.Two decades of array signal processing research:the parametric approach[J].IEEE Signal Processing Magazine,1996,13(4):67-94.

[2]Yaqoob M A,Mannesson A,Bemhardsson B,et al.On the performance of random antenna arrays for direction of arrival estimation[A].Proceedings of the 2014 IEEE International Conference on Communications Workshops (ICC)[C].USA:IEEE Press,2014.193-199.

[3]Qian C,Huang L,So H C.Improved unitary root-MUSIC for DOA estimation based on pseudo-noise resampling[J].IEEE Signal Processing Letters,2014,21(2):140-144.

[4]Xi N,Li L.A computationally efficient subspace algorithm for 2-D DOA estimation with L-shaped array[J].IEEE Signal Processing Letters,2014,21(8):971-974.

[5]Gou X,Liu Z,Xu Y.Biquaternion cumulant-MUSIC for DOA estimation of noncircular signals[J].Signal Processing,2013,93(4):874-881.

[6]Manikas A,Ng J W.Cross-dipole arrays for asynchronous DS-CDMA systems[J].IEEE Transactions on Antennas and Propagation,2004,52(1):122-131.

[7]Giuli D.Polarization diversity in radars[J].Proceedings of the IEEE,1986,74(2):245-269.

[8]Nehorai N,Paldi E.Vector-sensor array processing for electro-magnetic source localization[J].IEEE Transactions on Signal Processing,1994,42(2):376-398.

[9]Li J,Compton Jr R T.Angle and polarization estimation using ESPRIT with a polarization sensitive array[J].IEEE Transactions on Antennas and Propagation，1991,39(9):1376-1383.

[10]Lin C,Fang W,Yang W,et al.SPS-ESPRIT for joint DOA and polarization estimation with a COLD array[A].Proceedings of the IEEE International Symposium on Antennas and Propagation[C].USA:IEEE Press,2007.1136-1139.

[11]齊子森,郭英,姬偉峰,等.錐面共形陣列天線盲極化DOA估計(jì)算法[J].電子學(xué)報(bào),2009,37(9):1919-1925.

Qi Zi-sen,Guo Ying,Ji Wei-feng,et al.Blind DOA estimation algorithm for conical conformal array antenna with respect to polarization diversity[J].Acta Electronica Sinica,2009,37(9):1919-1925.(in Chinese)

[12]Miron S,Le B N,Mars J I.Quaternion-MUSIC for vector-sensor array processing[J].IEEE Transactions on Signal Processing,2006,54(4):1218-1229.

[13]Le B N,Maron S,Mars J I.MUSIC algorithm for vector-sensors array using biquaternions[J].IEEE Transactions on Signal Processing,2007,55(9):4523-4533.

[14]龔曉峰,等.電磁矢量傳感器陣列信號(hào)波達(dá)方向估計(jì)：雙模MUSIC[J].電子學(xué)報(bào),2008,36(9):1698-1703.

Gong Xiao-feng,et al.Direction-finding with electromagnetic vector-sensor array:dual-mode MUSIC[J].Acta Electronica Sinica,2008,36(9):1698-1703.(in Chinese)

[15]Li J,Stoica P,Zheng D.Efficient direction and polarization estimation with a COLD array[J].IEEE Transactions on Antennas and Propagation,1996,44(4):539-547.

[16]Malioutov D,et al.A sparse signal reconstruction perspective for source localization with sensor arrays[J].IEEE Transactions on Signal Processing,2005,53(8):3010-3022.

[17]Xu X,Wei X H,Ye Z F.DOA estimation based on sparse signal recovery utilizing weighted l1-norm penalty[J].IEEE Signal Processing Letters,2012,19(3):155-158.

[18]田野,孫曉穎,秦宇鏑.基于兩步加權(quán)L1范數(shù)約束的高分辨率波達(dá)方向和功率估計(jì)[J].電子與信息學(xué)報(bào),2014,36(7):1637-1641.

Tian Ye,Sun Xiao-ying,Qin Yu-di.High-resolution direction-of-arrival and power estimation using two-stage weighted L1penalty[J].Journal of Electronics & Information Technology,2014,36(7):1637-1641.(in Chinese)

[19]Hyder M M,Mahata K.Direction-of-arrival estimation using a mixed l2,0norm approximation[J].IEEE Transactions on Signal Processing,2010,58(9):4646-4655.

[20]Tian Y,et al.DOA and power estimation using a sparse representation of second-order statistics vector and l0-norm approximation[J].Signal Processing,2014,105:98-108.

[21]Gorodnitsky I F,Rao B D.Sparse signal reconstruction from limited data using FOCUSS:A re-weighted minimum norm algorithm[J].IEEE Transactions on Signal Processing,1997,45(3):600-616.

[22]Gasso G,Rakotomamonjy A,Canu S.Recovering sparse signals with a certain family of non-convex penalties and DC programming[J].IEEE Transactions on signal processing,2009,57(12):4686-4698.

[23]Wax M,Kailath T.Detection of signals by information theoretic criteria[J].IEEE Transactions on Acoustics,Speech and Signal Processing,1985,33(2):387-392.

[24]Nemirovski A,Ben T A.Lectures on Modern Convex Optimization:Analysis,Algorithms and Engineering Applications[M].Philadelphia,PA:SIAM,2001.

[25]Zhang T.Some sharp performance bounds for least squares regression with l1regularization[J].Annals of Statistics,2009,37(5A):2109-2144.

[26]Shao J.Linear model selection by cross-validation[J].Journal of the American Statistical Association,1993,88(422):486-494.

田野男,1985年12月生于河北省平泉縣,2014年于吉林大學(xué)通信與信息系統(tǒng)專業(yè)獲得博士學(xué)位，現(xiàn)為燕山大學(xué)信息科學(xué)與工程學(xué)院講師，主要研究方向?yàn)殛嚵行盘?hào)處理、稀疏信號(hào)重構(gòu)等.

E-mail:tianye@ysu.edu.cn

練秋生男,1969年8月生于江西遂川,工學(xué)博士，現(xiàn)為燕山大學(xué)信息科學(xué)與工程學(xué)院教授、博士生導(dǎo)師,主要研究方向?yàn)橄∈璞硎?、壓縮感知等.

E-mail:lianqs@ysu.edu.cn

徐鶴女,1989年12月生于遼寧省丹東市,2013年于吉林大學(xué)通信與信息系統(tǒng)專業(yè)獲碩士學(xué)位,主要研究方向?yàn)殛嚵行盘?hào)處理、空間譜估計(jì)等.

E-mail:xuhebest@sina.com

DOA and Polarization Angle Estimation AlgorithmBased on Sparse Signal Reconstruction

TIAN Ye,LIAN Qiu-sheng,XU He

(School of Information Science and Engineering,Yanshan University,Qinhuangdao,Hebei 066004,China)

Existing direction-of-arrival and polarization estimation methods mostly rely on subspace technique.This paper proposes a novel DOA and polarization angle estimation algorithm from sparse signal reconstruction perspective.The algorithm first constructs a cumulant matrix model which is only related to DOA parameter,and then obtains DOA estimation using the weighted l1-norm minimization.Further,this paper constructs another three cumulant vector models by sum-average arithmetic,and enforces sparsity by Zhang penalty,which leads to almost unbiased polarization angle estimation.Meanwhile,this paper also demonstrates how to identify two sources with same DOA using their polarization characteristics.Computer simulation results validate the effectiveness of the proposed algorithm.

DOA and polarization angle estimation；sparse signal reconstruction；weighted l1-norm；Zhang penalty

2015-01-28；

2015-04-08；責(zé)任編輯：孫瑤

國(guó)家自然科學(xué)基金(No.61471313)；河北省高等學(xué)校青年拔尖人才計(jì)劃(No.BJ2016051);燕山大學(xué)校內(nèi)自主研究課題(No.14LGA012)；燕山大學(xué)博士基金(No.B887)

TN911.7

A

0372-2112 (2016)07-1548-07

??學(xué)報(bào)URL:http://www.ejournal.org.cn

10.3969/j.issn.0372-2112.2016.07.004