李 恪,姚崇斌,徐紅新,謝寶蓉,尚吉揚(yáng)

(上海航天電子技術(shù)研究所,上海 201109)

?

一種改進(jìn)的星上高光譜異常檢測(cè)算法研究

李 恪,姚崇斌,徐紅新,謝寶蓉,尚吉揚(yáng)

(上海航天電子技術(shù)研究所,上海 201109)

為解決傳統(tǒng)RX異常檢測(cè)算法導(dǎo)致的星上高光譜圖像異常檢測(cè)準(zhǔn)確率低、成本高等問題，對(duì)一種改進(jìn)的星上高光譜異常檢測(cè)算法進(jìn)行了研究。在傳統(tǒng)RX異常檢測(cè)算法的基礎(chǔ)上，用噪聲自適應(yīng)主成分變換對(duì)原始圖像進(jìn)行變換，選擇變換后大于設(shè)定的合適信噪比閾值的數(shù)據(jù)；用局部正交子空間投影(LOSP)算法將相應(yīng)數(shù)據(jù)投影到正交子空間獲得圖像殘差數(shù)據(jù)，通過抑制背景等強(qiáng)干擾信息而突出小概率的異常目標(biāo)信息；用空域?yàn)V波方法提取殘差數(shù)據(jù)的特征，將大部分圖像信息以某種標(biāo)準(zhǔn)集中于少數(shù)的主成分；對(duì)獲得的波段子集用中值濾波器濾波，消除噪聲干擾；用RX異常檢測(cè)算子對(duì)濾波后的波段子集進(jìn)行異常目標(biāo)檢測(cè)。試驗(yàn)結(jié)果證明了算法的有效性，算法通過數(shù)據(jù)量壓縮，降低維度、抑制干擾信息，減少了異常檢測(cè)處理的數(shù)據(jù)量以提高異常檢測(cè)效率和精度。

星載探測(cè)； 高光譜； 異常檢測(cè)； 局部正交子空間投影； RX異常檢測(cè)算法； 噪聲自適應(yīng)主成分變換； 信噪比； 閾值

0 引言

星載高光譜探測(cè)技術(shù)是一種新型的具“圖譜合一”特性的遙感技術(shù)，其連續(xù)的光譜曲線可更好地表達(dá)物質(zhì)間的細(xì)微差異，在物質(zhì)的分類、解混和目標(biāo)探測(cè)等領(lǐng)域獲得了廣泛應(yīng)用[1-3]。近年來，隨著衛(wèi)星遙感技術(shù)的發(fā)展，高光譜圖像的應(yīng)用越來越廣泛。由于高光譜遙感圖像包含豐富的空間、輻射和光譜信息三重信息，其圖像數(shù)據(jù)是一個(gè)融合了圖像二維空間信息與光譜維信息的三維立方體，具有較高的數(shù)據(jù)維、冗余性及較強(qiáng)的相關(guān)性[4-6]。與傳統(tǒng)的異常檢測(cè)手段相比，高光譜圖像異常目標(biāo)檢測(cè)的步驟更多，內(nèi)容更復(fù)雜。隨著圖像光譜分辨率的不斷提高，高光譜遙感圖像的數(shù)據(jù)量也日益增加，不僅給星上圖像處理增加了難度，而且導(dǎo)致異常檢測(cè)準(zhǔn)確率較低[7-9]。因此，改進(jìn)異常檢測(cè)算法對(duì)解決高光譜圖像的背景噪聲干擾、數(shù)據(jù)間的冗余性和高數(shù)據(jù)維等問題有重要的現(xiàn)實(shí)意義。

