代 林,崔 琛,余 劍,梁 浩
(電子工程學院通信對抗系,安徽合肥230037)
基于貝葉斯壓縮感知的CSR穩(wěn)健參數(shù)估計方法
代 林,崔 琛,余 劍,梁 浩
(電子工程學院通信對抗系,安徽合肥230037)
針對“完全擾動”情況下壓縮感知雷達(compressed sensing radar,CSR)觀測矢量和感知矩陣嚴重失配,進而引起參數(shù)估計性能急劇下降的問題,提出了一種基于貝葉斯壓縮感知(Bayesian compressed sensing,BCS)的穩(wěn)健參數(shù)估計方法。首先構造“完全擾動”情況下CSR參數(shù)估計的稀疏線性模型,并從稀疏矢量的最大后驗概率(maximum a posteriori,MAP)出發(fā),推導了完全擾動矩陣服從柯西分布時的優(yōu)化目標函數(shù);隨后通過稀疏矢量和尺度參數(shù)的交替迭代,求得稀疏矢量的最優(yōu)解。與現(xiàn)有重構算法及其改進算法相比,該方法能夠有效改善CSR系統(tǒng)應對失配誤差的穩(wěn)健性,提高目標成功檢測的概率和參數(shù)估計的精度。計算機仿真實驗驗證了該方法的有效性和魯棒性。
壓縮感知雷達;完全擾動;柯西分布;Lorentzian范數(shù);交替迭代
壓縮感知雷達(compressed sensing radar,CSR)是近年來提出的一種新型體制雷達[1],它將壓縮感知(compressed sensing,CS)[2]技術應用到雷達信號處理(radar signal processing,RSP)[34]中,在硬件要求、旁瓣抑制、目標檢測和參數(shù)估計等方面相比于傳統(tǒng)雷達具有顯著的優(yōu)勢。CSR在接收端以信息率進行測量感知,無需匹配濾波,降低了系統(tǒng)對硬件的性能要求,經(jīng)稀疏目標場景的重構進行參數(shù)估計,能夠同時獲得時延和多普勒頻率的超分辨估計。
文獻[5]指出,測量感知和重構算法的設計構成了CSR研究的主要內(nèi)容,其中測量感知與測量方式及發(fā)射波形有關,重構算法依賴于目標函數(shù)的構造和優(yōu)化求解。如果已知測量方式和發(fā)射波形,重構算法的選擇或者優(yōu)化設計就顯得尤為重要。然而,現(xiàn)有大多數(shù)文獻主要針對“加性擾動”(觀測噪聲)下CSR的參數(shù)估計問題進行研究,如迭代重加權最小二乘算法(iteratively re-weighted least squares,IRLS)[6]等,但在實際應用中還存在著“乘性擾動”,它將造成觀測矢量和感知矩陣嚴重失配,導致成功重構的概率顯著降低、重構的誤差明顯增大。文獻[7]構造了感知矩陣中存在噪聲的“完全擾動”稀疏線性模型,從理論上分析了該模型解的界的問題,并由基追蹤(basis pursuit,BP)驗證了“乘性擾動”對重構算法的影響。文獻[8]推導了上述模型下稀疏信號重構的克拉美羅界(Cramér-Rao bound,CRB)。但M.Herman和Y.Tang都只是針對擾動模型本身展開研究,并沒有給出適合求解的重構算法。文獻[9]將總體最小二乘(total least squares,TLS)方法應用到重構算法的設計中,提出了適合上述模型優(yōu)化求解的稀疏TLS(sparse TLS,S-TLS)方法,通過稀疏矢量和擾動矩陣的交替迭代求得稀疏矢量。但是稀疏信號的重構沒有限制偽峰的措施,擾動矩陣的求解是已知矢量求矩陣的欠定過程,難以保證求解結果的有效性。類似地,文獻[10]采用貝葉斯框架和TLS方法對欠定系統(tǒng)聚焦求解算法(focal underdetermined system solver,F(xiàn)OCUSS)進行擴展,提出了基于總體最小二乘的欠定系統(tǒng)聚焦求解算法(TLS-FOCUSS),在保持FOCUSS收斂性的同時能夠處理稀疏線性模型“完全擾動”的情況。文獻[11]也針對測量矩陣中存在噪聲的情況進行了研究。但是,上述文獻都是在完全擾動矩陣服從高斯正態(tài)分布的假設下開展研究的,而在實際應用中感知矩陣中包含的噪聲更加復雜,更適合采用對稱α穩(wěn)定(symmetric alpha stable,SαS)分布進行描述[12]。這方面的研究見諸報導的還十分有限,主要針對的也是觀測噪聲服從SαS分布的情況。文獻[13]基于Hampel函數(shù)和l1范數(shù)稀疏正則化,提出了基于M估計的最小化絕對收縮和選擇算子(M-estimation-based LASSO,M-LASSO),也稱基于M估計的IRLS(M-estimation-based IRLS,Me-IRLS)。當脈沖能量較小時,該算法重構精度較高,但隨著脈沖能量增大,估計值的稀疏性急劇惡化。文獻[14]選擇由Lorentzian范數(shù)代替Hampel函數(shù)對觀測矢量中的異值元素進行回降壓縮,提出了基于Lorentzian范數(shù)的BP(Lorentzian-based BP,LBP)算法[14]和基于Lorentzian范數(shù)的迭代硬閾值(Lorentzian-based iterative hard thresholding,LIHT)[15]算法等兩種方法。其中LBP未對重構值中的偽峰進行限制,在較低信噪比條件下的重構誤差較大,且計算復雜度較高,無法滿足求解效率的要求;LIHT雖然求解效率滿足要求,但與IHT一樣,需要預知待求變量的稀疏度,且不能從隨機投影的觀測序列中恢復出待求變量。
針對“完全擾動”情況下CSR參數(shù)估計性能急劇下降的問題,本文提出了一種基于貝葉斯壓縮感知(Bayesian compressed sensing,BCS)的穩(wěn)健參數(shù)估計方法,與已有算法相比,能夠有效改善CSR應對失配誤差的穩(wěn)健性,提高目標成功檢測的概率和參數(shù)估計的精度。
假設遠場照射區(qū)域內(nèi)包含Nτ個距離分辨單元和ND個多普勒分辨單元,其中分辨單元(nτ,nD)的時延、多普勒頻率和復散射系數(shù)分別為,和,不考慮信道衰減,各分辨單元的回波信號在接收端疊加構成CSR的回波信號,有
式中,s(t)為發(fā)射信號,s=[s(0),s(1),…,s(Ns-1)]T為離散發(fā)射信號矢量,Ns為采樣長度。
由式(4)可知,ΔΦx是擾動矩陣與回波信號的乘積,將對θ的重構產(chǎn)生“乘性擾動”。令B=[ΔΦ e]表征CSR系統(tǒng)面臨的“加性擾動”和“乘性擾動”,稱之為完全擾動矩陣,本文假設B中的每個元素都是零均值獨立同分布的SαS隨機變量,其特征函數(shù)為
式中,α∈(0,2]為特征指數(shù),決定SαS分布的沖擊程度;μ∈(-∞,+∞)為位置參數(shù),當0<α≤1時表征SαS分布的中值,當1<α≤2時表征SαS分布的均值;γ∈(0,+∞)為分散系數(shù),表征穩(wěn)定分布中各變量偏離中值或均值的程度。一般情況下,SαS隨機變量沒有封閉的概率密度函數(shù),無法建立回波信號的貝葉斯模型,難以獲得θ的最大后驗概率估計。
根據(jù)貝葉斯理論,稀疏信號的重構要求模型具有封閉的概率密度函數(shù),而SαS分布只有在α的值為2、1和0.5時才能滿足要求,分別對應高斯正態(tài)分布、柯西分布和Levy分布。在高斯分布的假設下,觀測矢量中的異值元素將被線性放大,由此獲得的估計值可能出現(xiàn)無窮大的偏差;Levy分布的概率密度函數(shù)中,指數(shù)項和被除項都包含未知量,目標函數(shù)難以優(yōu)化求解;而柯西分布可以大致反映SαS分布在1≤α≤2范圍內(nèi)統(tǒng)計特征,由此建立的目標函數(shù)能夠有效地避免上述兩個問題。因此,下面將以柯西分布為假設前提,構建式(4)中稀疏感知模型的穩(wěn)健優(yōu)化目標函數(shù)并進行求解。
2.1 目標參數(shù)的貝葉斯模型
文獻[16]指出,可以認為稀疏矢量θ是獨立同分布的廣義高斯隨機矢量,其最大后驗概率估計可以表示為
式中,C為常數(shù)。將式(7)帶入式(6),可得
2.2 模型求解
然而,式(8)中的尺度參數(shù)包含待求變量,不能直接求解。為此,本文首先利用變量分裂技術[17]進行變換,得到等價的約束優(yōu)化問題
然后,通過交替迭代進行求解,即先固定尺度參數(shù)尋找稀疏矢量的估計值,再固定稀疏矢量計算尺度參數(shù)的估計值,循環(huán)往復,直至收斂。
式中,Wk+1=W()和Qk+1=Q()為對角陣,且
由定點法原理,可以得到式(12)的迭代形式
將式(16)中uk+1帶入式(11)以更新尺度參數(shù)ζk+1,重復式(11)和式(16)直到滿足終止條件,得到的解即為式(8)中欠定方程的最優(yōu)解。
3.1 收斂性證明
定理在迭代過程中,優(yōu)化目標函數(shù)L(θ)是隨著迭代次數(shù)單調(diào)遞減的,即L(θk+1)<L(θk),其中k為迭代次數(shù)。
證明由前面的描述,將L(θ)改寫為稀疏矢量u和尺度參數(shù)ζ交替迭代的形式,即
由于L(u,ζ)中包含兩個未知變量,必須保證L(u,ζ)關于u和ζ都是減小的,才能證明L(u,ζ)是隨著迭代次數(shù)單調(diào)遞減的,即L(θ)單調(diào)遞減。
經(jīng)過l+1和l次內(nèi)循環(huán)迭代后,可得L1(uk,l+1)和L1(uk,l),令二者相減,由文獻[19]中的推論2和文獻[20]中的引理1可得
證畢
3.2 算法的計算復雜度
為便于描述,在本小節(jié)中記M=Np和N=NτND。本文算法的計算量主要在于交替迭代過程中稀疏矢量uk+1的計算,包括6次矩陣相乘運算,一次矩陣求逆運算,一次矩陣與向量相乘運算和一次標量與矩陣相乘運算,單次迭代需要O(M3+2NM2+M2+4NM+M)次乘法。此外,式(16)中的閾值收縮有助于提高算法的收斂速率,將進一步降低計算復雜度。S-TLS算法的求解過程包含擾動矩陣的估計,其計算復雜度為O(N3+N2+NM2+4NM);TLSFOCUSS求解過程包含高維矩陣的特征值分解,計算復雜度為O((N+1)3+(M+1)(N+1)2+2M(N+1)+2M2)。
首先作如下定義:①若估計值中前K個峰值對應著遠場照射區(qū)域內(nèi)目標的散射系數(shù),且第K和K+1個峰值之比大于2.5,則認為目標成功檢測,反之不能檢測;②估計值的均方根誤差(root mean square error,RMSE)定義為
在仿真中,假設CSR載波頻率f0=10 GHz,工作波長λ=0.03 m,發(fā)射信號s為立方碼信號,碼元長度為Ns=50;遠場照射區(qū)內(nèi)的距離、多普勒分辨單元分別為Nτ=25和ND=31,利用(|x/ND|,x%ND)可以隨機產(chǎn)生分辨單元不同的目標,其散射系數(shù)由復高斯分布CN(0,1)產(chǎn)生,其中x?[1,NτND]為互不相鄰的隨機數(shù)列,數(shù)列中元素的個數(shù)代表目標的個數(shù);測量矩陣Φ中的元素都是獨立同分布的零均值高斯隨機變量且Np=50。
實驗1本文方法的有效性驗證。
假設CSR探測區(qū)域存在5個目標,位置坐標和散射系數(shù)由上文方法產(chǎn)生。圖1和圖2為完全擾動矩陣服從高斯分布(擾動方差σ2=0.062 5)和柯西分布(擾動離差γ=0.025)時,IRLS、LBP、S-TLS、TLS-FOCUSS和本文方法參數(shù)估計的結果,圖中橫坐標“索引值”表示不同的時延-多普勒分辨單元。
由圖1和圖2可以看出,CSR的參數(shù)估計在“完全擾動”情況下受到嚴重影響。當完全擾動矩陣服從高斯隨機分布時,IRLS和LBP的估計值中出現(xiàn)很多偽峰,只能準確地區(qū)分出少數(shù)強散射目標,如圖1(a)和圖1(b);S-TLS、TLS-FOCUSS以及本文方法都可以比較準確地估計出目標的參數(shù),如圖1(c)、圖1(d)和圖1(e)。這是因為,IRLS和LBP在建立優(yōu)化目標函數(shù)時就沒有考慮“完全擾動”的情況,不能克服觀測矢量的“乘性擾動”;S-TLS、TLS-FOCUSS以及本文方法都是將觀測矢量中的“乘性擾動”也作為優(yōu)化參量,通過各優(yōu)化參量的交替迭代,降低“乘性擾動”對目標參數(shù)估計的影響。
圖1 高斯擾動(σ2=0.25)下,CSR的參數(shù)估計結果
當完全擾動矩陣服從柯西隨機分布時,IRLS和LBP參數(shù)估計的性能進一步惡化,如圖2(a)和圖2(b);由于S-TLS和TLS-FOCUSS是在高斯擾動的假設下將TLS思想應用到Lasso和FOCUSS中,對于柯西擾動不具備魯棒性,目標大多湮沒在偽峰中,難以有效區(qū)分,如圖2(c)和圖2(d);而本文方法則是在柯西分布的假設下,由稀疏矢量的最大后驗概率建立優(yōu)化目標函數(shù),形成類似Lorentzian范數(shù)的誤差擬合函數(shù),可以有效地抑制觀測矢量中的模型噪聲和異值元素,如圖2(e)。
圖2 柯西擾動(γ=0.025)下,CSR的參數(shù)估計結果
實驗2不同的完全擾動矩陣下,各算法參數(shù)估計性能的比較。
假設CSR探測區(qū)域存在5個隨機分布的點目標,其他參數(shù)設置與實驗1相同。圖3為完全擾動矩陣服從高斯正態(tài)分布時,各算法目標成功檢測的概率和估計值的RMSE隨著擾動方差變化的曲線,不同擾動方差條件下蒙特卡羅實驗次數(shù)為500。從圖中可以看出,在仿真實驗限定的范圍內(nèi)IRLS和LBP目標成功檢測的概率和估計值的RMSE隨著擾動方差的增大而急劇下降;S-TLS、TLS-FOCUSS以及本文方法都具有較好的參數(shù)估計性能,其中本文方法目標成功檢測的概率介于S-TLS和TLS-FOCUSS之間,估計結果的RMSE略高。這是因為,由于IRLS和LBP未考慮“乘性擾動”,目標的參數(shù)估計性能將隨擾動方差急劇惡化;與本文方法中的Lorentzian范數(shù)相比,S-TLS和TLS-FOCUSS中的l2范數(shù)更適合去擬合高斯擾動帶來的誤差,可以有效克服模型噪聲對參數(shù)估計的影響,此外,S-TLS由l1范數(shù)來衡量稀疏度,僅適用于復散射系數(shù)矢量服從拉普拉斯先驗分布的情況,而TLS-FOCUSS由廣義高斯分布推導構造優(yōu)化目標函數(shù),利用lp范數(shù)可以更準確地刻畫復散射系數(shù)矢量的稀疏性。
圖3 CSR參數(shù)估計性能隨著擾動方差的變化曲線
圖4為完全擾動矩陣服從柯西隨機分布時,各算法目標成功檢測的概率和估計值的RMSE隨著擾動方差變化的曲線,不同擾動方差條件下蒙特卡羅實驗次數(shù)為500。如果目標不能成功檢測,估計值的RMSE視為無窮大。從圖中可以看出,LBP可以有效克服離差很小的柯西擾動,當擾動離差γ≥0.02時,IRLS、LBP、S-TLS以及TLS-FOCUSS都難以有效地進行目標的檢測,而本文方法在γ≥0.04時仍能保持較高的成功檢測概率和估計精度。這是因為,Lorentzian范數(shù)是軟回降的誤差擬合函數(shù),可以有效克服異值元素的影響,由此LBP可以有效克服較小的擾動離差引起的“乘性擾動”,具有一定的魯棒性;IRLS、S-TLS以及TLS-FOCUSS則在原理上就不能克服觀測矢量中的乘性擾動,只有在更小的擾動離差下才具有一定的成功檢測概率;本文方法由目標參數(shù)的最大后驗概率估計建立基于Lorentzian范數(shù)的穩(wěn)健目標函數(shù),通過尺度參數(shù)和目標參數(shù)的交替尋優(yōu),有效去除柯西擾動對參數(shù)估計的影響。
圖4 CSR參數(shù)估計性能隨著擾動離差的變化曲線
圖5為完全擾動矩陣服從SαS隨機分布(擾動離差γ=0.025)時,各算法目標成功檢測的概率和估計值的RMSE隨著特征指數(shù)變化的曲線,不同擾動方差條件下蒙特卡羅實驗次數(shù)為500。從圖中可以看出,本文方法在仿真實驗限定的范圍內(nèi)(1.0≤α≤2.0)都可以得到較好的成功檢測概率和估計值RMSE,而IRLS、LBP、S-TLS和TLS-FOCUSS等算法要求更大的α值;隨著α值不斷增大,各算法參數(shù)估計的性能越來越好,在α=2.0時可以得到最優(yōu)的成功檢測概率和估計值RMSE。
圖5 CSR參數(shù)估計性能隨著特征指數(shù)的變化曲線
這是因為,1.0≤α≤2.0范圍內(nèi)的SαS隨機分布可由柯西分布近似表征,而本文方法是在柯西擾動的假設下推導而來,理應可以有效克服SαS隨機擾動的影響;隨著α值逐漸增大,觀測矢量中異值元素不斷減小,IRLS、S-TLS以及TLS-FOCUSS下被線性放大的元素不斷減小,乘性擾動造成的影響越來越小,目標參數(shù)的估計越來越準確;當α=2.0時,完全擾動矩陣服從高斯正態(tài)分布,其方差σ=γ/2,幾乎不能對觀測矢量產(chǎn)生擾動。
針對“完全擾動”情況下CSR參數(shù)估計性能急劇下降的問題,本文提出了一種基于貝葉斯壓縮感知的穩(wěn)健參數(shù)估計方法,能夠有效改善CSR系統(tǒng)應對失配誤差的穩(wěn)健性,提高目標成功檢測的概率和參數(shù)估計的精度。與IRLS和LBP等相比,能夠有效降低觀測矢量和感知矩陣失配對參數(shù)估計性能的影響;與S-TLS和TLS-FOCUSS等相比,適用范圍可以擴展到由SαS分布描述的完全擾動矩陣。
[1]Baraniuk R,Steeghs P.Compressive radar imaging[C]∥Proc. of the IEEE Radar Conference,2007:128- 133.
[2]Donoho D L.Compressed sensing[J].IEEE Trans.on Information Theory,2006,52(4):1289- 1306.
[3]Ender J H G.On compressive sensing applied to radar[J].Signal Processing,2010,90(5):1402- 1414.
[4]Strohmer T,F(xiàn)riedlander B.Analysis of sparse MIMO radar[J].Applied and Computational Harmonic Analysis,2014,37(3):361- 388.
[5]Zhang J D,Zhang G,Pan H,et al.Optimized sensing matrix design of filter structure based compressed sensing radar[J].Acta Aeronautica et Astronautica Sinica,2013,34(4):864- 872.(張勁東,張弓,潘匯,等.基于濾波器結構的壓縮感知雷達感知矩陣優(yōu)化[J].航空學報,2013,34(4):864- 872.)
[6]Wipf D,Nagarajan S.Iterative reweighted l1and l2methods for finding sparse solutions[J].IEEE Journal of Selected Topics in Signal Processing,2010,4(2):317- 329.
[7]Matthew H,Thomas S.General deviants:an analysis of perturbations in compressed sensing[C]∥Proc.of the IEEE Journal of Selected Topics in Signal Processing,Special Issue on Compressed Sensing,2010,4(2):342- 349.
[8]Tang Y J,Chen L M,Gu Y T.On the performance bound of sparse estimation with sensing matrix perturbation[J].IEEE Trans.on Signal Processing,2013,61(17):4372- 4386.
[9]Zhu H,Geert L,Georgios G.Sparsity-cognizant total leastsquares for perturbed compressive sampling[J].IEEE Trans.on Signal Processing,2011,59(5):2002- 2016.
[10]Han X,Zhang H,Meng H.TLS-FOCUSS for sparse recovery with perturbed dictionary[C]∥Proc.of the International Conference on Acoustics Speech and Signal Processing,2011:3952- 3955.
[11]Jiang L B,Huang T.New compressive sampling matching pursuit algorithm[J].Application Research of Computers,2013,30(2):402- 404.(蔣留兵,黃韜.一種新的壓縮采樣匹配追蹤算法[J].計算機應用研究,2013,30(2):402- 404.)
[12]Shao M,Nikias C L.Signal processing with fractional lower order moments:stable processes and their applications[J].Pro-ceedings of the IEEE,1993,81(7):986- 1010.
[13]Zhang Z G,Chan S C,Zhou Y,et al.Robust linear estimation using M-estimation and weighted L1regularization:model selection and recursive implementation[C]∥Proc.of the IEEE International Symposium on Circuits and Systems,2009:1193- 1196.
[14]Carrillo R E,Barner K E,Aysal T C.Robust sampling and reconstruction methods for sparse signals in the presence of impulsive noise[J].IEEE Trans.on Selected Topics in Signal Processing,2010,4(2):392- 408.
[15]Carrillo R E,Barner K E.Lorentzian iterative hard thresholding:robust compressed sensing with prior information[J].IEEE Trans.on Signal Processing,2013,61(17/20):4822- 4833.
[16]Rao B D,Engan K,Cotter S F,et al.Subset selection in noise based on diversity measure minimization[J].IEEE Trans.on Signal Processing,2003,51(3):760- 770.
[17]Li Y F,F(xiàn)eng X C.The split Bregman method for L1projection problems[J].Acta Electronica Sinica,2010,38(11):2471-2475.(李亞峰,馮象初.L1投影問題的分裂Bregman方法[J].電子學報,2010,38(11):2471- 2475.)
[18]Petersen K B,Pedersen M S.The matrix cookbook[EB/OL].http:∥matrixcookbook.com.
[19]Carrillo R E.Robust methods for sensing and reconstructing sparse signals[D].USA:University of Delaware,2011.
[20]Rao B,Engan K,Cotter S,et al.Subset selection in noise based on diversity measure minimization[J].IEEE Trans.on Signal Processing,2003,51(3):760- 770.
代 林(198-6- ),男,博士研究生,主要研究方向為壓縮感知雷達、自適應信號處理。
E-mail:dailin513@126.com
崔 ?。?962- ),男,教授,碩士,主要研究方向為雷達信號處理、通信信號處理以及通信系統(tǒng)仿真。
E-mail:dailinincui@126.com
余 劍(1980- ),男,講師,碩士,主要研究方向為雷達信號處理、射頻信號處理。
E-mail:dailininyu@126.com
梁 浩(198-7- ),男,博士研究生,主要研究方向為陣列信號處理、MIMO雷達信號處理。
E-mail:lhmailhappy@163.com
Robust parameter estimation method for CSR based on Bayesian compressed sensing
DAI Lin,CUI Chen,YU Jian,LIANG Hao
(Department of Communication Countermeasure,Institute of Elctronic Engineering,Hefei 230037,China)
In practical application,the mismatch between measurement vector and sensing matrix caused by completely perturbed observations will result in a sharp decline in the performance of parameter estimation for compressed sensing radar(CSR).A novel robust parameter estimation algorithm is proposed based on Bayesian compressed sensing(BCS).The completely perturbed sparse linear model is firstly built,and the robust target function is derived with the maximum a posteriori(MAP)of the sparse vector when the completely perturbed matrix obeys Cauchy distribution.Then the optimal solution is achieved through the alternate iteration between the sparse vector and the scale parameter.Compared with most existing recovery algorithm and their derivants,the proposed method effectively improves the robustness against the foregoing mismatch,increases the target detection probability and reduces the estimation error.The effectiveness of the proposed algorithms is demonstrated by computer simulations.
compressed sensing radar(CSR);completely perturbation;Cauchy distribution;Lorentzian norm;alternate iteration
TN 958.8
A
10.3969/j.issn.1001-506X.2015.11.09
1001-506X(2015)11-2480-07
2014- 11- 24;
2015- 03- 01;網(wǎng)絡優(yōu)先出版日期:2015- 04- 29。
網(wǎng)絡優(yōu)先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20150429.1106.002.html
國家自然科學基金(60702015)資助課題