王 釗顧 紅蘇衛(wèi)民陳志仁陳金立

①（南京理工大學(xué)電子工程與光電技術(shù)學(xué)院 南京 210094）

②（南京信息工程大學(xué)電子與信息工程學(xué)院 南京 210044）

基于共軛噪聲組的寬帶噪聲雷達(dá)機(jī)動(dòng)目標(biāo)參數(shù)估計(jì)

王 釗*①顧 紅①蘇衛(wèi)民①陳志仁①陳金立②

①（南京理工大學(xué)電子工程與光電技術(shù)學(xué)院 南京 210094）

②（南京信息工程大學(xué)電子與信息工程學(xué)院 南京 210044）

寬帶噪聲雷達(dá)參數(shù)估計(jì)時(shí)通常會(huì)采用寬帶互模糊函數(shù)的方法，但是在處理機(jī)動(dòng)目標(biāo)時(shí)這種方法需要在距離、速度及加速度3維搜索而導(dǎo)致運(yùn)算量巨大，為此該文提出一種基于共軛噪聲組的機(jī)動(dòng)目標(biāo)參數(shù)估計(jì)算法。該算法首先根據(jù)回波伸縮效應(yīng)預(yù)設(shè)多路通道，每路通道截取固定長(zhǎng)度的噪聲組內(nèi)信號(hào)進(jìn)行混頻，然后利用分?jǐn)?shù)階傅里葉變換（FrFT）估計(jì)混頻信號(hào)的多普勒相位，根據(jù)相位信息構(gòu)造補(bǔ)償函數(shù)，并對(duì)補(bǔ)償后的噪聲組信號(hào)利用頻域尺度相關(guān)（FSC）算法估計(jì)回波的時(shí)延，最后聯(lián)立多普勒相位及時(shí)延信息獲取目標(biāo)的距離、速度和加速度。該算法避免了目標(biāo)參數(shù)3維搜索的過程，無需時(shí)域重構(gòu)回波信號(hào)，較寬帶互模糊函數(shù)方法極大地降低了運(yùn)算量，整個(gè)算法都可通過快速傅里葉變換（FFT）實(shí)現(xiàn)，便于系統(tǒng)實(shí)時(shí)處理。仿真結(jié)果驗(yàn)證了該算法的有效性及優(yōu)勢(shì)。

寬帶噪聲雷達(dá)； 共軛噪聲組； 分?jǐn)?shù)階傅里葉變換； 頻域尺度相關(guān)

寬帶噪聲雷達(dá)由于具有測(cè)距精度高、距離速度無模糊、截獲概率低、抗干擾能力強(qiáng)等優(yōu)點(diǎn)，受到廣泛的關(guān)注［1，2］。隨著隱身戰(zhàn)斗機(jī)和超高速戰(zhàn)斗機(jī)的發(fā)展，對(duì)機(jī)動(dòng)目標(biāo)的有效探測(cè)是當(dāng)今雷達(dá)的研究熱點(diǎn)。為了提高探測(cè)性能，通常會(huì)采用長(zhǎng)時(shí)間積累的方法。對(duì)于寬帶噪聲雷達(dá)而言，會(huì)存在以下問題：（1）與匹配濾波器失配相關(guān)的多普勒敏感問題；（2）由回波包絡(luò)的伸縮變換而導(dǎo)致的多普勒色散現(xiàn)象；（3）長(zhǎng)時(shí)間積累過程中由于加速度引起的多普勒展寬；（4）由于多普勒域的采樣頻率較低而產(chǎn)生多普勒模糊。

