秦立龍，王振宇（1.國防科學(xué)技術(shù)大學(xué)電子科學(xué)與工程學(xué)院，湖南長沙410073；.解放軍電子工程學(xué)院通信對抗工程系，安徽合肥30037）

邊際譜和多重分形在調(diào)制模式識別中的應(yīng)用

秦立龍1，2，王振宇2
（1.國防科學(xué)技術(shù)大學(xué)電子科學(xué)與工程學(xué)院，湖南長沙410073；2.解放軍電子工程學(xué)院通信對抗工程系，安徽合肥230037）

為了提高數(shù)字信號調(diào)制模式識別在低信噪比下的正確率，根據(jù)對邊際譜和多重分形理論原理的分析，提出了一種新的基于多重分形理論的特征提取方法。該方法首先引入HHT變換求得樣本的邊際譜，不同調(diào)制模式的邊際譜具有明顯的差異性，可以利用分形的方法提取邊際譜的分形維數(shù)作為調(diào)制識別的特征參數(shù)。最后利用支持向量機分類器進行信號的分類識別。并在求解支持向量機優(yōu)化問題中，利用通用的粒子群算法確定了最優(yōu)系數(shù)。計算機仿真研究證明，新方法提取的特征能有效地提高識別正確率，具有較好的工程應(yīng)用性。

調(diào)制識別；邊際譜；分形理論；支持向量機

數(shù)字信號調(diào)制模式識別能夠在未知調(diào)制信息或有干擾的條件下，正確識別出通信信號調(diào)制的模式，并進一步為解調(diào)器選擇相應(yīng)的解調(diào)算法提供依據(jù)。調(diào)制模式識別在軍事和民用中有著重要的研究前景和應(yīng)用價值。在軍事領(lǐng)域［1］，數(shù)字通信信號調(diào)制模式識別是敵對雙方進行通信偵察和干擾的前提，一旦明確了敵方通信系統(tǒng)的調(diào)制模式，就可以解調(diào)出敵方信號，獲得有用的情報信息，從而為制定偵察與反偵察、干擾與反干擾策略提供有力依據(jù)，最終實現(xiàn)通信對抗；在民用領(lǐng)域［2］，政府有關(guān)部門可以利用調(diào)制模式識別進行信號確認、干擾識別和頻譜監(jiān)測等無線電管理工作，以防止用戶對無線頻譜的非法利用和干擾，保證合法通信的正常進行。

調(diào)制模式識別一般采用統(tǒng)計模式識別方法和決策論方法。目前應(yīng)用較多的為統(tǒng)計模式識別方法［3］?；诮y(tǒng)計模式識別的方法可以分成2個部分：特征提取和分類器設(shè)計，其中好的特征是決定分類正確率的關(guān)鍵因素。現(xiàn)有的特征提取方法有基于時頻分析的方法［4］、基于累積量的方法［5］、基于譜相關(guān)法方法［6］等。然而以上很多通用算法都是在信噪比較高時才能對信號進行識別，在信噪比降低到10 dB以下時識別率明顯下降，算法有效性會大大降低。王建新等［7］利用星座圖的方法得到了識別MASK、MPSK、MQAM信號的特征參數(shù)，此方法提取的特征對噪聲不敏感，在低信噪比下依然有著很高的識別率，但是識別的信號種類較少。為了滿足識別要求，研究抗噪聲性能好、識別信號類型多、分類性能優(yōu)的調(diào)制模式自動識別方法將很有意義。

近些年來，隨著分形理論研究的深入，分形方法在許多基礎(chǔ)的信號和圖像處理中取得了目前最好的效果［8］。不同信號的分形維數(shù)不同，因此分形維數(shù)中包含了一定的分類信息，可以應(yīng)用于模式識別領(lǐng)域中。呂鐵軍等［9］通過對信號調(diào)制規(guī)律的研究，首次提出從信號的分形維數(shù)中提取特征參數(shù)，并在識別CW、MPSK、MFSK中取得了很好的效果。郭雙冰［10］通過對調(diào)制信號進行小波變換，利用變換后的小波系數(shù)提取了2種分形特征，取得了較高的識別率。但是經(jīng)典的小波分析是基于傅里葉分析的，因此具有傅里葉分析的不足：一方面只能給出線性情況的解釋，另一方面表現(xiàn)在小波基被選定后，必須用它來分析所有數(shù)據(jù)，因此不具有自適應(yīng)性。針對通信信號調(diào)制模式的特點，本文提出一種新的特征提取方法。該方法首先利用Hilbert?Huang Transform（HHT）提取數(shù)字信號的邊際譜，利用分形的方法提取特征，最后利用SVM進行調(diào)制模式的分類識別。

