頻率估計(jì)的多段正弦信號(hào)快速頻譜融合算法

肖瑋，涂亞慶，沈艷林，蘇丹，張磊

1.后勤工程學(xué)院后勤信息工程系，重慶 401311

2.涿州綜合倉庫，河北保定 611730

◎信號(hào)處理◎

頻率估計(jì)的多段正弦信號(hào)快速頻譜融合算法

肖瑋1，涂亞慶1，沈艷林1，蘇丹1，張磊2

1.后勤工程學(xué)院后勤信息工程系，重慶 401311

2.涿州綜合倉庫，河北保定 611730

多段正弦信號(hào)頻譜融合法（簡稱“原融合算法”）是提高低信噪比條件下正弦信號(hào)頻率估計(jì)精度的一條有效途徑，具有重要研究意義和應(yīng)用價(jià)值。為滿足雷達(dá)、聲納、電子對抗等實(shí)時(shí)性要求較高的頻率估計(jì)應(yīng)用需求，提出多段正弦信號(hào)快速頻譜融合算法。該方法通過設(shè)計(jì)離散時(shí)間傅里葉變換（Discrete Time Fourier Transform，DTFT）快速算法、降維處理加權(quán)融合頻譜矩陣和1/3主瓣相關(guān)性分析處理等措施來降低算法計(jì)算量，提高實(shí)時(shí)性。重點(diǎn)對上述三項(xiàng)措施的原理進(jìn)行了闡述與分析。計(jì)算量對比和仿真實(shí)驗(yàn)表明，多段正弦信號(hào)快速頻譜融合算法在精度損失極小的前提下，能夠大幅降低計(jì)算量；在信噪比極低的情況下(SNR≤-13 dB)，其性能略優(yōu)于原融合算法。

快速頻譜融合；多段正弦信號(hào)；離散時(shí)間傅里葉變換（DTFT）快速算法；降維；1/3主瓣相關(guān)

1 引言

多段正弦信號(hào)頻譜融合法（簡稱“原融合算法”）是一類新興的正弦信號(hào)頻率估計(jì)法。該方法通過對多段信號(hào)頻譜進(jìn)行加權(quán)融合實(shí)現(xiàn)多段信號(hào)的信息積累，能夠近似成倍延長被測信號(hào)時(shí)寬，從而解決單段信號(hào)包含信息量不足、抗噪性差的問題，從源頭上為提高低信噪比條件下正弦信號(hào)頻率估計(jì)精度創(chuàng)造條件[1-4]。因此其頻率估計(jì)精度較高，除應(yīng)用于頻率估計(jì)外，還可用于LFM（Linear Frequency Modulated）信號(hào)參數(shù)估計(jì)[1]、VCO（Voltage Controlled Oscillator）非線性度校正[5]等方面，具有重要的理論研究意義和工程實(shí)用價(jià)值。

多段正弦信號(hào)頻譜融合法的計(jì)算量主要集中在多段正弦信號(hào)頻譜計(jì)算、最優(yōu)加權(quán)融合頻譜生成和頻域相關(guān)性分析三部分[4]。為滿足雷達(dá)、聲納、電子對抗等實(shí)時(shí)性要求較高的應(yīng)用需求[2，6]，本文通過設(shè)計(jì)離散時(shí)間傅里葉變換（Discrete Time Fourier Transform，DTFT）快速算法、降維處理加權(quán)融合頻譜矩陣和1/3主瓣相關(guān)性分析處理等措施來分別降低上述三個(gè)環(huán)節(jié)的計(jì)算量，提出頻率估計(jì)的多段正弦信號(hào)快速頻譜融合算法（簡稱“快速融合算法”）。

2 算法原理

2.1 DTFT快速算法設(shè)計(jì)

2.1.1 DTFT快速算法計(jì)算量分析

在多段正弦信號(hào)頻譜融合法中，利用DTFT計(jì)算M (M≥2)段單段信號(hào)頻譜Xm(fP)所需的計(jì)算量較大[4]。文獻(xiàn)[7]提出了一種滑動(dòng)DTFT快速算法，該方法的實(shí)質(zhì)是通過不斷增加計(jì)算序列的長度來實(shí)現(xiàn)指定頻率處傅里葉系數(shù)的快速計(jì)算，并未對DTFT算法本身進(jìn)行改進(jìn)，因此不能適用于多段正弦信號(hào)頻譜融合法中序列長度不能增加的情況。