近年來，很多學(xué)者在改善噪聲雷達(dá)動(dòng)目標(biāo)探測(cè)性能方面提出了解決方法。針對(duì)噪聲雷達(dá)的多普勒敏感問題，文獻(xiàn)［3］設(shè)計(jì)了一種寬帶相關(guān)噪聲雷達(dá)系統(tǒng)，通過可變的延遲線來控制參考信號(hào)的時(shí)延，并進(jìn)一步估計(jì)出回波的多普勒頻率。文獻(xiàn)［4］將多重信號(hào)分類（MU ltiple SIgnal Classification， MUSIC）算法運(yùn)用到該系統(tǒng)中來，從而提高了多普勒頻率的分辨力。文獻(xiàn)［5］將互模糊函數(shù)用到噪聲雷達(dá)的時(shí)延和多普勒聯(lián)合估計(jì)中來，從而有效地解決噪聲雷達(dá)的多普勒敏感問題。考慮到多普勒色散現(xiàn)象，文獻(xiàn)［6］利用分段短時(shí)相關(guān)的方法，將包絡(luò)的伸縮變換理解為時(shí)變的距離走動(dòng)，但是這種方法僅適用于多普勒色散積較小的情況。文獻(xiàn)［7］采用寬帶互模糊函數(shù)來聯(lián)合估計(jì)時(shí)延及速度，但是這種方法在用于噪聲信號(hào)時(shí)通常需要構(gòu)造不同尺度和時(shí)延的參考信號(hào)，運(yùn)算量較大，不利于工程實(shí)現(xiàn)。文獻(xiàn)［8］將寬帶互模糊函數(shù)擴(kuò)展到寬帶距離-速度-加速度互模糊函數(shù)來處理機(jī)動(dòng)目標(biāo)，但是這種方法需要3維參數(shù)搜索，運(yùn)算量巨大。

與之前基于互模糊函數(shù)算法的研究不同，本文給出一種基于共軛噪聲組的雷達(dá)信號(hào)處理算法，用于實(shí)現(xiàn)機(jī)動(dòng)目標(biāo)距離、速度及加速度的聯(lián)合估計(jì)。本算法首先根據(jù)伸縮效應(yīng)截取回波中的有效數(shù)據(jù)，然后將截取后的噪聲組內(nèi)信號(hào)混頻，并對(duì)混頻信號(hào)作分?jǐn)?shù)階傅里葉變換（Fractional Fourier Transform，F(xiàn)rFT）以估計(jì)回波多普勒信息，接著構(gòu)造相位補(bǔ)償函數(shù)并利用頻域尺度相關(guān)（Frequency-domain Scale Correlation， FSC）估計(jì)整段噪聲組信號(hào)的時(shí)延信息，最后通過已得的多普勒信息與時(shí)延信息聯(lián)合估計(jì)出目標(biāo)的距離、速度及加速度。本文算法采用先多普勒估計(jì)后時(shí)延估計(jì)的方式避免了互模糊函數(shù)的3維搜索過程，通過回波預(yù)處理及FSC方法補(bǔ)償回波包絡(luò)的伸縮效應(yīng)來替代互模糊函數(shù)的時(shí)域插值運(yùn)算，整個(gè)算法都可通過快速傅里葉變換（Fast Fourier Transform， FFT）實(shí)現(xiàn)，運(yùn)算復(fù)雜度較低。計(jì)算機(jī)仿真驗(yàn)證了本文算法的有效性。

2 寬帶機(jī)動(dòng)目標(biāo)建模

假設(shè)存在勻加速運(yùn)動(dòng)的單個(gè)點(diǎn)目標(biāo)，信號(hào)發(fā)射時(shí)目標(biāo)初始距離為R0，徑向速度為 v?c，加速度為a，那么雷達(dá)接收的基帶回波信號(hào)可以表示為［9］

其中U表示信號(hào)能量，R（τ）為u（ t）的歸一化自相關(guān)函數(shù)，ε=1-α。將式（4）代入式（3）可得

針對(duì)這種情況，文獻(xiàn)［8］給出了一種寬帶距離-速度-加速度互模糊函數(shù)的算法，該算法通過構(gòu)造與目標(biāo)速度及加速度對(duì)應(yīng)的參考信號(hào)，來達(dá)到匹配濾波的效果。由于通常情況下，目標(biāo)的距離，速度和加速度都是未知的，該算法需要進(jìn)行3維的參數(shù)空間搜索，并且在構(gòu)造參考信號(hào)時(shí)重采樣的運(yùn)算復(fù)雜度較高。為了解決這一問題，根據(jù)噪聲信號(hào)的特性，本文提出一種共軛噪聲交替發(fā)射的信號(hào)形式及對(duì)應(yīng)的參數(shù)估計(jì)算法。

3 共軛噪聲組機(jī)動(dòng)目標(biāo)參數(shù)估計(jì)

3.1 共軛噪聲組回波預(yù)處理

假設(shè)雷達(dá)發(fā)射信號(hào)S（ t）為連續(xù)的兩段寬帶共軛噪聲信號(hào)，其表達(dá)式為

