張海燕，李海云

首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，北京，100069

核磁共振圖像腦組織自動提取方法

張海燕，李海云

首都醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，北京，100069

核磁共振圖像的腦組織提取是神經(jīng)圖像處理研究中的一個重要步驟。將傳統(tǒng)的幾何活動輪廓模型與二值水平集函數(shù)相結(jié)合，提出了一種新型的二值水平集活動輪廓模型，并基于該模型提出了一種能夠自動、準(zhǔn)確實現(xiàn)MRI腦組織提取的方法。該方法在腦組織內(nèi)部自動設(shè)定最優(yōu)初始輪廓曲線，將該演化曲線隱含地表示成一個高維函數(shù)的零水平集，零水平集在基于區(qū)域的圖像力驅(qū)動下不斷演化并達(dá)到待分割腦部圖像的邊緣。將基于該方法的腦組織提取結(jié)果與作為金標(biāo)準(zhǔn)的專家手動分割結(jié)果和其他流行算法相比較，結(jié)果表明提出的腦組織提取方法能夠自動、準(zhǔn)確和快速地提取MRI腦組織，是一種魯棒性較好的MRI腦組織提取方法。

腦組織提??；灰度直方圖；活動輪廓模型；水平集；質(zhì)量評估

1 引言

核磁共振圖像（M agnetic Resonance Imaging，MRI）腦組織提取是將核磁共振序列腦圖像中的腦部組織與非腦組織分離，去除腦外組織，也稱為頭骨剝離。腦組織的準(zhǔn)確提取是神經(jīng)圖像處理中一個重要的步驟，也是其他圖像處理算法例如配準(zhǔn)、腦組織分類或者灰度不均勻性校正等的預(yù)處理步驟[1-3]，例如，只有準(zhǔn)確地去除非腦組織才能夠?qū)崿F(xiàn)腦內(nèi)組織如腦白質(zhì)、腦灰質(zhì)和腦脊液的準(zhǔn)確分割[4]。腦組織的準(zhǔn)確提取不僅是腦組織體積測量以及實現(xiàn)三維重建的關(guān)鍵技術(shù)，在神經(jīng)圖像處理和分析中也是最耗時的預(yù)處理步驟之一，因此國際上提出了許多的腦組織提取算法（Brain Extraction A lgorithm s，BEAs），目前許多研究致力于開發(fā)準(zhǔn)確、自動的腦組織提取算法，這些算法在不同程度上提高了腦組織提取的精度或速度[5-6]，但是各種算法的分割質(zhì)量也大不相同，從而會影響后續(xù)的圖像分析[7]。Fennema-Notestine等人[8]比較了最常用的幾種頭骨剝離算法的性能，認(rèn)為每種算法都有自己的優(yōu)點和缺點，但是還沒有一種算法能夠適用于大規(guī)模圖像的處理與分析。

目前MRI腦組織提取方法大體可分為三類：基于圖像灰度的方法、基于形態(tài)學(xué)的方法和基于活動輪廓模型的方法?；诨顒虞喞Ｐ偷姆椒ㄊ悄壳皥D像分割中應(yīng)用最廣泛的一種方法。Smith等人[9]提出了一個用于頭骨剝離的Brain Extraction Tool（BET）算法，該算法將一組基于形態(tài)學(xué)和基于圖像特征的力應(yīng)用于演化曲面的切向和法向，驅(qū)動演化曲線到達(dá)待分割圖像邊緣。Zhuang等人[10]提出了一個用于頭骨分離的基于模型的水平集方法（Model-based Level Set，M LS），該模型將演化曲線作為高維函數(shù)的零水平集，利用基于腦皮層灰度特征的圖像力和曲率力控制曲線的演化。形變模型的主要優(yōu)點是能夠直接產(chǎn)生閉合的參數(shù)曲線或曲面，并對噪聲和偽邊界有較強的魯棒性。與其他兩種方法相比，基于形變模型的方法更容易實現(xiàn)自動的頭骨剝離，但對初始輪廓線的位置較敏感、計算量大造成分割時間較長。