高光譜異常目標(biāo)檢測(cè)是在無需先驗(yàn)信息的情況下，從信息量的角度對(duì)圖像中的“小目標(biāo)”的“盲”檢測(cè)，此處的“小目標(biāo)”在處理過程中表現(xiàn)為一種數(shù)據(jù)異常。RX 檢測(cè)算法是一種經(jīng)典的異常檢測(cè)算法，由REED，YU兩人提出[10]。該算法是在局部正態(tài)模型的假設(shè)下，用最大似然比檢測(cè)得到的恒虛警檢測(cè)算法。由于無需先驗(yàn)信息而只需計(jì)算數(shù)據(jù)的協(xié)方差矩陣，其實(shí)用性較強(qiáng)。RX算法依據(jù)恒虛警檢測(cè)準(zhǔn)則分離出目標(biāo)，優(yōu)點(diǎn)是算法原理簡(jiǎn)單和計(jì)算復(fù)雜度低，適于背景分布較簡(jiǎn)單的目標(biāo)檢測(cè)，該法是基于假設(shè)數(shù)據(jù)空間白化的，未利用高光譜圖像波段間的非線性特性，對(duì)全局范圍內(nèi)的異常及相對(duì)背景差異不明顯的弱異常的檢測(cè)效果受到很大限制。這是因?yàn)椋菏紫龋琑X檢測(cè)算法是一個(gè)局部檢測(cè)算法，在局部的劃分過程中可能將全局非異常的點(diǎn)檢測(cè)為異常，這種異?？僧?dāng)作光譜不存在異常，但其空間分布異常；其次，RX 檢測(cè)是在局部正態(tài)模型的假設(shè)前提下推導(dǎo)出的一種檢測(cè)算法，對(duì)不滿足該模型的圖像數(shù)據(jù)檢測(cè)效果則難以保證，而且對(duì)大圖像來說，地物分布較復(fù)雜、地物種類較多，在無任何先驗(yàn)信息的情況下，異常的可能性非常大。針對(duì)高光譜圖像的大數(shù)據(jù)量和高維特性，在傳統(tǒng)RX異常檢測(cè)算法的基礎(chǔ)上，本文對(duì)一種改進(jìn)的星上高光譜異常檢測(cè)算法進(jìn)行了研究，采用噪聲自適應(yīng)主成分變換及局部正交子空間投影算法以降低維度，抑制干擾信息減少異常檢測(cè)處理的數(shù)據(jù)量以提高異常檢測(cè)效率和精度。

1 噪聲自適應(yīng)主成分變換

針對(duì)高光譜圖像的高數(shù)據(jù)維和大數(shù)據(jù)量導(dǎo)致的異常檢測(cè)處理的特點(diǎn)，噪聲自適應(yīng)主成分變換是一種有效的高光譜圖像光譜維變換方法[11-13]。它通過變換使大部分圖像信息以某種標(biāo)準(zhǔn)集中于少數(shù)幾個(gè)波段上，實(shí)現(xiàn)數(shù)據(jù)量壓縮。變換后的圖像會(huì)使部分在原始圖像中不可識(shí)別的特征變得明顯，從而利于目標(biāo)的快速提取。

令高光譜數(shù)據(jù)X為K×L維矩陣(K為像元數(shù)，L為波段數(shù))，且均值向量E(X)=m。m為常量，便于分析一般取其值為0，協(xié)方差矩陣C(X)=Σ。設(shè)

(1)

式中：S，N分別為信號(hào)與噪聲分量，且S，N不相關(guān)。則有

(2)

式中：ΣS，ΣN分別為S，N的協(xié)方差矩陣。設(shè)噪聲為加性噪聲，以噪聲方差與該波段總方差的比表征噪聲比例，即

(3)

由文獻(xiàn)[14]，定義信噪比(SNR)為

(4)

由式(3)、(4)可得：當(dāng)信號(hào)能量不變時(shí)，噪聲比例與圖像信噪比成反比。為提高計(jì)算效率，將噪聲自適應(yīng)主成分變換定義為線性變換，則有

(5)

Zi的噪聲比例在所有正交于Zj(j=1，2，…，i-1)的成分中最大，同時(shí)將ai標(biāo)準(zhǔn)化為(ai)TΣai=1，則噪聲自適應(yīng)主成分變換可表示為

Z=ATX

(6)

式中：A為線性變換系數(shù)陣，且A=[a1a2…ap]，是Σ-1ΣN的特征向量矩陣，則有

(7)

式中：對(duì)角線矩陣Λ為特征值矩陣。令第i個(gè)元素為特征值λi，方程求解后可得第i個(gè)波段的 SNR為

(8)

