基于貝葉斯序貫推理的自適應(yīng)調(diào)制識別算法

付俊強(qiáng)，李 蓉，趙成林，李 斌

（1.北京郵電大學(xué)信息與通信工程學(xué)院，北京100876；2.國家無線電監(jiān)測中心，北京100037）

基于貝葉斯序貫推理的自適應(yīng)調(diào)制識別算法

付俊強(qiáng)1，李 蓉2，趙成林1，李 斌1

（1.北京郵電大學(xué)信息與通信工程學(xué)院，北京100876；2.國家無線電監(jiān)測中心，北京100037）

提出了一種時(shí)變衰落信道下的自適應(yīng)調(diào)制識別算法，設(shè)計(jì)出一種新的動態(tài)狀態(tài)空間模型，來刻畫信號調(diào)制方式與時(shí)變信道增益的時(shí)變特性，并引入一階有限狀態(tài)馬爾可夫（finite state Markov channel，F(xiàn)SMC）模型來描述衰落信道；基于上述，新算法采用貝葉斯序貫推理法，充分發(fā)掘利用了其中所隱含的信道動態(tài)相關(guān)特性，實(shí)現(xiàn)對調(diào)制方式和時(shí)變信道增益的聯(lián)合估計(jì)。仿真結(jié)果表明，新算法性能相比于傳統(tǒng)ALRT算法有極大提升，且增加采樣點(diǎn)數(shù)或者降低多普勒頻移值都會使算法性能得到改善。

自適應(yīng)調(diào)制識別；時(shí)變衰落信道；動態(tài)轉(zhuǎn)移；貝葉斯序貫推理；聯(lián)合估計(jì)

0 引 言

調(diào)制識別作為下一代通信技術(shù)［1］中鏈路自適應(yīng)傳輸?shù)年P(guān)鍵，目標(biāo)是識別當(dāng)前接收信號所采用調(diào)制方式，以及一些細(xì)節(jié)調(diào)制特征參數(shù)。調(diào)制識別技術(shù)可以大致分為兩類。一類是基于特征的，可以歸納為兩步，第一步是特征提取，例如高階累積量［2］、小波變換［3］等；第二步是分類器，諸如盲聚類［4］、神經(jīng)網(wǎng)絡(luò)［5］和支持向量機(jī)［6］。另一類是基于最大似然［79］的，利用接收信號的統(tǒng)計(jì)似然特征，在假設(shè)信號先驗(yàn)信息符合均勻分布的前提下，降低了信號錯(cuò)誤分類概率。

最大似然估計(jì)的一個(gè)關(guān)鍵問題需要引起重視，它嚴(yán)重依賴完整的信道狀態(tài)信息。為了解決未知信道參數(shù)問題，文獻(xiàn)［10］提出了一種平均似然檢測（average likelihood ratio test，ALRT）算法。

考慮到移動通信環(huán)境，自適應(yīng)調(diào)制編碼需應(yīng)用到時(shí)變衰落信道［11-12］中，因此時(shí)變衰落信道下的自適應(yīng)調(diào)制識別也將帶來兩個(gè)挑戰(zhàn)。第一，調(diào)制方式是隨著信道的變化而不斷改變的，它并不是完全獨(dú)立的，因此調(diào)制方式與信道增益是相關(guān)的，先驗(yàn)信息不再符合均勻分布而是隨著時(shí)間變化。第二，時(shí)變衰落信道也給接收信號帶來很多不確定性，使得似然函數(shù)變得不精確。ALRT算法的假設(shè)前提顯然不適合時(shí)變信道，而且不能實(shí)時(shí)跟蹤時(shí)變衰落信道。

作為另一種克服時(shí)變衰落信道的方案，亦可借助于先驗(yàn)性的導(dǎo)頻序列估計(jì)并跟蹤時(shí)變衰落信道，上述方案已為LTE網(wǎng)絡(luò)與IEEE 802.11 WLAN系統(tǒng)所廣泛采用。然而，進(jìn)一步考慮未來更加復(fù)雜的應(yīng)用場景，譬如認(rèn)知無線電和D2D異構(gòu)形式網(wǎng)絡(luò)，則很難在兩個(gè)非協(xié)作設(shè)備間配置先驗(yàn)導(dǎo)頻序列執(zhí)行信道估計(jì)。此類場景下，只能采用無數(shù)據(jù)輔助或者盲估計(jì)方式來應(yīng)對時(shí)變衰落信道。