其中，u（ t）=exp［jθ（ t ）］， θ（t）為隨機(jī)相位，f0為載頻，T表示每段信號(hào)的時(shí)寬，信號(hào)形式如圖1所示。

該共軛噪聲信號(hào)的基帶回波信號(hào)Sr（ t）可以表示為

其中n=1，2，…， L ， m=1，2，…， L， L為目標(biāo)數(shù)，φ（t）表示回波信號(hào)的多普勒信息。從回波數(shù)據(jù)中截取相鄰的兩段時(shí)寬為T的信號(hào)作混頻處理，為了方便論述，將混頻結(jié)果表示為回波包絡(luò)與多普勒信號(hào)乘積的形式

圖1 共軛噪聲基帶信號(hào)波形示意圖

當(dāng)T= T/αn時(shí)，兩段信號(hào)的混頻結(jié)果可以表示為

3.2 多普勒頻率估計(jì)

在回波預(yù)處理之后，當(dāng)混頻信號(hào)中的包絡(luò)部分近似為1時(shí)，剩下的只有多普勒頻率的部分Ffb（ t， n， m），根據(jù)式（1）回波的相位信息，φ（t）應(yīng)為

由式（13）可知，此時(shí)的多普勒頻率部分Ffb（ t， n， m）實(shí)際上為線性調(diào)頻信號(hào)，其一次相位項(xiàng)P1和二次相位項(xiàng)P2包含了目標(biāo)的時(shí)延，速度及加速度信息，通過對(duì)P1和P2的估計(jì)，進(jìn)而可得目標(biāo)的運(yùn)動(dòng)參數(shù)。分?jǐn)?shù)階傅里葉變換是近年來使用較多的線性調(diào)頻信號(hào)參數(shù)的時(shí)頻分析方法［11，12］，對(duì)式（13）給出的多普勒頻率部分Ffb（ t， n， m）進(jìn)行分?jǐn)?shù)階傅里葉變換，可得

3.3 頻域尺度相關(guān)運(yùn)算

噪聲信號(hào)的多普勒敏感特性可以通過多普勒補(bǔ)償來抑制，忽略常數(shù)相位項(xiàng)，由式（13）可得需要構(gòu)造的相位補(bǔ)償函數(shù)為

將式（16）代入，相位補(bǔ)償函數(shù)可以用P1和P2表示為

為了獲得最大的信噪比增益，這里拼接Sr（ t）的兩段信號(hào)到一起進(jìn)行聯(lián)合處理，并將構(gòu)造的相位補(bǔ)償函數(shù)與Sr（ t）相乘，可得

將相位補(bǔ)償以后的信號(hào)Src（ t）與ubi（ t）進(jìn)行尺度為α的相關(guān)運(yùn)算，可得

式中F［·］表示傅里葉變換，F(xiàn)-1［·］表示逆傅里葉變換，Src（ f）為Src（ t）的頻譜，R［·，α］表示以α為尺度的頻譜重構(gòu)運(yùn)算，具體重構(gòu)運(yùn)算的方法見文獻(xiàn)［15］，這里不再贅述。這種在頻域?qū)崿F(xiàn)信號(hào)重構(gòu)及相關(guān)運(yùn)算的過程本文稱之為頻域尺度相關(guān)（FSC）運(yùn)算。

3.4 時(shí)延估計(jì)誤差分析

在寬帶情況下，匹配濾波器的時(shí)寬仍然為2T，但回波信號(hào)由于伸縮效應(yīng)時(shí)寬變?yōu)?T/α，致使寬帶相關(guān)函數(shù)的時(shí)寬為2T+ 2 T/α，此時(shí)目標(biāo)實(shí)際時(shí)延應(yīng)為

其中τM表示匹配濾波器的時(shí)延測(cè)量值。式（23）中的α為與目標(biāo)速度相關(guān)的尺度因子，可見對(duì)于寬帶噪聲信號(hào)來說，時(shí)延的測(cè)量值與速度是耦合的，將α?n代入可得FSC運(yùn)算后時(shí)延的真實(shí)值與測(cè)量值之間的關(guān)系為

3.5 多普勒解模糊方法

當(dāng)FrFT的分?jǐn)?shù)階域欠采樣時(shí)，式（16）用目標(biāo)參數(shù)來表示為