將噪聲自適應(yīng)主成分變換的成分按信噪比的大小排序，則噪聲自適應(yīng)主成分變換受不同的噪聲影響將產(chǎn)生不同的變換結(jié)果。試驗(yàn)證明：該方法的優(yōu)點(diǎn)是經(jīng)光譜維或空間維變換直接將高維數(shù)據(jù)降低到所需維數(shù)，整個(gè)降維過程時(shí)間很短；缺點(diǎn)是對(duì)原始高光譜圖像進(jìn)行了若干變換，圖像原有的部分特性發(fā)生改變。

2 局部正交子空間投影算法

局部正交子空間投影方法是將像元矢量投影到正交于干擾特征的子空間，從而實(shí)現(xiàn)最小均方誤差意義的最優(yōu)干擾特征抑制。一旦干擾特征被消除，將殘差信號(hào)投影到感興趣的特征上，使信噪比達(dá)到最大，并產(chǎn)生一個(gè)單一圖像，該圖像便可作為目標(biāo)檢測(cè)和識(shí)別的依據(jù)。

本文算法首先對(duì)局部數(shù)據(jù)建模，然后通過比較中心點(diǎn)與其鄰域內(nèi)像元光譜的差異性，得到被檢測(cè)點(diǎn)的異常度。因在小范圍鄰域內(nèi)地物類型相對(duì)單一，高光譜異常檢測(cè)常用于自然背景中的目標(biāo)提取，故以被檢測(cè)點(diǎn)作為異常目標(biāo)，用被檢測(cè)點(diǎn)鄰域內(nèi)樣本點(diǎn)均值構(gòu)造局部投影算子。設(shè)被檢測(cè)點(diǎn)為d，鄰域內(nèi)樣本點(diǎn)的某種統(tǒng)計(jì)量為U，可通過增強(qiáng)被檢測(cè)點(diǎn)特征和抑制背景實(shí)現(xiàn)異常目標(biāo)檢測(cè)。若被檢測(cè)點(diǎn)與背景相似，則在抑制背景的過程中，被檢測(cè)點(diǎn)特征同時(shí)被削弱，輸出灰度值較低；否則，經(jīng) OSP 算子投影后被檢測(cè)點(diǎn)與背景輸出值差異較大，判定為異常。

基于上述分析，可將描述觀測(cè)信息r的線性混合模型簡(jiǎn)化為

(9)

式中：α，β分別為目標(biāo)端元和背景端元的系數(shù)矩陣[15]。

(10)

將被檢測(cè)點(diǎn)光譜向量乘以投影算子，其異常度

(11)

用內(nèi)積符號(hào)表示向量乘法運(yùn)算，為提高計(jì)算效率，令系數(shù) ακ =1，式(11)可化簡(jiǎn)為

(12)

(13)

由此可知：被檢測(cè)點(diǎn)異常度正比于被檢測(cè)點(diǎn)與背景均值向量的相似程度。兩者相似度越小，θ越大，經(jīng) LOSP 投影后異常度E越大；反之，E則越小。由此，從原理上證明了LOSP 異常檢測(cè)算子的有效性。

通過試驗(yàn)證明：與OSP算法相比，本文算法的優(yōu)點(diǎn)是無需提前獲取待檢測(cè)目標(biāo)的特征信息，缺點(diǎn)是采用了較理想的線性混合模型表示地物特征信息，不能完全準(zhǔn)確反映局部的地物特征，從而會(huì)導(dǎo)致部分虛警的產(chǎn)生。

3 基于RX改進(jìn)的異常檢測(cè)算法

基于RX改進(jìn)的異常檢測(cè)算法流程如圖1所示。其過程如下。

首先，對(duì)原始的高光譜圖像X，用噪聲自適應(yīng)主成分變換方法進(jìn)行變換，獲得變換矩陣和變換后的高光譜圖像數(shù)據(jù)。噪聲自適應(yīng)主成分變換能同時(shí)對(duì)角化高光譜數(shù)據(jù)的噪聲協(xié)方差和總的協(xié)方差。設(shè)VT(i)，VN(i)分別為第i個(gè)波段數(shù)據(jù)的總協(xié)方差和噪聲協(xié)方差；L為總的波段數(shù)。假定高光譜數(shù)據(jù)僅由信號(hào)和噪聲兩種成分組成，且這些成分間相互獨(dú)立，則由式(8)可計(jì)算出第i個(gè)波段的信噪比。

圖1 算法實(shí)現(xiàn)流程Fig.1 Flowchart of proposed algorithm