鑒于上述時(shí)變衰落信道給調(diào)制識別所帶來的挑戰(zhàn)，本文提出了一種基于貝葉斯序貫推理的自適應(yīng)調(diào)制識別算法，旨在解決傳統(tǒng)ALRT算法與導(dǎo)頻輔助估計(jì)方案在未來復(fù)雜應(yīng)用場景中所面臨的問題。基于新提出的動態(tài)狀態(tài)空間模型［13］（dynamic state-space model，DSM），充分利用時(shí)變信道的先驗(yàn)轉(zhuǎn)移特性，設(shè)計(jì)了一種新的貝葉斯序貫推理算法，通過對信道狀態(tài)增益和調(diào)制方式聯(lián)合估計(jì)來達(dá)到自適應(yīng)調(diào)制識別的效果。通過實(shí)驗(yàn)仿真，將新方案與現(xiàn)有ALRT方案和導(dǎo)頻輔助信道估計(jì)方案進(jìn)行了對比，驗(yàn)證了所提方案的有效性。

1 系統(tǒng)模型

1.1 自適應(yīng)調(diào)制

假設(shè)時(shí)變衰落信道服從萊斯衰落信道特性，信道衰落增益服從萊斯分布，如式（1）所示：

式中，x表示信道衰落增益；σ2R表示萊斯分布的方差；a表示主徑幅度；I0（·）為零階貝塞爾函數(shù)。為了實(shí)現(xiàn)自適應(yīng)調(diào)制，將信道增益劃分成K個(gè)離散狀態(tài)，須滿足K≥3，根據(jù)信道增益來進(jìn)行不同調(diào)制方式的選擇，如圖1所示，其中圖1右側(cè)為萊斯概率密度分布圖。

圖1 時(shí)變衰落信道下的自適應(yīng)調(diào)制模型

1.2 DSM模型

為了更好地應(yīng)對時(shí)變衰落信道下的自適應(yīng)調(diào)制帶來的挑戰(zhàn)，建立一種全新的時(shí)變衰落信道下自適應(yīng)調(diào)制識別的DSM模型如式（2）～式（4）所示：

式中，αn表示第n個(gè)時(shí)刻的信道衰落幅度增益，也按照特定狀態(tài)轉(zhuǎn)移函數(shù)H（·）進(jìn)行更新；cn表示第n個(gè)時(shí)刻的信號調(diào)制方式，sn表示第n個(gè)時(shí)刻發(fā)送的原始星座圖信號點(diǎn)，按照一定的映射函數(shù)T（·）進(jìn)行確定。下列分析對以上3個(gè)動態(tài)方程式（2）～式（4）作更進(jìn)一步解釋：

（1）狀態(tài)轉(zhuǎn)移方程H（·）：R1→R1，它基于一階馬爾科夫鏈建模，描述第n個(gè)時(shí)隙下的時(shí)變衰落信道增益αn→R1的動態(tài)轉(zhuǎn)移過程。

（2）映射函數(shù)T（·）：R1→R1，它給出了未知調(diào)制方式cn∈R1與信道衰落增益αn之間的映射曲線關(guān)系。

（3）觀測方程Z（·）［14］：RM→RM，它描述了兩個(gè)隱含狀態(tài)αn與cn和觀測向量zn→RM之間的關(guān)系，而觀測量由M個(gè)離散采樣觀測值zn，m（m＝0，1，…，M－1）組成。

基于DSM模型，進(jìn)行如下假設(shè)：

假設(shè)1 觀測過程中的噪聲wn，m∈R1是n時(shí)刻下第m個(gè)采樣點(diǎn)的獨(dú)立同分布的零均值方差為σ2加性高斯白噪聲，并和αn與cn相互獨(dú)立。

假設(shè)2 時(shí)變衰落信道是慢衰落過程，其信道增益與傳輸時(shí)隙La有關(guān)，而La與最大多普勒頻移值fd成反比，表示為La∞1／fd。La與感知周期Ts之間須滿足La是Ts的整數(shù)倍，如圖2所示。另外，它與采樣點(diǎn)數(shù)之間是相互獨(dú)立。

圖2 信道感知模型框圖

1.3 時(shí)變衰落信道模型