隨著MR圖像數(shù)量的增長及具有臨床應(yīng)用意義的實時圖像處理的需要，開發(fā)全自動、準(zhǔn)確和快速的腦組織提取算法變得尤為重要。本文基于二值水平集活動輪廓模型提出了一種新的MRI腦組織提取算法，使用該算法對臨床收集的MR真實腦圖像和網(wǎng)絡(luò)上下載的的MR腦圖像數(shù)據(jù)自動分割腦組織與非腦組織邊界，并將分割結(jié)果與專家手動分割結(jié)果相比較，實驗結(jié)果表明該算法能自動、準(zhǔn)確地進(jìn)行腦組織提取而且有很好的魯棒性。

2 基于活動輪廓模型的圖像分割的研究

2.1 幾何活動輪廓模型

傳統(tǒng)幾何活動輪廓模型（Geometric Active Contour Model，GACM）是由Caselles等人[11]與Malladi等人[12]分別獨立地提出的基于曲線演化和水平集方法的活動輪廓模型，該模型利用水平集函數(shù)來隱含地表示模型輪廓線，從而成功地解決了輪廓線拓?fù)渥兓膯栴}。傳統(tǒng)的幾何活動輪廓模型的曲線演化方程可表示成以下形式：

傳統(tǒng)幾何活動輪廓模型是通過檢測圖像邊緣來提取目標(biāo)邊界的一種方法，對有較好對比度的圖像分割效果較好。

2.2 二值水平集活動輪廓模型

傳統(tǒng)水平集活動輪廓模型雖然具有自動處理輪廓線拓?fù)渥兓膬?yōu)點，但為了保證水平集方法在數(shù)值實現(xiàn)過程中的穩(wěn)定性和精確性，水平集函數(shù)在演化過程中需要反復(fù)地被重新初始化為其零水平集的符號距離函數(shù)，而這個初始化過程是個非常耗時的運算過程。Zhang等人[13]提出了一種二值水平集活動輪廓模型。在該活動輪廓模型中，水平集函數(shù)僅取-1和1這兩個數(shù)值。在二值水平集函數(shù)更新過程中，將該二值水平集函數(shù)與一個二值圖像相互進(jìn)行轉(zhuǎn)化，并使用形態(tài)學(xué)中的腐蝕與膨脹算子對該二值圖像進(jìn)行操作，使其模擬曲線的演化過程。與傳統(tǒng)水平集活動輪廓模型相比，二值水平集活動輪廓模型在運行效率上得到了極大的提高，同時保持了自動處理輪廓線的拓?fù)渥兓哪芰?，但卻喪失了曲線演化的漸進(jìn)性，從而使得分割效果類似于簡單的閾值分割法，不具有針對性。

2.3 一種新型的二值水平集活動輪廓模型

由于MR圖像上腦組織與非腦組織的邊緣較弱，腦組織之間邊界模糊以及圖像噪聲的影響，導(dǎo)致MR圖像的分割問題非常復(fù)雜和困難?；趥鹘y(tǒng)幾何輪廓模型的分割方法使用圖像的邊界梯度信息，邊緣定位較精確，性能上優(yōu)于傳統(tǒng)的基于邊緣檢測的分割方法，但是也存在一定的局限性，即對噪聲和初始條件很敏感，檢測腦組織和非腦組織的邊界時輪廓曲線容易發(fā)生邊緣泄露或者過分割，而且在輪廓線的初始化問題上，初始輪廓線必須被完全地設(shè)置在目標(biāo)邊界的內(nèi)部或者外部，否則無法得到正確的分割結(jié)果。本文在傳統(tǒng)幾何活動輪廓模型的框架下，將圖像的邊界信息和區(qū)域統(tǒng)計信息相結(jié)合，提出了一種新型的二值水平集活動輪廓模型。首先構(gòu)造具有圖像統(tǒng)計特征信息的SPF函數(shù)如下：

其中Gσ是標(biāo)準(zhǔn)偏差為σ的高斯核函數(shù)，*代表卷積算子。f1和f2為兩個平滑常數(shù)，定義如下：

其中是正則化的Heaviside函數(shù)，Gσ是方差為σ的高斯核函數(shù)。

方程（3）中的SPF函數(shù)也可以稱為區(qū)域函數(shù)[14]，它的取值范圍為[-1，1]。區(qū)域函數(shù)可以利用區(qū)域統(tǒng)計信息調(diào)節(jié)壓力的符號，使得當(dāng)演化曲線位于感興趣區(qū)域內(nèi)部時演化曲線膨脹或位于外部時收縮，因此SPF函數(shù)可以解決所謂的由弱邊緣引起的邊緣泄露問題。用方程（3）中的SPF函數(shù)代替方程（2）的邊緣檢測函數(shù)g，得到本文提出的用于腦組織提取的水平集演化方程如下：