其次，設(shè)定一個(gè)信噪比閾值RSNRth。閾值選取的原則為：大于閾值的噪聲自適應(yīng)主成分變換波段含有較高的信號(hào)和較低的噪聲。選取大于閾值的信噪比對(duì)應(yīng)波段的噪聲自適應(yīng)主成分變換矩陣U代表圖像中的背景信息。然后，用U構(gòu)造正交投影算子Q。U實(shí)際上包含了原始高光譜圖像數(shù)據(jù)中信息的絕大部分，用U得到的正交投影算子Q構(gòu)成了一個(gè)正交子空間，大概率背景信息等強(qiáng)干擾信息在這個(gè)子空間中被抑制，從而使小概率異常目標(biāo)信息變得突出。將噪聲自適應(yīng)主成分變換后的高光譜圖像數(shù)據(jù)X投影到該子空間。經(jīng)局部正交子空間投影變換，小概率的異常目標(biāo)在高光譜圖像中將表現(xiàn)為較亮，大概率的背景干擾信息則受到一定的削弱而變暗，使小目標(biāo)信息易被檢測(cè)到。

對(duì)經(jīng)噪聲自適應(yīng)主成分變換和OSP后得到高光譜殘差數(shù)據(jù)，再次用空域?yàn)V波方法進(jìn)行變換，使殘差數(shù)據(jù)中包含的信號(hào)能量進(jìn)一步集中于前幾個(gè)主成分，同時(shí)削弱噪聲等干擾信息的影響。因空域?yàn)V波方法得到的各主成分與各波段的方差相關(guān)聯(lián)，故此處需設(shè)定合適的方差閾值，選定大于閾值的主成分用于后續(xù)的高光譜圖像數(shù)據(jù)處理。由空域?yàn)V波方法得到的波段子集是一個(gè)異常目標(biāo)更突出的數(shù)據(jù)集，將該數(shù)據(jù)集作為中值濾波器的輸入，進(jìn)行濾波變換處理。在中值濾波器對(duì)高光譜圖像的處理過程中，從圖像空間上看，高光譜圖像數(shù)據(jù)是地物在不同波長(zhǎng)條件下的輻射值圖像。最后，將經(jīng)中值濾波處理后的高光譜圖像數(shù)據(jù)輸入，求出RX異常檢測(cè)算子，通過設(shè)定閾值判決后，就可得所需的異常目標(biāo)檢測(cè)結(jié)果。

基于前文分析的各類高光譜圖像處理技術(shù)、算法實(shí)現(xiàn)原理和流程，基于RX改進(jìn)的異常檢測(cè)算法的步驟如下。

a)輸入需處理的高光譜圖像。

b)用噪聲自適應(yīng)主成分變換對(duì)原始高光譜圖像進(jìn)行變換，根據(jù)變換后產(chǎn)生的各波段的信噪比值，設(shè)定合適的信噪比閾值，選擇大于閾值的噪聲自適應(yīng)主成分變換后的波段對(duì)應(yīng)的噪聲自適應(yīng)主成分變換矩陣，構(gòu)成正交子空間投影算子Q。

c)將噪聲自適應(yīng)主成分變換后的高光譜圖像數(shù)據(jù)構(gòu)成的X投影到正交子空間，使包括背景等強(qiáng)干擾信息在此子空間中被抑制，從而使小概率的異常目標(biāo)信息變得突出。得到高光譜圖像殘差數(shù)據(jù)。

d)用空域?yàn)V波方法對(duì)高光譜圖像殘差數(shù)據(jù)進(jìn)行特征提取，使異常目標(biāo)的能量集中于前面幾個(gè)主成分，并設(shè)定合適的閾值，確定需選擇的主成分?jǐn)?shù)量。

e)對(duì)步驟d)得到的波段子集用中值濾波器濾波，盡量消除數(shù)據(jù)中含有的各類噪聲干擾。

f)結(jié)合RX異常檢測(cè)算子對(duì)濾波后的波段子集進(jìn)行異常目標(biāo)檢測(cè)。

g)輸出的檢測(cè)結(jié)果為一灰度圖像，為使結(jié)果更可視化，用設(shè)定閥值的方法將灰度圖像轉(zhuǎn)為二值圖像，異常目標(biāo)對(duì)應(yīng)白色像元點(diǎn)。