為了有效描述信道的動態(tài)轉(zhuǎn)移特性，主要考慮一階有限狀態(tài)馬爾可夫信道（finite state Markov channel，F(xiàn)SMC）［15］模型。首先假設(shè)狀態(tài)概率向量π＝［π0，π1，…，πK－1］T，其中πk＝Pr（αn＝Ak），Ak∈A∈R1，k∈｛0，1，…，K－1｝表示第k個(gè)信道狀態(tài)下的平均增益，A表示各信道平均增益集合，穩(wěn)定狀態(tài)下需要滿足PTπ＝π［16］，PK×K＝｛Pk1→k2，i，j∈0，1，…，K－1｝可以由離散信道狀態(tài)轉(zhuǎn)移矩陣（probability transfer matrix，PTM）可得

矩陣中的各個(gè)元素Pk1→k2表示從狀態(tài)k1到狀態(tài)k2的轉(zhuǎn)移概率。

在FSMC模型中，連續(xù)信道需要被劃分成多個(gè)無重疊的離散區(qū)間，Ak由一組信噪比門限值s＝｛s1，s2，…，sK｝決定，或者是一組信道增益值g＝｛g1，g2，…，gK｝s。

信道狀態(tài)增益αn在第k個(gè)信道區(qū)間里的穩(wěn)態(tài)概率由下式可得：

式中，f（αn－1；αn）表示時(shí)變信道的聯(lián)合概率密度分布。

對于一階FSMC，信道狀態(tài)轉(zhuǎn)移的可能性局限于自身和它相鄰的狀態(tài)，因此對于｜k1－k2｜＞1可以得到Pk1→k2＝0。實(shí)際上，可以通過計(jì)算區(qū)間跨越率，得到Pk1→k2≈Nk2／Rk1，其中Rk1＝πk2／T表示在狀態(tài)k1下每秒內(nèi)所含的感知周期數(shù)，其中T為符號周期。

1.4 似然觀測

接收端對自適應(yīng)調(diào)制信號進(jìn)行似然觀測，它的原理就是在特定觀測時(shí)間窗內(nèi)，對接收信號進(jìn)行采樣，將觀測量代入似然函數(shù)計(jì)算似然值?；鶐盘柨梢员硎緸?/p>

式中，i為虛部單位；Es為發(fā)射信號的平均符號能量；sn，m（l）和φn，m（l）分別為信號的幅度和相位；g（·）為成型濾波器；Cl為第l種調(diào)制方式。經(jīng)過模數(shù)轉(zhuǎn)換為

文獻(xiàn)［7］中提出的最大似然算法（maximum likelihood，ML）通常假設(shè)先驗(yàn)信息服從均勻分布，其似然函數(shù)的對數(shù)形式如下所示：

式中，M為一個(gè)感知周期內(nèi)的采樣點(diǎn)數(shù)；zn，m表示第n時(shí)刻對應(yīng)的采樣值；｜Cl｜表示第l種調(diào)制方式的星座點(diǎn)總數(shù)。

1.5 映射規(guī)則

發(fā)射端根據(jù)信道增益Ak∈A→R1動態(tài)地調(diào)整調(diào)制方式cn（l）＝Cl→R1（l＝0，2，…，L－1），所有調(diào)制方式的總數(shù)為L。而接收端對信道增益進(jìn)行聯(lián)合估計(jì)，依據(jù)特定的映射關(guān)系來確定調(diào)制方式。采取一對一映射，其映射關(guān)系式T（·）：R1→R1，它將信道增益與調(diào)制方式進(jìn)行一一對應(yīng)，即L＝K，具體關(guān)系式如下：

2 貝葉斯序貫推理

基于上述系統(tǒng)模型，進(jìn)一步對兩個(gè)隱含系統(tǒng)狀態(tài)（即信道增益和未知調(diào)制方式）進(jìn)行聯(lián)合估計(jì)。估計(jì)檢測算法主要依賴于最大后驗(yàn)概率和序貫推理思想，主要包含以下3個(gè)步驟：

步驟1根據(jù)特定時(shí)間窗內(nèi)的信號計(jì)算似然值；

步驟2聯(lián)合先驗(yàn)信息依據(jù)最大后驗(yàn)概率判決出信道增益；

步驟3最后由映射函數(shù)T（·）確定未知調(diào)制方式。

算法流程如圖3所示。