1 HHT變換和分形原理

1.1 HHT變換

1998年N.E.Huang等［11］提出了信號時頻分析方法HHT，由于其具有的超自適應(yīng)性，非常適用于非線性和非平穩(wěn)信號的處理，因此非常適于處理具有典型非平穩(wěn)性的數(shù)字調(diào)制信號。HHT理論中主要是經(jīng)驗?zāi)Ｊ椒纸猓‥MD）和Hilbert變換，具體實現(xiàn)過程如下：

不妨設(shè)待識別信號為x（t），首先通過EMD［12］方法將x（t）分解為n個內(nèi)蘊模式分量IMF和一個剩余分量之和：

式中：ci（t）是第i個內(nèi)蘊模式分量，rn（t）為剩余分量。EMD的目的是通過對非線性非平穩(wěn)信號的分解獲得一系列表征信號特征時間尺度的本征模態(tài)函數(shù)（IMF），使得各個IMF是窄帶信號，從而可以進行Hilbert分析。通過式（1）對每一個分量IMF進行Hilbert變換：

式中：P表示柯西主值。然后通過式（2）和式（3）計算每一層信號的瞬時頻率ωi（t）和瞬時幅值ai（t）。

可見，由Hilbert變換得到的幅值和頻率都是時間的函數(shù)，把振幅顯示在頻率－時間的平面上，就可以得到Hilbert幅值譜：

H（ω，t）精確地描述了信號的幅值隨時間和頻率的變化規(guī)律。將H（ω，t）對時間積分，就得到Hilbert邊際譜：

式中：h（ω）表明單位頻率內(nèi)的幅度分布（或者能量分布），代表著整個數(shù)據(jù)段幅度概率分布的累加。

圖1 不同信號的邊際譜示意圖Fig.1 Themarginal spectrum s of differentmodulation m odes

如圖1所示，不同的調(diào)制模式具有不同的邊際譜。因此可以利用分形原理提取目標信號邊際譜的分形特征。

1.2 分形原理

分形是一門嶄新的學(xué)科，是一個研究和處理自然與工程中不規(guī)則圖形的強有力的數(shù)學(xué)理論工具，為處理非線性的系統(tǒng)問題提供了新思路和新方法。分形維數(shù)是分形理論中的重要參數(shù)，它有多種定義和計算方法。Hausdorff維數(shù)是分形理論中一種最基本的分形維數(shù)，但Hausdorff維數(shù)難以計算，實際應(yīng)用中較多采用盒維數(shù)和信息維數(shù)［13］。

設(shè)待測集合A是Rn中的一個非空子集，對于每個正數(shù)ε，令N（A，ε）表示能夠覆蓋A的、邊長為ε的最小的盒子數(shù)目，定義集合A的上下盒維為

實際計算中可以取一系列不同的ε值計算得到對應(yīng)的N（A，ε）值，通過式（4）做最小二乘法線性回歸，得到直線的斜率為信號的盒維數(shù)：

由定義能夠得出，盒維數(shù)沒有反映出空間上的分布情況，只是體現(xiàn)了集合的尺度問題，而信息維數(shù)恰好彌補了這一點。信息維數(shù)的定義如下：

設(shè)｛Aj｝（j＝1，2，…，K）是集合A的一個有限ε?格形覆蓋，pj的含義是A的元素落在Aj中的概率。那么就有信息熵：

假設(shè)信息熵有如下關(guān)系：

那么就說DIA()是A的信息維數(shù)。

計算時，可通過粗視化變換得到信息維數(shù)：

式中：樣本間隔λ＝1／fs，fs是采樣頻率。

上述的盒維數(shù)和信息維數(shù)屬于單分形維數(shù)。有時面對自然界中復(fù)雜的信號，僅用單分形維數(shù)是不夠的，這時候就必須同時使用多個維數(shù)來描述它，才能全面地刻畫其特征，因此引入多重分形［14］的概念，對分形參數(shù)進行擴展。多重分形彌補了分形理論的不足，下面對多重分形的原理進行討論。