在多數(shù)數(shù)字信號(hào)處理的書籍和文獻(xiàn)中，DTFT是針對離散信號(hào)的[8-9]，一般僅考慮n=0，1，…，N-1的有限長序列，其定義如式（1）所示：

根據(jù)需要確定被測頻率f的取值范圍fscope= [fmin，fmax]、頻率分辨率Δf、采樣頻率fs和輸出點(diǎn)數(shù)K= (fmax-fmin)/Δf+1，代入式（1）得：

其中，k=0，1，…，K-1。

利用布魯斯坦（Bulestein）所提出的式（3）[8-9]為DTFT設(shè)計(jì)快速算法：

則式（2）可以表示為：

由式（4）可知，計(jì)算x(n)的K點(diǎn)DTFT可以通過計(jì)算序列g(shù)(n)和h(n)的線性卷積與表達(dá)式e-jπk2Δf/fs的乘積來實(shí)現(xiàn)。由于g(n)為長度為N的序列，h(n)為一無窮長的序列；為計(jì)算g(n)和h(n)的線性卷積，只需取h(n)在-N+1≤n≤K-1內(nèi)的值即可，因此可以把h(n)看成一個(gè)長度為L=K+N-1的有限長序列。由于時(shí)域計(jì)算線性卷積效率低，當(dāng)循環(huán)卷積的循環(huán)長度大于或等于L+N-1時(shí)，計(jì)算循環(huán)卷積和線性卷積的結(jié)果相同，且計(jì)算循環(huán)卷積可以利用FFT算法[8-9]，因此可將g(n)和h(n)的線性卷積轉(zhuǎn)換為循環(huán)卷積。又因?yàn)橛?jì)算K點(diǎn)DTFT只需要輸出0，1，…，K-1共K個(gè)點(diǎn)的卷積值，即使后面N-1個(gè)點(diǎn)發(fā)生混疊也不會(huì)影響前面的值，因此可將循環(huán)卷積的周期縮短為L。

綜上所述，DTFT快速算法的基本思想如圖1所示，實(shí)現(xiàn)步驟如下：

步驟1確定g(n)和h(n)循環(huán)卷積的周期L，L≥K+N-1，且L為2的整數(shù)次冪。

步驟2將h(n)按式（7）轉(zhuǎn)化為一個(gè)長度為L點(diǎn)的新序列hL(n)：

步驟3根據(jù)已知參數(shù)按式（5）計(jì)算序列g(shù)(n)。

步驟4將g(n)補(bǔ)零成為長度為L的序列g(shù)1(n)，利用FFT算法計(jì)算g1(n)的L點(diǎn)離散傅里葉變換（Discrete Fourier Transform，DFT），記為G(l)。

步驟5利用FFT算法計(jì)算hL(n)的L點(diǎn)的DFT，記為H(l)。

步驟6將G(l)與H(l)相乘，得到g(n)和h(n)的循環(huán)卷積，記為Y(l)，并對其進(jìn)行L點(diǎn)快速傅里葉逆變換（Inverse Fast Fourier Transform，IFFT）得到序列y′(n)。取y′(n)在0≤n≤K-1范圍內(nèi)的值，即為g(n)和h(n)的線性卷積，記為y″(n)。

步驟7將y″(n)與表達(dá)式e-jπk2Δf/fs相乘，得到K點(diǎn)DTFT的值。

圖1 DTFT快速算法的基本思想

2.1.2 DTFT快速算法計(jì)算量分析

（1）計(jì)算g(n)：根據(jù)式（5），需要N次復(fù)數(shù)乘法。

（2）3次計(jì)算L點(diǎn)FFT（步驟4、5和6）：共需要1.5L lbL次復(fù)數(shù)乘法和3L lbL次復(fù)數(shù)加法。

（3）計(jì)算g(n)和h(n)的循環(huán)卷積Y(l)：需要L次復(fù)數(shù)乘法。

（4）計(jì)算y″(n)與e-jπk2Δf/fs的乘積：需要K次復(fù)數(shù)乘法。

故在DTFT快速算法中，計(jì)算K點(diǎn)DTFT的計(jì)算量約為S′DTFT×次復(fù)數(shù)乘法和S′DTFT+次復(fù)數(shù)加法：