圖3 接收端調(diào)制識別流程圖

下面針對上述步驟進(jìn)行詳細(xì)介紹：

（1）似然觀測

根據(jù)第1.4節(jié)中已經(jīng)介紹的ML算法中似然函數(shù)式（13）計(jì)算似然值。

（2）對信道增益進(jìn)行貝葉斯序貫估計(jì)

為了有效描述兩個(gè)隱含狀態(tài)，用向量Xn＝｛αn，cn＝T（αn）｝∈R2表示，定義n－1時(shí)刻的后驗(yàn)概率為pn－1｜n－1（Xn－1｜z1∶n－1），那么pn｜n（Xn｜z1∶n）表示n時(shí)刻的后驗(yàn)概率。

根據(jù)最大后驗(yàn)概率準(zhǔn)則，通過式（15）估計(jì)出n－1時(shí)刻的信道狀態(tài)。

相位估計(jì)式（15a）［17］與幅度估計(jì)式（15b）是相互獨(dú)立的，其中Mk為相位轉(zhuǎn)移鍵控（phase shift keying，PSK）和正交幅度調(diào)制（quadrature amplitude modulation，QAM）的調(diào)制階數(shù)。利用n－1時(shí)刻的估計(jì)值通過概率轉(zhuǎn)移矩陣（probability transfer matrix，PTM）得到n時(shí)刻的先驗(yàn)轉(zhuǎn)移概率，如下式

結(jié)合之前所得到的n時(shí)刻的似然值，根據(jù)式（17）求出此時(shí)各信道狀態(tài)的后驗(yàn)概率值。

式中，似然函數(shù)φn（·）＝exp（Λ（zn｜sn，m，αn））。

（3）根據(jù)映射規(guī)則確定調(diào)制方式

映射信道狀態(tài)與未知調(diào)制方式呈一一對應(yīng)的關(guān)系，可以根據(jù)第1.5節(jié)中所述的映射規(guī)則得到下式：

3 仿真結(jié)果

通過Matlab仿真實(shí)驗(yàn)給出了萊斯時(shí)變衰落信道下的自適應(yīng)調(diào)制識別與估計(jì)性能。設(shè)萊斯信道均值為1，方差為0.1，信道狀態(tài)個(gè)數(shù)K＝6，發(fā)送端選用二相相移鍵控（binary PSK，BPSK）、正交相位偏移鍵控（quadrature PSK，QPSK）、8PSK、16QAM、32QAM、64QAM等，詳見表1，并給出了自適應(yīng)調(diào)制識別算法與現(xiàn)有的克服信道衰落的導(dǎo)頻輔助估計(jì)算法和ALRT算法性能對比，突出本算法的優(yōu)越性。

表1 映射關(guān)系表

3.1 不同的多普勒頻移信道

仿真實(shí)驗(yàn)過程中，感知周期Ts內(nèi)采樣點(diǎn)個(gè)數(shù)M＝100，多普勒頻移fd＝0.1，0.05，0.02，即La＝10，20，50。兩種算法性能仿真結(jié)果見圖4，在相同條件下，新算法比ALRT算法有更高的準(zhǔn)確識別率，性能得到顯著的提升。在識別率接近100%的信噪比區(qū)域內(nèi)，不同fd下的曲線幾乎重合。對于新算法，降低fd會提高識別率，而ALRT算法與fd沒有明顯的變化。圖5為新算法在不同fd下的信道估計(jì)性能歸一化均方誤差（mean square error，MSE）曲線，MSE值也是隨著fd的減小而不斷變小。

圖4 不同fd值下的識別率曲線圖

3.2 不同的采樣點(diǎn)數(shù)

在其他條件保持不變的前提下，固定fd＝0.05，改變采樣點(diǎn)數(shù)M，M＝60，100，200。實(shí)驗(yàn)仿真得到的算法識別性能結(jié)果如圖6所示，新算法仍比ALRT算法具有更好的性能。而且隨著M值的增大，識別率得到顯著地提高。圖7為新算法在不同M值下的信道估計(jì)性能歸一化MSE曲線，M值的增加同樣使得信道增益估計(jì)性能得到明顯改善。

圖5 新算法在不同fd值下的歸一化MSE曲線圖

圖6 不同M值下的識別率曲線圖

圖7 新算法在不同M值下的歸一化MSE曲線圖