如果χd表示d維歐氏空間，χd的任何子集都用X來標記。使用數(shù)學(xué)迭代或者遞推的方法來分割X，并將每個部分的測度μ設(shè)定為一個不變的值，如果α表示一個與劃分相關(guān)聯(lián)的參數(shù)，第n次劃分之后X的子集用Xn(α)來表示，如果Xn＝linm∞Xn(α)

→表示的是一個分形集，由此給出一個測度空間的二元組定義表示式（X，μ），那么就能夠把（X，μ）的其中一個子集用Xn來表示。如果在上述劃分方式下，能夠以幾個分形子集的并集來代替（X，μ）的分形子集，而且其各自都有屬于自己的分維數(shù)，那么就能夠把這個分形集叫做多重分形。假如把（X，μ）按照一定的尺度δ來進行分割，定義單元測度為μα，那么μα和δ之間的關(guān)系為

式中：α是h?lder指數(shù)。對于概率測度為μα的子集Xα的任何δ?覆蓋{Ui}(i∈N)，也就是0＜diam Ui＜δ（diam代表直徑）有

分析可得，當(dāng){Uj}是尺寸為δ的盒子時，若在[α，α＋dα]范圍內(nèi)概率測度為μ的單元個數(shù)（即覆蓋盒子數(shù)）為N(α)，則Xα，μα)＝N(α)δr。因此當(dāng)N(α)～δ－f(α)時，)有限。由此可以α，μα得到

稱fα()為多重分形譜。在α～fα()坐標系中，fα()是一單峰圖像，一個復(fù)雜的分形體，能夠看到它是由一組取值不同的α值所代表的子集構(gòu)成，而fα()所反映的是這些子集的一些局部的分形特性。

如圖2所示，不同調(diào)制方式邊際譜的多重分形譜的譜寬Δα＝αmax－αmin不同，可以作為調(diào)制識別的特征之一。

圖2 不同調(diào)制信號邊際譜的多重分形譜Fig.2 Themultifractal spectrums of differentmodulation modes

綜上所述，盒維數(shù)所體現(xiàn)的是集合的尺度的大小，而信息維數(shù)所體現(xiàn)的是分形集的分布情況，多重分形的譜寬能夠從局部出發(fā)來研究分形集的整體結(jié)構(gòu)和特征，因此特征向量（由信息維數(shù)、分形譜以及盒維數(shù)三者構(gòu)成）所體現(xiàn)的特征有頻率、相位以及幅度的變動情況，在進行相關(guān)分類的時候以上這些特征屬于有效特征集。綜合上述原因，選取向量T：

2 特征提取及識別算法

為了驗證算法在特征提取中的有效性，首先利用MATLAB產(chǎn)生了7種最常用數(shù)字信號的仿真數(shù)據(jù)，這7種數(shù)字調(diào)制信號為：2ASK、4ASK、2FSK、4FSK、2PSK、4PSK、MSK，主要調(diào)制參數(shù)為：載頻fc＝10 kHz，采樣頻率fs＝100 kHz，MSK碼元速率fb＝4× fc／5，其他信號碼元速率為fb＝2 000 B，頻率偏移量△f＝2 000 Hz。為了消除信道差異，還必須對仿真產(chǎn)生的數(shù)字調(diào)制信號進行功率歸一化處理。

通過上述步驟產(chǎn)生的仿真數(shù)據(jù)在使用前，先將連續(xù)的長數(shù)據(jù)分割成便于處理的數(shù)據(jù)段，每段有2 048個采樣數(shù)據(jù)，這一過程稱為分幀。

對7種調(diào)制信號加入信噪比為0～20 dB的噪聲，并且以間隔1 dB變化每種信號，各提取150個特征（50個特征用于后文支持向量機的學(xué)習(xí)，100個特征用于測試分類正確率），計算特征參數(shù)的均值。計算機仿真結(jié)果如圖3～5所示：

圖3 不同調(diào)制信號邊際譜盒維數(shù)均值Fig.3 Themean value of box dimensions

圖4 不同調(diào)制信號邊際譜信息維數(shù)均值Fig.4 Them ean value of information dimensions