直接K點(diǎn)DTFT需要復(fù)數(shù)乘法次數(shù)SDTFT×=KN，需要復(fù)數(shù)加法次數(shù)SDTFT+=(N-1)K。當(dāng)K，N取值變大時(shí)，DTFT快速算法的計(jì)算量將遠(yuǎn)小于常規(guī)DTFT算法。

2.2 加權(quán)融合頻譜矩陣降維處理

由文獻(xiàn)[4]可知，加權(quán)融合頻譜矩陣X′(fP)是由P×1階矩陣Xm(fP)經(jīng)M×P×Q(M≥2，P≥1，Q≥1)階加權(quán)因子e-jD加權(quán)融合生成，X′(fP)是一個(gè)P×Q階的矩陣，可用圖2所示的圖形表示。圖2是由P根橫線和Q根豎線編制的一個(gè)網(wǎng)格圖，X′(fP)中第(p，q)個(gè)元素X′q[fP(p)]用圖中第p條橫線和第q條豎線的交點(diǎn)表示，最優(yōu)加權(quán)融合頻譜Xq0'(fP)由第q0條豎線上所有交點(diǎn)表示。

圖2 加權(quán)融合頻譜矩陣降維原理示意圖

（1）當(dāng)P≠Q(mào)時(shí)，由DTFT頻率估計(jì)原理[10]可知，Xq0'(fP)可能位于X′(fP)任一列上。只能通過峰值搜索abs[Xq'(fP)]，根據(jù)其峰值元素所在列q0來確定Xq0'(fP)。因此必須計(jì)算出X′(fP)中P×Q個(gè)元素，計(jì)算量較大。

（2）當(dāng)P=Q時(shí)，①當(dāng)被測頻率f恰好等于序列fP中第p1(p1∈[1，P])個(gè)元素，即f=fP(p1)。其中，fP表示將fscope線性等分構(gòu)成的一個(gè)長度為P的序列。在無噪聲條件下，abs[Xq'(fP)]的峰值元素必定位于X′(fP)的 (p1，q1)處，此時(shí)p1=q1，即在X′(fP)的對角線元素上。②當(dāng)f等于fP中第p1個(gè)元素和第p2(p2∈[1，P])個(gè)元素之間的某個(gè)值時(shí)，假設(shè)f更接近fP(p1)，即f≈fP(p1)，在無噪聲情況下，abs[Xq'(fP)]的峰值元素也必定位于矩陣X′(fP)的(p1，q1)處，也在X′(fP)的對角線元素上。綜上所述，只需設(shè)置P=Q，abs[Xq'(fP)]的峰值元素必定會(huì)位于X′(fP)的對角線元素上。為確定最優(yōu)加權(quán)融合頻譜Xq0'(fP)只需計(jì)算由X′(fP)對角線上元素構(gòu)成的P×1階矩陣Xpq'(fP)即可。通過峰值搜索abs[Xpq'(fP)]即可確定Xq0'(fP)在X′(fP)的列數(shù)q0。按此方法可將原融合算法中計(jì)算P×Q階X′(fP)降維為計(jì)算P×1階的Xpq'(fP)，從而大幅度減少計(jì)算量。

Xpq'(fP)={X1'[fP(1)]，X2'[fP(2)]，…，XQ'[fP(Q)]}（10）

同時(shí)，對X′(fP)進(jìn)行降維處理還可提高算法的抗噪性。因?yàn)楫?dāng)P=Q時(shí)，雖然由理論分析Xq0'(fP)對應(yīng)的q0應(yīng)出現(xiàn)在X′(fP)的對角線元素上，但受噪聲干擾，若通過直接峰值搜索abs[Xq'(fP)]來確定的q0很可能不在X′(fP)的對角線上。但采用上述降維處理后，即使在噪聲干擾情況下，也能確保q0來源于X′(fP)對角線元素，因此從一定程度上提高了算法的抗噪性。

2.3 1/3主瓣相關(guān)性分析