3.3 現(xiàn)有方案對比

現(xiàn)有的導(dǎo)頻輔助估計(jì)方案，發(fā)送端每隔Ts周期發(fā)送固定復(fù)序列［1＋j，1－j，－1－j］，并且采用最小均方誤差估計(jì)（minimum mean square error，MMSE）算法對信道增益進(jìn)行估計(jì)。圖8給出了3種算法的性能對比圖，導(dǎo)頻估計(jì)性能要優(yōu)于新算法2～3 dB，ALRT算法性能最差。然而，需要注意導(dǎo)頻輔助估計(jì)方案需要額外的先驗(yàn)導(dǎo)頻序列；在某些分布式或異構(gòu)網(wǎng)絡(luò)場景中，導(dǎo)頻輔助方案存在一定的應(yīng)用局限性。相比之下，新方案不受先驗(yàn)性導(dǎo)頻序列的影響，因而也具有更廣泛的應(yīng)用場景。

圖8 fd＝0.05，m＝100下3種方案的性能對比圖

最后，借助于不同算法所需的乘法次數(shù)來評估期計(jì)算復(fù)雜度。由之前的算法實(shí)現(xiàn)流程分析得出，新算法的時(shí)間復(fù)雜度約為O（LM（Cl＋log2Mk））；而導(dǎo)頻輔助估計(jì)的時(shí)間復(fù)雜度約為O（LCl），ALRT算法的時(shí)間復(fù)雜度約為O（LM（KCl＋log2Mk））?？梢园l(fā)現(xiàn)，導(dǎo)頻估計(jì)方案具有最低實(shí)現(xiàn)復(fù)雜度，而新算法相對于ALRT方案將更適合實(shí)時(shí)處理要求。

4 結(jié) 論

針對時(shí)變衰落信道的自適應(yīng)調(diào)制識別問題，區(qū)別與傳統(tǒng)的靜態(tài)系統(tǒng)模型，將時(shí)變信道增益建模為一種FSMC模型，并提出一種全新的基于貝葉斯序貫推理的自適應(yīng)調(diào)制識別算法。新方案充分利用了時(shí)變信道狀態(tài)之間的轉(zhuǎn)移記憶特性，并評估特定時(shí)間窗內(nèi)接收信號的似然函數(shù)，最后利用最大后驗(yàn)概率準(zhǔn)則以及序貫推理思想對信道增益進(jìn)行聯(lián)合估計(jì)，進(jìn)而根據(jù)映射準(zhǔn)則確定調(diào)制方式。相比于ALRT算法，新算法可獲得更優(yōu)的識別性能，并為后續(xù)相關(guān)研究提供了一種全新理論框架，將有利于無線通信信號的信號解調(diào)、參數(shù)提取、干擾分析、信息對抗等各方面的研究。

［1］Andrews J G，Buzzi S，Choi W，et al.What will 5G be［J］.IEEE Journal on Selected Areas in Communications，2014，31（6）：1065- 1082.

［2］Liu M Z，Zhao Y，Shi L，et al.Research on recognition algorithm of digital modulation by higher order cumulants［C］∥Proc.of the 4th IEEE International Conference on Instrumentation and Measurement，Computer，Communication and Control（IMCCC），2014：686- 690.

［3］Yin SG，Liu W，Wang C.Based on the singular value decomposition fast wavelet ridge intra pulse modulation recognition［C］∥Proc.of the Applied Mechanics and Materials，2014：4933- 4940.

［4］Wang J X，Zhang L P.Blind modulation recognition algorithm for MQAM signals［J］.Journal of Electronics&Information Technology，2011，33（2）：332- 336.

［5］Albers H E.The regulation of social recognition，social communication and aggression：vasopressin in the social behavior neural network［J］.Hormones and Behavior，2012，61（3）：283- 292.

［6］Kannan R，Ravi S.Second-order statistical approach for digital modulation scheme classification in cognitive radio using support vector machine and k-nearest neighbor classifier［J］.Journal of Computer Science，2013，9（2）：235- 243.

［7］Bai D，Lee J，Kim S，et al.Near ML modulation classification［C］∥Proc.of the IEEE Vehicular Technology Conference（VTC Fall），2012：1- 5.