圖5 不同調(diào)制信號邊際譜多重分形譜譜寬均值Fig.5 Them ean value ofm ultifractal spectrum w idths

由圖3～5可得，總體上，待識別信號的特征均值向量差別較大（各參數(shù)曲線間距離較大），所以提取的分類特征能夠?qū)崿F(xiàn)信號的調(diào)制模式識別。為了更直觀地說明特征包含各類調(diào)制信號的差別信息，將各信號提取的150個特征作如圖6所示的特征空間分布圖。

圖6 不同調(diào)制信號特征分布圖（SNR＝8 dB）Fig.6 The distribu tion of the features（SNR＝8 dB）

由圖6可以看出，2ASK和4ASK信號的特征混疊在一起，理論上無法加以區(qū)別，這是因為2ASK和4ASK信號的邊際譜分布幾乎完全相同，但是通過觀察及仿真研究，由于MASK信號有相位以及振幅包絡(luò)的變化，而其他調(diào)制信號的波形只有相位的變化，振幅包絡(luò)是恒定的。因此直接提取數(shù)字調(diào)制信號的幅度方差：

圖7 調(diào)制識別模型Fig.7 Themode of the digitalmodulation recognition

除此之外還注意到僅2PSK信號與其他信號的特征略有重疊，其余各類信號特征聚集分界線非常明顯，進一步說明了本文算法的有效性。總之，以上提取的特征提供了待識別信號的差別信息，可作為本文原始特征集的元素。在提取上述特征的基礎(chǔ)上，可以進一步進行調(diào)制模式識別。

3 識別性能分析

本文使用支持向量機分類器（support vector machine，SVM）［15］測試算法的有效性。支持向量機是根據(jù)統(tǒng)計學(xué)習(xí)理論提出的一種機器學(xué)習(xí)方法。它不但結(jié)構(gòu)簡單，能夠很好地解決神經(jīng)網(wǎng)絡(luò)分類器中的傳統(tǒng)問題，而且其有著很強的魯棒性與很好的推廣性，能夠很好地處理高維模式識別、小樣本以及非線性等問題。SVM的本質(zhì)即為求式（6）的對偶問題：

式中：K（xi，xj）為核函數(shù)，實現(xiàn)了從輸入空間到高維空間的映射。由于徑向基支持向量機［16］通用性較好，仿真中使用此類向量機，其核函數(shù)表達式為

在實際應(yīng)用中，需要事先確立2個參數(shù)（懲罰因子c和式（7）中的核參數(shù)σ），不同參數(shù)值，向量機的分類性能不同。利用粒子群算法［17］尋找最優(yōu)參數(shù)的流程如下：

1）在搜索空間中產(chǎn)生m個粒子，隨機初始化粒子位置和初始速度；

2）根據(jù)不同的參數(shù)σ、c，計算各粒子的適應(yīng)度值，并保存各粒子的個體最優(yōu)值和粒子群全局最優(yōu)值；

3）根據(jù)式（8）計算自適應(yīng)慣性權(quán)重，wmax、wmin表示w的最大值和最小值，tmax表示最大迭代步數(shù)。再根據(jù)式（9）、（10）更新粒子的位置和速度：

4）迭代次數(shù)加1，判斷是否滿足最大迭代次數(shù)。若滿足，算法終止并輸出最優(yōu)參數(shù)值和正確率，否則，轉(zhuǎn)向2）。

最后得到最優(yōu)參數(shù)值下信號的正確率如表1。

表1 不同調(diào)制模式分類正確率Table 1 The classification accuracy of differentmodulation modes／%

如表1所示，計算機仿真結(jié)果表明：無需任何先驗知識，在信噪比為0～20 dB的情況下，基于邊際譜和分形理論的特征識別方法在低信噪比下就能夠得到較高的調(diào)制模式識別正確率。

由表1還可以直觀地觀察到，本文的方法對于高斯白噪聲具有很強的抗噪性，其有效性遠在傳統(tǒng)方法之上，其原因在于通信信號的不規(guī)則程度主要取決于調(diào)制類型，而噪聲對之影響較小，反映到特征上，就是分形特征對噪聲不敏感，具有很強的抗噪性［9］，使得最后得到的特征參數(shù)的魯棒性得到提高。

4 結(jié)束語