3 基于二值水平集活動輪廓模型的MR圖像腦組織提取方法

基于以上提到的新型二值水平集活動輪廓模型，本文提出的MR圖像腦組織提取方法主要包括三個步驟：首先根據(jù)圖像灰度直方圖估計圖像灰度參數(shù)以得到頭部的二值圖像；然后在腦區(qū)內(nèi)部自動設(shè)定自適應(yīng)的初始輪廓曲線，并將該輪廓曲線隱含地表示成高維水平集函數(shù)的零水平集，零水平集在圖像力的驅(qū)動下不斷演化并達(dá)到待分割腦部圖像的邊緣。最后介紹了該方法中的人腦MR圖像分割的簡化原則和兩頭分割方案。

3.1 圖像灰度參數(shù)估計及頭部的二值化圖像

采用Sm ith[9]提出的方法來估計灰度的有效范圍。通過分析人腦MR圖像的直方圖特征，選取和設(shè)定高低兩個閾值t1和t2構(gòu)成灰度的有效范圍[t1，t2]，其中較低的閾值t1用來去除代表背景，顱骨，空氣和其他較低灰度級噪點的立體像素；較高的t2閾值用來消去類似脂肪和血管系統(tǒng)的較亮的非腦組織。然后根據(jù)學(xué)習(xí)經(jīng)驗，確定腦部區(qū)域和非腦區(qū)域的分割閾值tm：

其中局部閾值Th為一常數(shù)，通常該常數(shù)取值由圖像數(shù)據(jù)的學(xué)習(xí)經(jīng)驗確定，且與圖像的掃描方向有關(guān)，即橫軸位，冠狀位和矢狀位三個方向的取值不同。分割閾值tm確定后，以該閾值將腦部組織和非腦組織進(jìn)行粗略分割，得到頭部二值化圖像。

3.2 自適應(yīng)的水平集初始輪廓曲線的自動設(shè)定

水平集初始輪廓曲線的自動設(shè)定不僅是保證腦組織提取完全自動化的首要步驟，同時由于分割結(jié)果的準(zhǔn)確性對初始輪廓曲線極為敏感，所以初始輪廓線應(yīng)盡可能地靠近實際的輪廓才能得到最好的分割結(jié)果。另外由于一些非腦組織（例如皮膚和肌肉組織）具有和腦組織相似的灰度特性，所以如果初始輪廓曲線包含這些非腦組織分割時就會產(chǎn)生錯誤?，F(xiàn)存的初始化方法通常存在準(zhǔn)確性不高或者運算效率不高的缺點。Zhuang等人[10]將腦部區(qū)域的中心作為初始化水平集的圓心，取一個較小的半徑作為初始化水平集的半徑，以保證初始化水平集完全在腦區(qū)內(nèi)部。這種水平集的初始化方法雖然能提高分割結(jié)果的正確性，但是由于初始水平集的半徑很小，降低了運算效率。本文提出了一種自適應(yīng)的初始化水平集的自動設(shè)定方法，取二值化頭部圖像腦區(qū)的中心為初始化水平集的圓心，但是在保證初始化水平集在腦區(qū)內(nèi)部的前提下半徑盡可能的大，從而提高了運算效率，如圖1所示。初始化方法如下：

（1）對原始圖像用閾值tm處理后得到腦部二值化圖像。

（2）對二值化后的圖像計算其中腦部區(qū)域的最左點和最右點，得到最左列和最右列的直線距離d。

（3）通過觀察頭顱MR圖像特點，發(fā)現(xiàn)頭顱大小在MR圖像的三個掃描方向上有所不同，因此根據(jù)學(xué)習(xí)經(jīng)驗，對于橫軸位的MR圖像，可以將顱腦形狀近似為邊長等于d的正方形；而對于冠狀位和矢狀位的MR圖像，顱腦形狀都可近似為矩形，矩形的寬度都等于d，但是寬與高之比近似為5∶4。