其中η表示分?jǐn)?shù)階域模糊數(shù)。將式（26）代入式（18）可得修正的相位補(bǔ)償函數(shù)為

當(dāng)相位補(bǔ)償函數(shù)不準(zhǔn)確時(shí)，Src（ t）與參考信號(hào)是失配的，因此可以構(gòu)造代價(jià)函數(shù)來搜索η值

由于每路通道的目標(biāo)速度范圍是已知的，代價(jià)函數(shù)的搜索次數(shù)可以有效地降低。

在確定η值以后，將式（28）帶入式（26），并與式（25）聯(lián)立就可以得到目標(biāo)參數(shù)的估計(jì)值為

式（29）中的參數(shù)求解過程是按照an→ vn→τn的順序進(jìn)行的，這里認(rèn)為在計(jì)算下一個(gè)參數(shù)時(shí)上一個(gè)參數(shù)為已知量，可以直接代入。

3.6 虛假目標(biāo)

在之前的討論中，都是假設(shè)混頻信號(hào)為同一目標(biāo)回波的情況。當(dāng)混頻信號(hào)為不同目標(biāo)回波時(shí)，根據(jù)3.1節(jié)的討論結(jié)果，若兩目標(biāo)速度剛好與相同的多普勒通道匹配，且兩目標(biāo)時(shí)延滿足一定關(guān)系時(shí)，對(duì)其混頻信號(hào)進(jìn)行FrFT亦可以估計(jì)出相應(yīng)多普勒相位，但是通過式（27）的方法構(gòu)造的相位補(bǔ)償函數(shù)不能與任何目標(biāo)回波匹配，在后續(xù)的FSC運(yùn)算中將無法檢測(cè)出目標(biāo)，可以在FSC以后設(shè)置門限將此虛假目標(biāo)剔除。

3.7 算法步驟

綜上所述，圖2給出共軛噪聲組信號(hào)估計(jì)機(jī)動(dòng)目標(biāo)參數(shù)的流程圖，具體步驟為：

步驟1 根據(jù)不同的包絡(luò)伸縮情況分多個(gè)通道對(duì)回波信號(hào)進(jìn)行預(yù)處理，每個(gè)通道按照不同的速度選取時(shí)寬為的回波信號(hào)，根據(jù)目標(biāo)最大速度 vmax可得通道數(shù)為。

步驟5 根據(jù)P1和P2以及時(shí)延τMn的估計(jì)結(jié)果，利用式（29）逐次計(jì)算出目標(biāo)的距離、速度及加速度。

4 運(yùn)算量分析

圖2 算法流程框圖

5 仿真分析

本節(jié)將通過距離-速度-加速度星座圖、參數(shù)估計(jì)均方根誤差曲線及運(yùn)算時(shí)間曲線這3個(gè)方面來驗(yàn)證所提寬帶噪聲雷達(dá)參數(shù)估計(jì)算法的有效性和優(yōu)勢(shì)。參數(shù)設(shè)置為：假設(shè)雷達(dá)信號(hào)為連續(xù)波寬帶共軛噪聲組信號(hào)，信號(hào)載頻f0= 10 GHz ，帶寬B=25 MHz，噪聲組信號(hào)時(shí)寬為T=100 ms，采樣頻率Fs= 60 MHz。

圖3給出了多目標(biāo)情況下的距離-速度-加速度星座圖的結(jié)果，圖中目標(biāo)距離分別為｛300，400，500｝ km，速度分別為｛0.75，1.72，1.502｝ km/s，加速度分別為｛200，300，500｝ m/s2，信噪比同為10 dB。從圖3可以看出，本文算法可以對(duì)目標(biāo)的3個(gè)參數(shù)進(jìn)行準(zhǔn)確的估計(jì)，且同時(shí)估計(jì)多個(gè)目標(biāo)。

圖3 距離-速度-加速度星座圖