［8］Xu J L，Su W，Zhou M.Likelihood-ratio approaches to automatic modulation classification［J］.IEEE Trans.on Systems，Man，and Cybernetics，Part C：Applications and Reviews，2011，41（4）：455- 469.

［9］Hameed F，Dobre O，Popescu D C.On the likelihood-based approach to modulation classification［J］.IEEE Trans.on Wireless Communications，2009，8（12）：5884- 5892.

［10］Dobre O A，Abdi A，Bar-Ness Y，et al.Survey of automatic modulation classification techniques：classical approaches and new trends［J］.Communications，2007，1（2）：137- 156.

［11］Abeida H，Al-Harthi M M，Merabtine N，et al.EM and MAP methods for joint path delay and complex gain estimation of a slowly varying fading channel for CPM signals［J］.International Journal of Engineering and Technology，2014，4（5）：261- 270.

［12］Sklar B.Rayleigh fading channels in mobile digital communication systems I：characterization［J］.IEEE Communications Magazine，2002，35（7）：90- 100.

［13］Wang H S，Moayeri N.Finite-state Markov channel-a useful model for radio communication channels［J］.IEEE Trans.on Vehicular Technology，1995，44（1）：163- 171.

［14］Rabiner L R，Juang B H.An introduction to hidden Markov models［J］.IEEE Trans.on Acoustics，Speech，Signal Processing Magazine，1986，3（1）：4- 16.

［15］Sadeghi P，Kennedy R A，Rapajic P B，et al.Finite-state Markov modeling of fading channels-a survey of principles and applications［J］.Signal Processing Magazine，2008，25（5）：57- 80.

［16］Li B，Zhao C L，Sun M W，et al.Spectrum sensing for cognitive radios in time-variant flat-fading channels：a joint estimation approach［J］.IEEE Trans.on Communications，2014，62（8）：2665- 2680.

［17］Dobre O，Hameed F.Likelihood-based algorithms for linear digital modulation classification in fading channels［C］∥Proc. of the IEEE Conference on Electrical and Computer Engineering，2006：1347- 1350.

Sequential Bayesian inference based adaptive modulation recognition algorithm

FU Jun-qiang1，LI Rong2，ZHAO Cheng-lin1，LI Bin1
（1.School of Information and Communication Engineering，Beijing University of Posts and Telecommunications，Beijing 100876，China；2.The State Radio Monitoring Center，Beijing 100037，China）

Under the time-varying fading channel，an adaptive modulation recognition algorithm is presented.A new dynamic state space model is designed to describe time-varying characteristics of modulation schemes and channel gain.A first-order finite state Markov channel（FSMC）model is introduced for the fading channel.On this basis，a new algorithm，which adopts the sequential Bayesian inference method and is proposed to fully exploit the dynamic transfer characteristics of the hidden channel state，achieves joint estimation of modulation and time-varying channel gain.The simulation results prove that performance of the algorithm compared to traditional ALRT algorithms greatly improves，and increasing the number of sampling points or reducing the Doppler shift value can make the performance better.

adaptive modulation recognition；time-varying fading channel；dynamic transfer characteristics；Bayesian sequential inference；joint estimation

TN 911.7

A

10.3969／j.issn.1001-506X.2015.12.30

付俊強(qiáng)（199-2- ），男，碩士研究生，主要研究方向?yàn)榻y(tǒng)計(jì)信號處理、調(diào)制識別。

E-mail：750450554＠qq.com

李 蓉（1984- ），女，工程師，博士，主要研究方向?yàn)闊o線通信、頻譜監(jiān)測與頻譜管理。

E-mail：lirong＠srrc.org.cn

趙成林（1964- ），男，教授，博士，主要研究方向?yàn)樾盘柼幚?、無線通信技術(shù)。

E-mail：clzhao＠bupt.edu.cn

李 斌（198-5- ），男，講師，博士，主要研究方向?yàn)榻y(tǒng)計(jì)信號估計(jì)與檢測、認(rèn)知無線電。

E-mail：binli＠bupt.edu.cn

1001-506X（2015）12-2860-05

2014- 01- 27；

2015- 06- 01；網(wǎng)絡(luò)優(yōu)先出版日期：2015- 07- 28。

網(wǎng)絡(luò)優(yōu)先出版地址：http：／／www.cnki.net／kcms／detail／11.2422.TN.20150728.0941.002.html

國家自然科學(xué)基金（61271180）資助課題