本文針對7種常用的數(shù)字調(diào)制信號進行了仿真，仿真結(jié)果表明：基于邊際譜和分形理論提取的特征值具有更好的特征均值和方差，提高了識別的正確率，充分說明新算法有一定的有效性。但是該算法也有一定的局限性，IMF分量的邊界處有些失真，這是分解中使用三次樣條插值擬合時，不穩(wěn)定的端點極值會出現(xiàn)擺動的情況，而且會在分解的過程中影響到其他的數(shù)據(jù)，進而影響了系統(tǒng)的識別率，這就是EMD方法中的端點效應(yīng)問題。如何克服端點效應(yīng)問題也是今后研究的方向之一。HHT還是一種新興的非平穩(wěn)信號處理技術(shù)，隨著研究的深入，會有越來越多的學(xué)者加入到這一新興方向的研究隊伍中，新算法也層出不窮，相信以上問題以后會得到解決。

［1］羅吉.數(shù)字調(diào)制信號識別算法研究［D］.成都：電子科技大學(xué)，2008：1?2.LUO Ji.The algorithm research of digitalmodulation recog?nition［D］.Chengdu：University of Electronic Science and Technology，2008：1?2.

［2］劉娟.數(shù)字通信信號自動調(diào)制識別算法研究［D］.重慶：重慶大學(xué)，2009：2?3.LIU Juan.The algorithm research ofmodulation recognition of digital communication signal［D］.Chongqing：Master Dissertation of Chongqing University，2009：2?3.

［3］DOBRE O A，ABDI A，BAR?NESS Y.Blind modulation classification：a concept whose time has come［C］／／IEEE Sarnoff Symphosium.Princeton，USA，2005：223?228.

［4］HASSAN K，DAYOUB I，HAMOUDAW.Automaticmodu?lation recognition usingwavelet transform and neuralnetwork［C］／／Proceedings of the 9th International Conference on Intelligent Transport Systems Telecommunications.Chicago，USA，2009：234?238.

［5］呂新正，魏平，肖先賜.利用高階累積量實現(xiàn)數(shù)字調(diào)制信號的自動識別［J］.電子對抗技術(shù)，2004，19（6）：3?6.LüXinzheng，WEIPing，XIAO Xianci.Automatic identifi?cation of digitalmodulation signals using high order cumula?nts［J］.Electronic Warfare Technology，2004，19（6）：3?6.

［6］郭黎利，齊琳，王東凱.軟件無線電中基于譜相關(guān)理論的調(diào)制模式識別技術(shù)［J］.哈爾濱工程大學(xué)學(xué)報，2004，125（6）：799?802.GUO Lili，QILin，WANG Dongkai.Modulation recognition technique using spectral correlation approach in software ra?dio［J］.Journal of Harbin Engineering University，2004，125（6）：799?802.

［7］王建新，宋輝.基于星座圖的數(shù)字調(diào)制方式識別［J］.通信學(xué)報，2004，25（6）：166?173.WANG Jianxin，SONG Hui.Digitalmodulation recognition based on constellation diagram［J］.Journal of China Institu?te of Communications，2004，25（6）：166?173.

［8］林俊，熊剛，王智學(xué).基于分形特征提取的調(diào)制信號識別［J］.彈箭與制導(dǎo)學(xué)報，2004，24（3）：358?363.LIN Jun，XIONG Gang，WANG Zhixue.Modulation classi?fication based on fractal features［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2004，24（3）：358?363.

［9］呂鐵軍，郭雙冰，肖先賜.調(diào)制信號的分形特征研究［J］.中國科學(xué)：E輯，2001，31（6）：508?513.LüTiejun，GUO Shuangbin，XIAO Xianci.Modulation sig?nals research based on fractal features［J］.Science in Chi?na：Ser E，2001，31（6）：508?513.

［10］郭雙冰.基于小波和分形理論的調(diào)制信號特征提取方法研究［J］.信號處理，2005，21（3）：111?115.GUO Shuangbing.A method of feature extraction formodu?lated signals based wavelet and fractal theory［J］.Signal Processing，2005，21（3）：111?115.

［11］HUANGN E，SHEN Z，LONG SR，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non?stationary time series analysis［J］.Proceedings of the Royal Society：Lond.A，1998，454：903?995.