圖5給出了本文方法與互模糊函數(shù)方法的運(yùn)算時(shí)間比值與采樣數(shù)的關(guān)系圖，其中通道數(shù)為16路，速度搜索1000次，加速度搜索1000次，將信號(hào)時(shí)寬設(shè)置為16.7～167.0 ms（樣本點(diǎn)數(shù)為1× 106～1× 107點(diǎn)）變化，步進(jìn)間隔為16.7 m s，Monte Carlo實(shí)驗(yàn)為100次。可以看出，本文方法的運(yùn)算時(shí)間要明顯低于互模糊函數(shù)的運(yùn)算時(shí)間，且這種優(yōu)勢(shì)隨著采樣點(diǎn)數(shù)的增加愈加顯著，其中當(dāng)時(shí)寬為100.0 m s（6× 106點(diǎn)）時(shí)，互模糊函數(shù)方法的運(yùn)算時(shí)間約為本文方法的3× 108倍。

圖4 參數(shù)估計(jì)的均方根誤差曲線

圖5 運(yùn)算復(fù)雜度對(duì)比

6 結(jié)束語

寬帶噪聲雷達(dá)在機(jī)動(dòng)目標(biāo)探測(cè)時(shí)存在多普勒敏感、多普勒色散、多普勒展寬及多普勒模糊的問題，互模糊函數(shù)的方法可以解決這些問題，但是其運(yùn)算復(fù)雜度高而難以滿足實(shí)時(shí)性要求。本文提出了一種基于共軛噪聲組信號(hào)的機(jī)動(dòng)目標(biāo)參數(shù)估計(jì)方法，該方法通過回波預(yù)處理及FSC運(yùn)算補(bǔ)償回波的伸縮效應(yīng)，無需高計(jì)算復(fù)雜度的時(shí)域重構(gòu)運(yùn)算；而且采用先多普勒估計(jì)后時(shí)延估計(jì)的方式，避免了互模糊函數(shù)的3維搜索過程；此外，所涉及的Zoom-FrFT及FSC運(yùn)算都可通過FFT實(shí)現(xiàn)，利于工程應(yīng)用。仿真結(jié)果表明，本文方法可以顯著減小運(yùn)算復(fù)雜度，并且在高信噪比的情況下?lián)碛信c互模糊函數(shù)方法相近的參數(shù)估計(jì)性能。

［1］ Sun Hong-bo， Lu Y i-long， and Liu Guo-sui. U ltra-w ideband technology and random signal radar: an ideal combination［J］. IEEE Aerospace and Electronic System s Magazine， 2003，18（11）: 3-7.

［2］ 劉國歲， 顧紅， 蘇衛(wèi)民. 隨機(jī)信號(hào)雷達(dá)［M］. 北京: 國防工業(yè)出版社， 2005: 26-30. Liu Guo-su i， Gu Hong， and Su W ei-m in. Random Signal Radar ［M］. Beijing: National Defense Industry Press， 2005: 26-30.

［3］ Narayanan R M and Dawood M. Doppler estimation using a coherent ultraw ide-band random noise radar［J］. IEEE Transactions on Antennas and Propagation， 2000， 48（6）: 868-878.

［4］ Dawood M， Quraishi N， and A lejos A V. Superresolution Doppler estimation using UWB random noise signals and MUSIC［J］. IEEE Transactions on Aerospace and Electronic Systems， 2013， 49（1）: 325-340.

［5］ M alanow ski M and Kulpa K. Detection of m oving targets w ith continuous-wave noise radar: theory and measurements［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2012， 50（9）: 3502-3509.

［6］ Axelsson S R J. Random noise radar/sodar w ith u ltraw ideband waveform s［J］. IEEE Transactions on Geoscience and Remote Sensing， 2007， 45（5）: 1099-1114.

［7］ Jin Q， W ong K M， and Luo Z Q. The estimation of time delay and Dopp ler stretch of w ideband signals［J］. IEEE Transactions on Signal Processing， 1995， 43（4）: 904-916.

［8］ M alanow ski M， Kulpa K， and Olsen K E. Extending the integration tim e in DVB-T-based passive radar［C］. P roced ings of European Radar Conference （EURAD），M anchester， UK， 2011: 190-193.

［9］ W eiss L G. W avelets and w ideband correlation processing［J］. IEEE Signal Processing Magazine， 1994， 11（1）: 13-32.

［10］ 王釗， 顧紅， 蘇衛(wèi)民， 等. 寬帶高斯型功率譜噪聲信號(hào)相關(guān)特性分析［J］. 電波科學(xué)學(xué)報(bào)， 2013， 28（6）: 1047-1053. W ang Zhao， Gu Hong， Su W ei-m in， et al.. Correlation analysis of noise signal w ith w ideband Gaussian power spectrum［J］. Chinese Journal of Radio Science， 2013， 28（6）: 1047-1053.