（4）基于步驟3得到顱腦的近似形狀，初始化零水平集曲線為一個圓，圓心位于正方形或矩形的中心，半徑分別等于正方形邊長的三分之一（相對于橫軸位MR圖像）和矩形寬度的五分之一（相對于冠狀位和矢狀位MR圖像）。這種初始輪廓曲線的自動設(shè)定過程稱之為自適應(yīng)的初始輪廓曲線自動設(shè)定方法，即根據(jù)輸入的顱腦MR圖像掃描方向或個體差異可自己調(diào)整初始輪廓線的位置及大小，從而保證初始水平集在腦區(qū)內(nèi)部并盡可能的靠近圖像待分割邊緣。

3.3 二值水平集函數(shù)設(shè)定及分割過程中簡化原則和兩頭分割方案

為簡化水平集函數(shù)的初始化過程，使用二值水平集函數(shù)取代傳統(tǒng)水平集函數(shù)。初始化水平集函數(shù)如下：

通過對MR腦圖像進(jìn)行研究，發(fā)現(xiàn)同一序列的前幾層和后幾層圖像上包含的腦組織區(qū)域很小甚至為零，更談不上分割腦組織的邊緣，所以前面和后面的幾層可以不予理會。根據(jù)學(xué)習(xí)經(jīng)驗提出以下簡化原則：算法分割從每一圖像序列的第十分之一層開始，到第十分之九層結(jié)束。又因為同一序列顱腦MR圖像的相鄰層具有較強的相干性，所以采取如下的初始輪廓曲線設(shè)定方案：（1）從開始分割的層數(shù)到序列的中間層，由于腦區(qū)面積呈現(xiàn)增大的趨勢，所以當(dāng)前層的最終演化輪廓線可以作為后面相鄰層的初始輪廓線，這樣既保證初始輪廓曲線在腦區(qū)內(nèi)部又使初始輪廓曲線盡可能地靠近待分割邊緣；（2）從序列的中間層到分割結(jié)束的層數(shù)，由于腦區(qū)面積呈現(xiàn)逐漸減小的趨勢，所以為了保證初始輪廓曲線在腦區(qū)內(nèi)部，采取反向分割順序，即先從最后一層開始分割，將其作為前面相鄰層的初始輪廓線，直到中間層結(jié)束。這種方案也可以稱為兩頭分割方案，即分割時分別從前、后兩頭向中間層過渡。實驗結(jié)果表明以上介紹的分割簡化原則和兩頭分割方案，不僅可以完全自動地在腦區(qū)內(nèi)部設(shè)定與待分割邊緣靠近的初始輪廓線，而且大大提高了算法效率和分割結(jié)果的準(zhǔn)確性。

4 實驗結(jié)果與質(zhì)量評估

本文采用M atlab 7.0作為軟件實驗平臺，算法測試的硬件配置為CPU：Intel雙核T6500 2.1 GHz，內(nèi)存2 GB。測試數(shù)據(jù)為臨床收集的10組正常人的T1 MR腦圖像，每位受試者的體數(shù)據(jù)共176層，空間分辨率為256×256，層厚為1 mm，無間隔。圖2為應(yīng)用本文算法對10例正常人腦MR圖像的部分腦組織提取結(jié)果。

圖1 自適應(yīng)的初始輪廓曲線自動設(shè)定方法

圖2 MR TIW I真實圖像橫軸位、冠狀位和矢狀位的腦組織提取結(jié)果

為了對本文算法的分割效果進(jìn)行質(zhì)量評估，從網(wǎng)絡(luò)腦圖像標(biāo)準(zhǔn)分割庫（Internet Brain Segmentation Repository，IBSR）[15]下載10組正常人體大腦MR圖像（包括待分割圖像數(shù)據(jù)和對應(yīng)的專家手動分割結(jié)果）。用提出的方法對這10組正常人體冠狀位MR腦圖像進(jìn)行腦組織提取，分割結(jié)果如圖3所示。將算法自動分割結(jié)果與專家手動分割結(jié)果進(jìn)行比較，采用靈敏度（sensitivity）、特異性（specificity）和FP_Rate系數(shù)來衡量算法自動分割結(jié)果與IBSR提供的專家手動分割結(jié)果的一致性：

靈敏度和特異性分別是算法識別出的腦組織和非腦組織的正確分類的像素的百分比：