［12］趙知勁，尹霆，郭殿武.一種基于HHT的數(shù)字調(diào)制信號識別方法［J］.微波學(xué)報，2006，22（3）：66?70.ZHAO Zhijin，YIN Ting，GUO Dianwu.A method to Iden?tify digital modulation signals based on the Hilbert?Huang transform［J］.Journal of Microwaves，2006，22（3）：66?70.

［13］張旻，鐘子發(fā)，王若冰.通信電臺個體識別技術(shù)研究［J］.電子學(xué)報，2009，37（10）：2125?2129.ZHANGMin，ZHONG Zifa，WANG Ruobing.Research on the technique of individual communication transmitter iden?tification［J］.Acta Electronica Sinica，2009，37（10）：2125?2129.

［14］李彤，商朋見.多重分形在掌紋識別中的研究［J］.物理學(xué)報，2007，56（8）：4393?4400. LI Tong，SHANG Pengjian.A multifractal approach to palmprint recognition［J］.Acta Physica Sinica，2007，56（8）：4393?4400.

［15］何曉琴，白勇，冉啟陽.基于支持向量機的Freegate軟件流量檢測研究［J］.計算機科學(xué)，2012，39（11）：89?90.HE Xiaoqin，BAI Yong，RAN Qiyang.Freegate software flow test research based on support vector machine［J］.Computer Science，2012，39（11）：89?90.

［16］張翔，田金文，肖曉玲，等.支持向量機及其在醫(yī)學(xué)圖像分類中的應(yīng)用［J］.信號處理，2004，20（2）：208?212.ZHANG Xiang，TIAN Jinwen，XIAO Xiaoling，et al.Sup?port vector machine and its application in medical images classification［J］.Signal Processing，2004，20（2）：208?212.

［17］KENNEDY J，EBERHART R C.Particle swarm optimiza?tion［C］／／Proceedings of IEEE InternationalConference on Neural Networks.Perth，Australia，1995：1942?1948.

［18］易蘇軍.Hilbert?Huang變換及其在信號處理中的應(yīng)用［D］.大連：大連理工大學(xué)，2006：9?10.YISujun.Hilbert?Huang transform and its application in signal processing［D］.Dalian：Dalian University of Tech?nology，2006：9?10.

秦立龍，男，1988年生，博士研究生，主要研究方向為調(diào)制模式識別。

王振宇，男，1956年生，副教授，主要研究方向為信號處理。

M arginal spectrum and multifractal theory and its app lication in m odulation recognition

QIN Lilong1，2，WANG Zhenyu2
（1.School of Electronic Science and Engineering，National University of Defence Technology，Changsha410073，China；2.Department of Communication Countermeasure Engineering，Electronic Engineering Institute，Hefei230037，China）

Through the analysis of the marginal spectrum and multifractal theory，a new feature extraction method based on multifractal theory was proposed to improve the accuracy of the digitalmodulation recognition under the low signal?to?noise ratio.First，the Hilbert?Huang transform was put forward to obtain themarginal spectrum of the samples.There are differences among differentmodulation modes.The fractal dimensions of the sample after Hil?bert?Huang transform were calculated by the fractalmethod.Next，the featurewas extracted.Finally，the identifica?tion task was solved by using SVM classification machine.In order to determine the optimal coefficient of the sup?port vectormachine，a universal particle swarm optimization algorithm was used.The computer simulation results showed that the performance of this feature extracted by the new algorithm efficiently improves the accuracy ofmod?ulation recognition and could be feasible to use in engineering applications.

modulation recognition；marginal spectrum；fractal theory；support vectormachine

TP18；TN911.7

A

1673?4785（2014）06?0756?07

秦立龍，王振宇.邊際譜和多重分形在調(diào)制模式識別中的應(yīng)用［J］.智能系統(tǒng)學(xué)報，2014，9（6）：756?762.

英文引用格式：QIN Lilong，WANG Zhenyu.M arginal spectrum and m ultifractal theory and its app lication in modulation recogni?tion［J］.CAAITransactions on Intelligent System s，2014，9（6）：756?762.

10.3969／j.issn.1673?4785.201301031

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1673－4785.201301031.htm l

2013?01?16.

日期：2014?09?30.

國家自然科學(xué)基金資助項目（61040007）.

秦立龍.E?mail：tank2908989＠163.com.