［11］ Chen Xiao-long， Guan Jian， Bao Zhong-hua， et al.. Detection and extraction of target w ith m icrom otion in spiky sea clutter via short-tim e fractional Fourier transform［J］. IEEE Transactions on Geoscience and Rem ote Sensing， 2014， 52（2）: 1002-1018.

［12］ Shen Li-ran， Yin Qing-bo， Lu M ing-yu， et al.. Linear FM signal parameter estimation using STFT and FrFT［J］. Chinese Journal of Electronics， 2013， 22（2）: 301-307.

［13］ Ozaktas H M， A rikan O， Ku tay M A， et al.. D igital com putation of the fractional Fourier transform［J］. IEEE Transactions on Signal Processing， 1996， 44（9）: 2141-2150.

［14］ Zhao Xing-hao， Tao Ran， and W ang Yue. A fast algorithm for fractional Fourier transform w ith zoom ing-in ability［C］. Proceedings of 4th International Conference on W ireless Comm unications， Networking and M obile Com puting， Dalian，Ch ina， 2008: 1-4.

［15］ Bi G， M itra S K， and Li S. Samp ling rate conversion based on DFT and DCT［J］. Signal Processing， 2013， 93（2）: 476-486.

王 釗： 男，1985年生，博士生，研究方向?yàn)殡S機(jī)信號(hào)雷達(dá)及寬帶信號(hào)處理.

顧 紅： 男，1967年生，博士，教授，博士生導(dǎo)師，研究方向?yàn)樵肼暲走_(dá)、高速目標(biāo)探測(cè)、M IMO雷達(dá)信號(hào)處理.

蘇衛(wèi)民： 男，1959年生，博士，教授，博士生導(dǎo)師，研究方向?yàn)殛嚵行盘?hào)處理和雷達(dá)成像.

Parameter Estimation of Maneuvering Targets in W ideband Noise Radar Based on Con jugate Noise Group

Wang Zhao①Gu Hong①Su Wei-m in①Chen Zhi-ren①Chen Jin-li②①（Schoo l of E lectronic Engineering & Optoelectronic Techno logy， Nanjing University of Science & Technology， Nanjing 210094， China）
②（College of Electronic and Information Engineering， Nanjing University of Information Science and Technology， Nanjing 210044， China）

The w ideband cross-ambiguity function method is commonly adopted to execute the parameter estimation of w ideband noise radar， but it needs three-dimensional search in distance， velocity and acceleration when dealing w ith maneuvering targets， which takes huge com putation burden. A novel method based on the conjugate noise group is p roposed for addressing the p rob lem of parameter estimation of maneuvering targets. Firstly， the multiple channel is set up according to the echo stretching effect， and the internal signals of the noise group is cut out in fixed length for m ixing in each channel. Then the Doppler phase is estimated with the m ixed signal by Fractional Fourier Transform （FrFT）. The Phase com pensation function is constructed by the Doppler phase and the delay is estimated by Frequency-domain Scale Correlation （FSC） algorithm w ith the com pensated noise group signal. Finally， the range， velocity and acceleration are obtained by the two simultaneous equations of the Dopp ler phase and delay. The proposed m ethod avoids th ree-dimensional search and reconstruction of the echo signal in tim e domain， w hich reduces a large amount of com pu tation com pared to the w ideband cross-ambiguity function method. The method is feasible for real time p rocessing as the whole algorithm can be accom plished by Fast Fourier Transform （FFT）. The effectiveness and superiority of the proposed method are demonstrated by the simulation resu lts.

W ideband noise radar； Conjugate noise group； Fractional Fourier Transform （FrFT）； Frequencydomain Scale Correlation （FSC）

TN957.51

： A

：1009-5896（2015）05-1071-07

10.11999/JEIT140737

2014-06-03收到 ，2014-09-28改回

國家部級(jí)基金，教育部博士點(diǎn)基金（20113219110018），中國航天科技集團(tuán)公司科技創(chuàng)新基金（CASC04-02），國家自然科學(xué)基金（61302188）和江蘇省自然科學(xué)基金（BK 20131005）資助課題

*通信作者：王釗 alife_numone@126.com