在原融合算法中，通過計(jì)算同一信號(hào)的兩個(gè)不同類型頻譜——多段正弦信號(hào)的最優(yōu)加權(quán)融合頻譜Xq0'(fP)和其累加頻譜X(fP)的頻域相關(guān)譜r(fP)，來抑制虛假譜峰和噪聲干擾，提高信號(hào)參數(shù)估計(jì)精度。由圖3可知，Xq0'(fP)僅在其1/3主瓣對應(yīng)的頻段內(nèi)與其X(fP)對應(yīng)部分具有最大的相似度，因此沒有必要分析Xq0'(fP)和X(fP)在整個(gè)頻段內(nèi)的相關(guān)性，僅需計(jì)算圖3中灰色區(qū)域的1/3主瓣相關(guān)譜r13(fP)即可，從一定程度上消減原融合算法的計(jì)算量，其數(shù)學(xué)表達(dá)式如式（11）和式（12）所示：

圖31 /3主瓣相關(guān)性分析示意圖

表1 多段正弦頻譜融合法計(jì)算量

表2 多段正弦快速頻譜融合法計(jì)算量

3 算法計(jì)算量分析

將算法所需的復(fù)數(shù)乘法和復(fù)數(shù)加法次數(shù)按如下關(guān)系轉(zhuǎn)化為實(shí)數(shù)乘法和實(shí)數(shù)加法：一次復(fù)數(shù)乘法需要四次實(shí)數(shù)乘法和兩次實(shí)數(shù)加法，一次復(fù)數(shù)加法需兩次實(shí)數(shù)加法。設(shè)置M=4，P=150，Q=150，Nm=50，N′=200，Lm=256(m∈[1，M])，Pl1l2=P/5，采用DSP2812[11]為核心處理器，設(shè)置其主頻為40 MHz（其峰值可達(dá)150 MHz）。以計(jì)算一次實(shí)數(shù)乘法運(yùn)算需要四個(gè)機(jī)器周期，計(jì)算一次實(shí)數(shù)加法運(yùn)算需要一個(gè)機(jī)器周期為例對算法計(jì)算量進(jìn)行分析，分析結(jié)果如表3所示。原融合法的總耗時(shí)均約為22.868 ms，能夠滿足普通應(yīng)用所需。本文提出的快速融合算法的總耗時(shí)約為3.575 ms，較原融合算法的運(yùn)算速度提高了約84.37%。當(dāng)提高DSP的主頻時(shí)，上述算法的耗時(shí)將進(jìn)一步降低。

表3 算法計(jì)算量和耗時(shí)說明表

4 仿真實(shí)驗(yàn)

（1）信噪比變化條件下的仿真實(shí)驗(yàn)

為測試快速融合算法在不同信噪比條件下的性能，針對信號(hào)b1和信號(hào)b2，分別進(jìn)行了11組實(shí)驗(yàn)，每組包括1 000次Monte_Carlo實(shí)驗(yàn)，各組實(shí)驗(yàn)的信噪比設(shè)置如表5第一列所示，其余參數(shù)見表4。

表4 實(shí)驗(yàn)參數(shù)設(shè)置

表5 不同信噪比條件下的頻率估計(jì)均方根誤差kHz

由表5所示的仿真結(jié)果可知，快速融合算法在估計(jì)信號(hào)b2頻率時(shí)的RMSE低于估計(jì)信號(hào)b1的頻率。在信噪比為-11～5 dB的范圍內(nèi)，快速融合算法頻率估計(jì)的RMSE略高于原融合算法；在信噪比為-15～-13 dB的實(shí)驗(yàn)范圍內(nèi)，快速融合算法頻率估計(jì)的RMSE略低于原融合算法；這是因?yàn)榭焖偃诤纤惴ㄖ袑訖?quán)融合頻譜矩陣的降維處理能夠進(jìn)一步提高算法抗噪性，因此在信噪比較低時(shí)其性能略優(yōu)于原融合算法。

（2）單段信號(hào)長度變化條件下的仿真實(shí)驗(yàn)

為測試快速融合算法在單段信號(hào)長度變化條件下的性能，針對信號(hào)b1和信號(hào)b2，分別進(jìn)行了10組實(shí)驗(yàn)，每組包括1 000次Monte_Carlo實(shí)驗(yàn)，每組b1和b2中第m段信號(hào)的長度Nb1m和Nb2m設(shè)置如表6所示，其余參數(shù)參見表4。

表6 多段正弦信號(hào)中單段信號(hào)長度設(shè)置表