FP_Rate系數(shù)是算法自動分割結(jié)果中把實際上是非腦組織的體素誤判為腦組織類的體素個數(shù)與專家手動分割結(jié)果中分類為腦組織體素個數(shù)的比值。因此FP_Rate系數(shù)的值越低，表明算法的分割結(jié)果越準(zhǔn)確。

圖3 本文算法的分割結(jié)果與專家手動分割結(jié)果

其中TP表示判斷正確的腦部組織的體素；TN表示判斷正確的非腦組織區(qū)域；FN表示分割結(jié)果中把實際上是腦組織的體素誤判為非腦組織類的體素個數(shù)；FP表示分割結(jié)果中把實際上是非腦組織的體素誤判為腦組織類的體素個數(shù)。靈敏度和特異性的取值在0～1之間，通常來講，靈敏度的值越大，分割結(jié)果越準(zhǔn)確。但是有一種特例，那就是如果分割算法比較保守（臨床應(yīng)用時通常要采取保守措施），即寧肯將許多非腦組織誤判為腦組織而盡量避免將任何腦組織分類為非腦組織，那么此時FN就等于0，所以這種情況下靈敏度總是等于1。所以一個分割算法不能只是單一的使用靈敏度系數(shù)，必須綜合其他系數(shù)進(jìn)行評估。

將本文提出的算法分割結(jié)果與專家手動分割結(jié)果比較，以專家手動分割結(jié)果為模板計算算法的靈敏度、特異性和FP_Rate系數(shù)，結(jié)果如表1所示。低靈敏度表示分割結(jié)果遺漏一些腦部區(qū)域，低特異性表示頭骨剝離后結(jié)果中腦部區(qū)域錯誤地包含非腦區(qū)域，而FP_Rate系數(shù)的值越低表明算法的誤判率越低。從表1可見，綜合靈敏度、特異性和FP_Rate系數(shù)三個參數(shù)，本文提出的方法與傳統(tǒng)的腦組織提取方法BET，BSE相比具有明顯的優(yōu)越性，分割精度明顯提高，因為BET算法在分割過程中總是采取保守措施，所以它的靈敏度最高但也導(dǎo)致特異性和FP_Rate系數(shù)性能最差。BET算法如此設(shè)計的初衷可能是因為從臨床應(yīng)用的角度出發(fā)，避免誤判任何的腦組織顯然要比去掉所有的非腦組織更重要。與M LS相比，本文的算法在分割精度上性能稍差，但FP_Rate和Specificity兩項系數(shù)性能略優(yōu)于M LS；分割效率方面，由于本文提出的初始化水平集自動設(shè)定方法在很大程度上提高了運算效率，因此與其他三種算法相比，分割速度具有明顯的優(yōu)越性，因此提出的腦組織提取方法可有效地實現(xiàn)MR圖像腦組織的自動和準(zhǔn)確提取。

表1 各種算法性能比較

5 結(jié)束語

本文提出了一種新型的二值水平集活動輪廓模型用于腦組織提取。與傳統(tǒng)的腦組織提取方法相比，該方法具有以下優(yōu)點：初始化方法使用二值水平集函數(shù)代替?zhèn)鹘y(tǒng)水平集函數(shù)，避免了傳統(tǒng)算法中費時的重新初始化步驟，而且自適應(yīng)的水平集輪廓曲線能夠盡可能地靠近待分割的邊界，從而大大提高了算法效率；第二，提出的活動模型利用圖像的統(tǒng)計信息作為約束條件，對具有弱邊界的腦組織能夠獲得較好的分割效果；第三，可以在腦區(qū)內(nèi)部自動設(shè)定初始水平集，實現(xiàn)了腦組織提取的全自動化；最重要的是，該方法非常簡單，對原始MR圖像直接分割，無需任何預(yù)處理步驟，方便使用。因此本文提出的腦組織提取方法可有效地實現(xiàn)腦組織和非腦組織的自動和準(zhǔn)確分割。然而，本文提出的模型對一些含有病變組織的MR腦圖像，分割精度降低，因此可繼續(xù)深入研究，使算法能夠適用于大規(guī)模MR腦圖像的處理與分析。

[1]Woods R P，Dapretto M，Sicotte N L，et al.Creation and use of a talairach-compatible atlas for accurate，automated，nonlinear intersubject registration，and analysis of functional imaging data[J].Human Brain Mapping，1999，8：73-79.