4 仿真與分析

為驗(yàn)證算法的有效性，將美國(guó)海軍某機(jī)場(chǎng)原始高光譜的偽彩色圖像作為原始圖像，與分別用RX算法和本文算法處理獲得的結(jié)果進(jìn)行比較。原始高光譜圖像空間大小400像元×400像元，截取其含有異常目標(biāo)較多的部分進(jìn)行實(shí)驗(yàn)，其中：飛機(jī)38架，波段數(shù)126個(gè)，如圖2所示。

圖2 原始圖像Fig.2 Original image

針對(duì)原始圖像的特點(diǎn)，選取合適的閾值，分別用傳統(tǒng)RX異常檢測(cè)算法和本文改進(jìn)的異常檢測(cè)算法處理后的結(jié)果如圖3所示。圖3統(tǒng)計(jì)結(jié)果如下。

圖3 異常檢測(cè)算法處理后的結(jié)果Fig.3 Results processed by anomaly detection algorithm

a)傳統(tǒng)RX異常檢測(cè)算法檢測(cè)出異常目標(biāo)23個(gè)，目標(biāo)所占像元數(shù)233，虛警所占像元數(shù)267，平均處理時(shí)間3 562 s。

b)本文改進(jìn)的異常檢測(cè)算法檢測(cè)出異常目標(biāo)38個(gè)，目標(biāo)所占像元數(shù)360，虛警所占像元數(shù)140，平均處理時(shí)間348 s。

由此可知：與傳統(tǒng)RX異常檢測(cè)算法相比，用本文的改進(jìn)異常檢測(cè)算法處理使高光譜數(shù)據(jù)的高維性、波段間的非線性和冗余性等不利條件得以改善，同時(shí)抑制了背景噪聲干擾對(duì)目標(biāo)的影響，最終可在虛警率較低情況下檢測(cè)到全部38個(gè)異常目標(biāo)，且計(jì)算時(shí)間也有大幅減少，提高了檢測(cè)效率和準(zhǔn)確率。

5 結(jié)論

本文針對(duì)傳統(tǒng)RX異常檢測(cè)算法導(dǎo)致的星上高光譜圖像異常檢測(cè)難度大、準(zhǔn)確率低、計(jì)算時(shí)間長(zhǎng)等問題，在傳統(tǒng)RX異常檢測(cè)算法的基礎(chǔ)上，采用噪聲自適應(yīng)主成分變換及局部正交子空間投影算法對(duì)傳統(tǒng)RX異常檢測(cè)算法進(jìn)行改進(jìn)，極大地提高了星上異常檢測(cè)效率和準(zhǔn)確性。與傳統(tǒng)RX算法相比，本文算法的優(yōu)點(diǎn)有：可更準(zhǔn)確地區(qū)分異常目標(biāo)和背景干擾信息，提高了算法檢測(cè)精度，降低了虛警率；計(jì)算復(fù)雜度極低，用極短的時(shí)間就可完成目標(biāo)檢測(cè)，算法的實(shí)時(shí)性較好，利于算法在時(shí)效性要求較高的情況下使用，非常適于工程應(yīng)用。但目前該算法還未進(jìn)行全面的試驗(yàn)驗(yàn)證，異常檢測(cè)效果僅在特定工程應(yīng)用中完成了驗(yàn)證，對(duì)邊界條件設(shè)定尚未進(jìn)行深入研究，具體的適用范圍尚需后續(xù)研究認(rèn)定。后續(xù)研究需探討方法的普適性，并可考慮對(duì)該算法的繼續(xù)完善，以進(jìn)一步提高高光譜異常檢測(cè)的自適應(yīng)性和實(shí)時(shí)性。

[1] 麻永平, 張煒, 劉東旭. 高光譜偵察技術(shù)特點(diǎn)及其對(duì)地面軍事目標(biāo)威脅分析[J]. 上海航天, 2012, 29(1): 37-41.

[2] 康倩, 于晉, 林軍. 星載高光譜成像儀數(shù)據(jù)地面預(yù)處理系統(tǒng)設(shè)計(jì)[J]. 航天器工程, 2011, 20(2): 97-101.