由表7所示的仿真結(jié)果可知，隨著單段信號(hào)長度在表6所示范圍內(nèi)變化?？焖偃诤纤惴ǖ念l率估計(jì)RMSE隨著單段信號(hào)長度的增加而遞減，略大于原融合算法頻率估計(jì)的RMSE，整體性能穩(wěn)定。

表7 不同單段信號(hào)長度條件下的頻率估計(jì)均方根誤差kHz

（3）被測頻率變化條件下的仿真實(shí)驗(yàn)

為測試快速融合算法在被測頻率f變化條件下的性能，進(jìn)行了10組實(shí)驗(yàn)，每組實(shí)驗(yàn)中f設(shè)置如表8第一列所示，每組包括1 000次Monte_Carlo實(shí)驗(yàn)，其余參數(shù)見表4。

表8 不同被測頻率f條件下的頻率估計(jì)均方根誤差kHz

由表8所示的仿真結(jié)果可知，隨著被測頻率f在7.5～12 MHz內(nèi)變化，在估計(jì)信號(hào)b1頻率時(shí)，快速融合算法和原融合算法的頻率估計(jì)平均RMSE分別為38.27 kHz，37.41 kHz；在估計(jì)信號(hào)b2頻率時(shí)，快速融合算法和原融合算法的頻率估計(jì)平均RMSE分別為37.35 kHz，36.22 kHz。上述仿真結(jié)果表明快速融合算法的性能較為接近原融合算法，較穩(wěn)定；在滿足奈特斯特采樣定理的前提下，能夠適用于不同的被測頻率對象，具有良好的普適性。

（4）采樣頻率變化條件下的仿真實(shí)驗(yàn)

為測試快速融合算法在采樣頻率fs變化條件下的性能，進(jìn)行了10組實(shí)驗(yàn)，每組實(shí)驗(yàn)中fs設(shè)置如表9第一列所示，每組包括1 000次Monte_Carlo實(shí)驗(yàn)，其余參數(shù)設(shè)置見表4。

表9 不同采樣頻率fs條件下的頻率估計(jì)均方根誤差kHz

由表9所示的仿真結(jié)果可知，隨著采樣頻率fs在38～42.5 MHz內(nèi)變化，在估計(jì)信號(hào)b1頻率時(shí)，快速融合算法和原融合算法的頻率估計(jì)平均RMSE分別為38.44 kHz，37.54 kHz；在估計(jì)信號(hào)b2頻率時(shí)，快速融合算法和原融合算法的頻率估計(jì)平均RMSE分別為36.99 kHz，36.04 kHz。上述仿真結(jié)果表明快速融合算法的性能較為接近原融合算法，在滿足奈特斯特采樣定理的前提下，其性能基本不依賴于采樣頻率。

5 總結(jié)

為進(jìn)一步提高多段正弦信號(hào)頻譜融合算法的實(shí)時(shí)性，滿足雷達(dá)、聲納、電子對抗等應(yīng)用需求，本文通過設(shè)計(jì)DTFT快速算法、降維處理加權(quán)融合頻譜矩陣和1/3主瓣相關(guān)性分析處理等措施來降低算法計(jì)算量，提出頻率估計(jì)的多段正弦信號(hào)快速頻譜融合算法。計(jì)算量對比和仿真實(shí)驗(yàn)表明，多段正弦信號(hào)快速頻譜融合算法在精度損失極小的前提下，能夠大幅降低計(jì)算量；在信噪比極低的情況下(SNR≤-13 dB)，其性能略優(yōu)于原融合算法。

[1]Xiao Wei，Tu Yaqing，Liu Lingbing.Parameters estimation of LFM signal based on fusion of signals with the same length and known frequency-difference[C]//IEEE Proceedings of 8th World Congress on Control and Automation（WCICA’2010），Jinan，China，2010：6776-6781.

[2]Wu Guibo，Zhang Lijun，Lu Hui，et al.A real-time andhigh-precision algorithm for frequency estimation byfusion multi-segment signals[J].Procedia Engineering，2011，15：2313-2317.

[3]劉良兵，涂亞慶.基于多段分頻等長信號(hào)融合的頻率估計(jì)方法[J].計(jì)算機(jī)工程與應(yīng)用，2008，44（11）：170-175.

[4]肖瑋，涂亞慶，劉良兵，等.頻率估計(jì)的多段同頻正弦信號(hào)頻譜相關(guān)算法[J].電子與信息學(xué)報(bào)，2012，34（3）：564-570.