[2]Shattuck D W，Sandor-leahy S R，Schaper K A，et al.Magnetic resonance image tissue classification using a partial volume model[J].Neuro Image，2001，13：856-876.

[3]Strother S，La Conte S，Kai Hansen L，et al.Optimizing the fMRI dataprocessing pipeline using prediction and reproducibility performance metrics：I.A preliminary group analysis[J].Neuro Image，2004，23：196-207.

[4]Zhang Y Y，Brady M，Smith S.Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm[J].IEEE Trans on Med Imaging，2001，20：45-57.

[5]M arroquin J，Vemuri B.An accurate and efficient Bayesian method for automatic segmentation of brain MRI[J].IEEE Trans on M ed Imaging，2002，21：934-945.

[6]Liu J X，Chen Y S，Chen L F.Accurate and robust extraction of brain regions using a deformable model based on radial basis functions[J].Journal of Neuroscience Methods，2009，183：255-266.

[7]Boesen K，Rehm K，Schaper K，et al.Quantitative comparison of four brain extraction algorithms[J].Neuro Image，2004，22（3）：1255-1261.

[8]Fennema-notestine C，Ozyurt I B，Clark C P，et al.Quantitative evaluation of automated skull-stripping method applied to contemporary and legacy images，effects of diagnosis，bias correction，and slice location[J].Hum Brain Mapp，2006，27（2）：99-113.

[9]Sm ith S M.Fast robust automated brain extraction[J]. Hum Brain Mapp，2002，17：143-155.

[10]Zhuang A H，Valentino D J，Toga A W.Skull-stripping magnetic resonance brain images using a model-based level set[J].Neuro Image，2006，32：79-92.

[11]Caselles V，Catte F，Coll T，et al.A geometric model for active contours in image processing[J].Numerische Mathematik，1993，66（1）：1-31.

[12]Malladi R，Sethian J，Vemuri B.Shape modeling with front propagation：a level set approach[J].IEEE Transactions on Pattern Analysis and Machine Intelligence，1995，17（2）：158-175.

[13]Zhang S.A fast level set implementation method for image segmentation and object tracking[C]//Proc SPIE，Applications of Digital Image Processing XXIX，63121R， doi：10.1117/12.682228，2006-08-24.

[14]Zhang K H，Zhang L，Song H H.Active contours with selective local or global segmentation：a new formulation and level set method[J].Image and Vision Computing，2010，28：668-676.

[15]IBSR was developed by the center for Morphometric Analysis at Massachusetts General Hospital.Internet Brain Segmentation Repository[DB/OL].[2012-05-12].http：//www. cma.mgh.harvard.edu/ibsr.

ZHANG Haiyan,LI Haiyun

College of Biomedical Engineering,Capital Medical University,Beijing 100069,China

The segmentation of brain region from non-brain region in magnetic resonance images,also referred to brain tissue extraction,is an important and difficult task due to the convoluted brain surface and weak boundaries between brain and non-brain tissue.This paper proposes an accurate and automatic method of brain tissue extraction based on a novel binary level set active contour model.An initial contour,which is automatically set inside the brain,is driven by a region-based image force until it well fit the brain and non-brain boundaries.Quantitative evaluation of the new method is performed by comparing the result of the new method to those obtained using manual segmentation in 10 sets normal adult MR brain images.Experimental results demonstrate that the proposed method can provide accurate and automatic brain tissue extraction results.

brain tissue extraction;intensity histogram;active contour model;level set;quality evaluation

A

TP391

10.3778/j.issn.1002-8331.1209-0273

ZHANG Haiyan,LI Haiyun.Automated MRI brain tissue extraction method.Computer Engineering and Applications,2014,50（16）：168-172.

北京市自然科學(xué)基金（No.4122018）。

張海燕（1977—），女，博士，副教授，研究方向為醫(yī)學(xué)圖像處理，生物醫(yī)學(xué)工程。E-mail：zhangxyz@ccmu.edu.cn

2012-09-24

2013-01-16

1002-8331（2014）16-168-05

CNKI網(wǎng)絡(luò)優(yōu)先出版：2013-01-25,http://www.cnki.net/kcms/detail/11.2127.TP.20130125.1658.002.htm l