[3] 濮建福, 裴加軍, 張寧. 基于CCSDS的高光譜壓縮空譜聯(lián)合FPGA設(shè)計(jì)與實(shí)現(xiàn)[J]. 上海航天, 2015, 32(6): 53-57.

[4] 相里斌, 王忠厚, 劉學(xué)斌. 環(huán)境減災(zāi)-1A衛(wèi)星空間調(diào)制型干涉光譜成像儀技術(shù) [J]. 航天器工程， 2009, 18(6): 43-49.

[5] 袁崇謙, 周建勛. 衛(wèi)星遙感技術(shù)在油氣勘探中的應(yīng)用[J]. 海相油氣地質(zhì), 2010, 15(2): 69-75.

[6] 穆欣侃, 陳永紅, 羅海波. 星載高光譜圖像實(shí)時(shí)處理系統(tǒng)[J]. 火力與指揮控制, 2016, 41(1): 22-27.

[7] 王彩玲, 胡柄樑, 王洪偉. 約束最大相關(guān)系數(shù)的高光譜影像目標(biāo)探測(cè)研究[J]. 激光與紅外, 2016, 46(1): 98-102.

[8] 林穎, 徐衛(wèi)明, 袁立銀, 等. 熱紅外高光譜系統(tǒng)信號(hào)成分分析及處理[J]. 激光與紅外, 2010, 40(12): 1324-1329.

[9] 廖佳俊, 劉志剛, 姜江軍. 基于稀疏表示分步重構(gòu)算法的高光譜目標(biāo)檢測(cè)[J]. 紅外技術(shù), 2016, 38(8): 699-704.

[10] REED I S, YU X. Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution[J]. IEEE Transactions on Acoustic Speech and Signal Process, 1990, 38(10): 1760-1770.

[11] ROGER R E. A fast way to compute the noise-adjusted principal components transform matrix[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994, 32(11): 1194-1196.

[12] CHANG C I, DU Q. Interference and noise-adjusted principal components analysis[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(5): 2387-2396.

[13] CHANG C I, DU Q. A joint band prioritization and band decorrelation approach to band selection[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(6): 2631-2641.

[14] SD W J, BS G, HL E. Anomaly detection from hyperspectral imagery[J]. IEEE Signal Processing Magazine, 2002, 19(1): 58-69.

[15] 李智勇. 高光譜圖像異常檢測(cè)方法研究[D]. 長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué), 2004.

Study on Improved Algorithm of Hyperspectral Anomaly Detection on Satellite

LI Ke, YAO Chong-bin, XU Hong-xin, XIE Bao-rong, SHANG Ji-yang

(Shanghai Aerospace Electronic Technology Institute, Shanghai 201109, China)

To solve the problems of the low accuracy and high expense of hyperspectral anomaly detection on-borne using traditional RX anomaly detection algorithm, an improved algorithm of hyperspectral anomaly detection was studied in this paper. Based on traditional RX anomaly detection algorithm, the original image was transferred by noise adaptive principal component transform. The data which single-to-noise ratio was bigger than threshold set would be selected. The selected data space was projected to orthogonal subspace to obtain the characteristics of the residual error of the image by local orthogonal subspace projection algorithm, so the anomaly target information with low probability could the outstand by suppressing the background and other strong interference. The characteristics of the residual error data were picked up by space-domain filtering method which could make majority information of the original image concentrate on few wave bands under certain standard. The subspace obtained was filtered to eliminate the noise by mean value filter. The anomaly detection of the filtered subspace was carried on using RX anomaly detection algorithm. An experiment showed the effectiveness of the method proposed. The method proposed can reduce the information in order to improve the accuracy and efficiency of anomaly detection onborne by data compressing, dimension reducing and disturbance suppressing.

on-borne detection; hyperspectral; anomaly detection; local orthogonal subspace projection (LOSP); RX anomaly detection algorithm; noise adaptive principal component transform; signal-to-noise ratio; threshold

1006-1630(2017)03-0047-05

2016-12-25；

2017-03-12

國(guó)家青年科學(xué)基金項(xiàng)目資助(F050405)

李 恪(1980—)，男，博士，主要從事光電遙感技術(shù)研究。

TP731

A

10.19328/j.cnki.1006-1630.2017.03.006