[5]Liu Liangbing，Tu Yaqing，Xu Baosong.A division ratiovariable delay method for VCO FM nonlinearity correction[C]//IEEE Proceedings of 7th World Congress on Control and Automation（WCICA’2008），Chongqing，China，2008：2692-2695.

[6]閻振華，黃建國，韓晶.改進(jìn)的低信噪比短序列快速頻率估計(jì)算法[J].系統(tǒng)工程與電子技術(shù)，2009，31（8）：1990-1992.

[7]張海濤，涂亞慶.計(jì)及負(fù)頻率影響的科里奧利質(zhì)量流量計(jì)信號(hào)處理方法[J].儀器儀表學(xué)報(bào)，2007，28（3）：539-544.

[8]Prioakis J G，Manolakis D G.Digital signal processing：principle，algorithms，and application[M].New Jersey：Prentice Hall，2006.

[9]胡廣書.數(shù)字信號(hào)處理——理論、算法與實(shí)現(xiàn)[M].北京：清華大學(xué)出版社，2003.

[10]Ayhan B，Trussell H J，Chow M Y，et al.On the use of a lower sampling rate for broken rotor bar detection with DTFT and AR-based spectrum methods[J].IEEE Transactions on Industrial Electronics，2008，55（3）：1421-1434.

[11]Li Fengpan，Lu Wenhua.Turbine governor system frequency measurement based on DSP[C]//2011 2nd International Conference on Control，Instrumentation and Automation（ICCIA），2011：459-462.

XIAO Wei1,TU Yaqing1,SHEN Yanlin1,SU Dan1,ZHANG Lei2

1.Department of Logistical Information Engineering,Logistical Engineering University,Chongqing 401311,China
2.Zhuozhou Comprehensive Storehouse,Baoding,Hebei 611730,China

The spectra fusion method for multi-sections sinusoids（called“the fore-fusion method”for short）is an effective way of estimation frequency for the sinusoids in low Signal-to-Noise Ratio（SNR）,which has an important theoretical significance and practical value.In order to meet the high real-time demand in some fields such as radar,sonar and electronic countermeasures,a fast spectra fusion algorithm for multi-section sinusoids is put forward.This proposed algorithm can reduce the calculation and improve the real-time characteristic by the following techniques：design a fast DTFT algorithm,reduce dimensions of the weighted fusion spectrum matrix,and analyse the correlation of the 1/3 main-lodes of the optimization weighted-accumulation spectrum and the accumulation spectrum.The principles of the above techniques are expatiated.Calculation analyses and simulations demonstrate that the proposed algorithm can reduce most calculation of the fore-fusion method with lower little precision,and it works better in very low SNR(SNR≤-13 dB).

fast spectra fusion;multi-sections sinusoids;fast Discrete Time Fourier Transform（DTFT）algorithm;dimensions reduction;correlation of 1/3 main-lodes

A

TN957

10.3778/j.issn.1002-8331.1207-0435

XIAO Wei,TU Yaqing,SHEN Yanlin,et al.Fast frequency estimation algorithm based on spectra fusion of multi-section sinusoids.Computer Engineering and Applications,2014,50（6）：186-191.

國家自然科學(xué)基金（No.61271449，No.61302175）；重慶市自然科學(xué)基金（No.cstc2012jjA0877，No.cstc2011BA2015）；后勤工程學(xué)院青年科研基金項(xiàng)目資助課題。

肖瑋（1982—），女，博士，研究領(lǐng)域?yàn)樾盘?hào)處理，嵌入式系統(tǒng)；涂亞慶（1963—），男，博士，教授，博士生導(dǎo)師，研究領(lǐng)域?yàn)樾盘?hào)處理；沈艷林（1987—），男，碩士研究生，研究領(lǐng)域?yàn)樾盘?hào)處理；蘇丹（1983—），女，博士研究生，研究領(lǐng)域?yàn)樾盘?hào)處理；張磊（1980—），男，博士，研究領(lǐng)域?yàn)檐娛挛锪髋c信號(hào)處理。E-mail：WZWRY@163.com

2012-07-29

2012-11-23

1002-8331（2014）06-0186-06

CNKI網(wǎng)絡(luò)優(yōu)先出版：2012-12-13，http://www.cnki.net/kcms/detail/11.2127.TP.20121213.1629.003.html