隋中山，李俊山，張 姣，隋曉斐

(1.火箭軍工程大學(xué)，信息工程系，陜西 西安 710025；2.96618部隊，北京 100085)

?

基于張量低秩分解和稀疏表示的紅外微小氣體泄漏檢測

隋中山1*，李俊山1，張 姣1，隋曉斐2

(1.火箭軍工程大學(xué)，信息工程系，陜西 西安 710025；2.96618部隊，北京 100085)

為了檢測石化工業(yè)生產(chǎn)過程中微小氣體的泄漏，提出了一種應(yīng)用紅外成像技術(shù)的單幀紅外小目標(biāo)檢測方法。研究了低秩稀疏分解理論和稀疏表示理論，并提出了一種新的基于張量低秩分解和稀疏表示的小目標(biāo)檢測方法。該方法基于張量分解的形式充分發(fā)掘背景矩陣所包含的信息；利用先驗知識構(gòu)造微小氣體泄漏的目標(biāo)字典；同時利用背景的低秩約束和小目標(biāo)的稀疏表示約束分解出微小氣體的泄漏目標(biāo)。最后基于非精確增廣拉格朗日乘子法(IALM)，對本文算法進行最優(yōu)化求解，并通過實驗分析比較了本文方法和已有方法的優(yōu)缺點。結(jié)果表明：本文方法的檢測效果優(yōu)于其他已有方法，并且具有較好的ROC(受試者工作特征)曲線， 可以滿足工業(yè)生產(chǎn)中對微小氣體泄漏檢測的要求。

計算機視覺；紅外檢測；泄漏檢測；張量低秩分解；稀疏表示；紅外成像

1 引 言

石化工業(yè)中很多中間或最終產(chǎn)品都是氣體，這些氣體大部分是無色透明的且易發(fā)生泄漏，易造成重大安全事故。目前，傳統(tǒng)的石化工業(yè)泄漏氣體的檢測方式大都基于化學(xué)方法。近年來，紅外成像技術(shù)不斷進步，紅外成像設(shè)備成本逐漸降低，使得利用紅外圖像處理和計算機視覺原理進行氣體泄漏的自動檢測已經(jīng)成為可能。自然界中一切溫度高于絕對零度的物體都在不停地發(fā)射紅外輻射，而且基于此可進行熱紅外成像，在熱紅外圖像中不同的亮度區(qū)域代表不同的溫度區(qū)域。在石化工業(yè)的生產(chǎn)過程中，工業(yè)氣體的溫度一般與周邊自然環(huán)境有較大區(qū)別，則當(dāng)管道或者容器存在微小泄漏時，泄漏點的溫度將與周邊溫度明顯不同。在紅外圖像中，微小氣體本身的形狀和內(nèi)部結(jié)構(gòu)很難體現(xiàn)出來，但是它與周圍環(huán)境之間的溫度差而產(chǎn)生的亮度差則很明顯。因此可以將微小氣體檢測問題歸為紅外小目標(biāo)檢測問題，小目標(biāo)是指在圖像中所占像素較少，內(nèi)部結(jié)構(gòu)和外部輪廓均不明顯的目標(biāo)[1-2]。

目前紅外小目標(biāo)檢測主要分為兩類：基于單幀圖像的方法和基于圖像序列的方法。就本文研究背景而言，氣體泄漏位置在視頻序列中不會發(fā)生變化，所以基于圖像序列的方法并不適用。因此本文將重點研究基于單幀圖像的紅外小目標(biāo)檢測算法。該算法主要分為3種：濾波方法[3-4]、基于人眼視覺注意機制的顯著性檢測方法[5-6]、基于成分表達和字典學(xué)習(xí)的方法[7-9]。濾波方法在局部信噪較低、頻率特性不顯著的情況下誤警率比較高；顯著性檢測方法的理論基礎(chǔ)尚待進一步完善且結(jié)果相對不穩(wěn)定；基于成分表達和字典學(xué)習(xí)的方法結(jié)果較為穩(wěn)定、效果較好，近年來越來越受到重視。

2 相關(guān)方法和理論

2.1 低秩稀疏分解

低秩稀疏分解方法在前景檢測領(lǐng)域已經(jīng)取得了重要進展[10-12]。圖像構(gòu)成的矩陣可以分解成低秩矩陣代表的背景和稀疏矩陣代表的前景對象。它可以表示如下：

D=L+S，

(1)

其中：D∈Rm×n是所觀察到的圖像，當(dāng)然D必須滿足非局部自相似性的條件，L和S分別表示背景和前景?；谥鞒煞肿粉?Principal Component Pursuit，PCP)的魯棒主成分分析(Robust Principal Component Andysis，RPCA)首次提出用l1-范數(shù)來約束前景矩陣，因為它們是包含一小部分非零元素的稀疏矩陣，并且假設(shè)背景圖像線性相關(guān)，從而形成低秩矩陣L。這個矩陣分解可以通過以下凸優(yōu)化解決：

(2)

其中：‖L‖*是矩陣L的核范數(shù)，是其奇異值的總和，‖S‖1表示S的l1-范數(shù)。λ為一個正的系數(shù)。

以上為低秩稀疏分解的基本形式，隨后又發(fā)展出多種不同形式：IRPCA(Inductive Robust Principal Component Analysis)、低秩表示(Low-Rank Representation，LRR)、LatLRR(Latent Low-Rank Representation)等方法。其中LRR方法應(yīng)用較為廣泛，因為它可以處理背景數(shù)據(jù)來自多個不同子空間的情況，其核心思想是找到一個低秩表達矩陣V=[v1，v2，…，vn]∈Rn×n，它可以將數(shù)據(jù)向量用數(shù)據(jù)字典A的基的線性組合表達出來，如果將原始數(shù)據(jù)D直接作為數(shù)據(jù)字典(即A=D)，則LRR方法可以用下列最優(yōu)化問題表示：

(3)

當(dāng)然LRR方法也存在一些缺點如：LRR方法在恢復(fù)或分割數(shù)據(jù)時會丟失原始數(shù)據(jù)矩陣列方向上包含的信息；LRR方法要求原始數(shù)據(jù)D必須包含足夠多的純凈信息，但這在現(xiàn)實生活中是較難滿足的。

2.2 稀疏表示

生物學(xué)家的研究成果顯示，哺乳類動物通過視覺神經(jīng)能快速、準(zhǔn)確、低耗地對自然景物成像，其核心在于它們感知視覺信息時，只使用視覺皮層中少量神經(jīng)元即可完成任務(wù)，即視覺信息可以用少量的神經(jīng)元進行稀疏表示。從信號處理角度看，即信號可以表示為字典中少量原子的線性組合。在圖像處理應(yīng)用的過程中，其表達形式如下：

T=OH，

(4)

其中：O={o1,o2,…,op}∈Rq×p為字典矩陣，每一列oi為字典原子，若滿足p?q，則O為超完備字典；H為稀疏表示系數(shù)矩陣，且必須滿足‖H‖0<γ，γ為一較小常數(shù)。其結(jié)果可以通過將‖H‖0松弛為‖H‖1來求解。

2.3 低秩稀疏分解和稀疏表示的融合

目前，有研究者嘗試在一定條件下對兩種方法進行融合。他們利用背景的低秩特性，同時構(gòu)造字典對前景的稀疏表達。針對紅外小目標(biāo)檢測，將LRR和稀疏表示的方法相融合，如下所示：

D=DV+OH，

(5)

其中：D為原始數(shù)據(jù)，O為紅外小目標(biāo)超完備字典，需要先行構(gòu)造，且對算法性能影響較大。這個矩陣分解可以通過以下凸優(yōu)化解決：

(6)

雖然融合方法在性能和效果上有一定提升，但是原始LRR方法的一些缺陷仍然沒有解決，因此本文將對其進行改進。

3 張量低秩分解和稀疏表示

3.1 張量低秩分解和稀疏表示

LRR方法是沿著原始數(shù)據(jù)矩陣的行方向恢復(fù)數(shù)據(jù)的，所以列方向上的信息將丟失。為了能同時利用行信息和列信息，本文引入張量的表示形式。尋求一個三階張量的表達形式UDV來代替LRR方法中的低秩部分DV，那么張量低秩分解和稀疏表示融合的表達形式為：

D=UDV+OH，

(7)

其中：D為原始數(shù)據(jù)，U為一個低秩投影矩陣，V為一個低秩表示矩陣，O為紅外小目標(biāo)超完備字典，H為一個稀疏的系數(shù)矩陣，即須滿足‖H‖0<γ，γ為一較小常數(shù)。為了完成式(7)所示的分解，結(jié)合所須條件，需要實現(xiàn)下列目標(biāo)：

s.t.D=UDV+OH.

(8)

對式(8)進行松弛，結(jié)果如下：

s.t.D=UDV+OH.

(9)

3.2 最優(yōu)化求解

本文采用非精確增廣拉格朗日乘子(IALM)方法進行最優(yōu)化求解，首先引入兩個輔助變量J，K，將式(9)轉(zhuǎn)變?yōu)橄率剑?/p>

s.t.D=XV+OH，X=UD，U=J，V=K.

(10)

其增廣拉格朗日方程為：

L(U，V，H，J，K，Y1，Y2，Y3)=

‖J‖*+‖K‖*+λ‖H‖1+

〈Y1，D-UDV-OH〉+〈Y2，U-J〉+〈Y3，V-K〉+

(11)

其中：Y1，Y2，Y3是拉格朗日乘子，μ為懲罰因子(μ>0)。

其具體求解步驟為：

(1)初始化：輸入原始數(shù)據(jù)，對各變量賦初值，設(shè)定誤差。

(2)固定其他變量更新J：

(12)

SVTμ(Y)=Pdiag[(δ-μ)+]Q

(13)

其中：P，Q，δ由Y的奇異值分解獲得，即Y=PΣQ，Σ=diag(δ)。

(3)固定其他變量更新K：

(14)

式(14)的求解方法與式(12)相同。

(4)固定其他變量，更新H：

H=argminβ‖H‖1+〈Y1，D-UDV-OH〉+

(15)

式(15)的解可以通過二次最小代價流的方法獲得。

(5)固定其他變量更新U：

令?L(U，V，H，J，K，Y1，Y2，Y3)/?U=0，求得：

(DVVTDT+I)-1.

(16)

(6)固定其他變量更新V：

令?L(U，V，S，E，J，K，Y1，Y2，Y3)/?V=0，求得：

V=(DTUTUD+I)-1·

(17)

(7)更新拉格朗日乘子Y1，Y2，Y3：

Y1=Y1+μ(D-UDV-OH)，

(18)

Y2=Y2+μ(U-J)，

(19)

Y3=Y3+μ(V-K).

(20)

(8)更新參數(shù)μ：

μ=min(κμ，μmax)，

(21)

其中：κ為更新步長且κ>1，μmax為指定的μ的上限。

(9)判斷是否達到收斂條件：

‖D-UDV-OH‖∞<ε，

‖U-J‖∞<ε，‖V-K‖∞<ε，

(22)

其中ε為誤差。

當(dāng)達到收斂條件時算法結(jié)束，否則重復(fù)步驟(2)～(9)。

3.3 算法分析

本文方法相對于第2節(jié)介紹的已有方法的最大優(yōu)點在于算法在恢復(fù)背景數(shù)據(jù)矩陣(UDV)是從行(V)和列(U)兩個方向同時進行計算的，因此本文算法可以同時利用行方向和列方向上的信息進行恢復(fù)；其次，本文算法將背景的低秩約束和小目標(biāo)的稀疏表示約束相結(jié)合，傳統(tǒng)的低秩分解重點對背景進行約束，本文算法則可以更好地突出前景小目標(biāo)的重要性，也更適合于本文應(yīng)用背景。本文算法的復(fù)雜度為O(m2n+n3)，較為復(fù)雜，在實時性要求較高的場景中應(yīng)用還存在一定困難。

4 張量低秩分解和稀疏表示的紅外

微小氣體泄漏檢測算法

利用第3節(jié)提出的方法進行泄漏檢測，其基本步驟如圖1所示。

圖1 張量低秩分解和稀疏表示的紅外微小氣體泄漏檢測算法

Fig.1 Micro gas leakage detection algorithm based on tensor low rank decomposition and sparse representation

4.1 構(gòu)造微小氣體泄漏超完備字典

超完備字典的選擇與構(gòu)造是一個基礎(chǔ)且關(guān)鍵的問題，影響著算法的最終效果。利用大量真實樣本進行字典學(xué)習(xí)和更新能取得最佳效果，但是在本文實驗條件下很難獲得足夠數(shù)量的泄漏樣本。現(xiàn)實中發(fā)生泄漏時，泄漏中心點往往是溫度的極值點(低溫氣體泄漏對應(yīng)溫度最低點，高溫氣體泄漏對應(yīng)溫度最高點)。以高溫氣體為例，當(dāng)泄漏發(fā)生時，泄漏中心點的溫度最高，往四周擴散時逐漸降低。體現(xiàn)在紅外圖像上即是泄漏中心點亮度最高，往四周逐漸降低。因此本文選擇二維高斯模型(如式(23)所示)對微小氣體泄漏進行建模，通過調(diào)節(jié)模型中的參數(shù)，生成一系列樣本圖像，繼而構(gòu)成微小氣體泄漏超完備目標(biāo)字典：

(23)

其中：(i，j)為圖像點的坐標(biāo)；Imax是泄漏中心點的亮度值；(x0，y0)是泄漏中心點坐標(biāo)；σx和σy分別為水平散布參數(shù)和垂直散布參數(shù)。通過調(diào)節(jié)Imax，(x0，y0)，σx和σy4個參數(shù)，即可獲得不同形狀，不同中心點位置，不同亮度的氣體泄漏樣本圖像。α，β表示Imax和σx，σy的約束關(guān)系。最后將得到的k×k氣體泄漏樣本圖像按列展開成k2維向量并組合成字典矩陣O={o1,o2,…,op}∈Rk2*p，其中的每一個列向量oi是超完備字典的基，總共有p個原子。圖2為高溫氣體泄漏檢測所使用的超完備字典。

圖2 微小氣體泄漏超完備字典

4.2 構(gòu)造數(shù)據(jù)矩陣

由于在算法中需要使用超完備字典O，所以需要將原始圖像重新構(gòu)造成與超完備字典O相匹配的數(shù)據(jù)矩陣。其方法為使用滑動窗口進行掃描，窗口大小一般與氣體泄漏目標(biāo)大小匹配，為k×k，移動步長為k/2，共獲得N個圖像子塊，相鄰子塊之間有重疊部分可以強化其低秩屬性，將每個子塊按列展開成k2維向量，組成數(shù)據(jù)矩陣D∈Rk2×N。圖3為原始圖像與數(shù)據(jù)矩陣對比圖。

圖3 原始圖像與數(shù)據(jù)矩陣

4.3 運行算法獲得結(jié)果

將重構(gòu)后的數(shù)據(jù)矩陣D與超完備字典O作為輸入，通過第3.2節(jié)中提出的最優(yōu)化方法進行求解，最終得到最優(yōu)估計U*、V*和H*，計算OH*得到微小氣體泄漏檢測的初步結(jié)果，然后對OH*構(gòu)造數(shù)據(jù)矩陣過程的逆處理(如4.2節(jié)所述)，由于滑動窗口的步長小于滑動窗口本身大小，所以在原始圖像中某一個位置點有可能會對應(yīng)多個子窗口中不同的亮度值，當(dāng)出現(xiàn)這種情況時就取各個不同亮度值的平均值，最終即可獲得檢測結(jié)果。

5 實驗與結(jié)果

設(shè)計了兩組實驗以驗證本文算法的有效性，一組為有效性實驗，另一組為對比性實驗。所有程序運行在Windows 7+Intel Core i5-3470 3.2 GHz+8G RAM配置下，采用Matlab 2010b編程。

5.1 有效性實驗

在本文應(yīng)用背景下，對算法的有效性進行驗證。

圖4 中間結(jié)果

圖4展示了算法在運行過程中的中間結(jié)果，(1)為原始圖像，(2)為原始圖像轉(zhuǎn)置后的數(shù)據(jù)矩陣，(3)為U*DV*，(4)為OH*，最終結(jié)果如圖5所示。

圖5中，A，B，D為高溫氣體發(fā)生泄漏的場景，C為低溫氣體發(fā)生泄漏的場景。第(1)列為原始圖像，第(2)列為算法得到的背景成分，第(3)列為檢出的氣體泄漏信息。需要指出的是對高溫氣體泄漏和低溫氣體泄漏檢測兩種情況，需要構(gòu)造不同的字典。通過實驗發(fā)現(xiàn)，本文方法可以有效檢出高溫和低溫微小氣體泄漏。

圖5 最終檢測結(jié)果

5.2 對比性實驗

圖6中，第(1)列為原始圖像。第(2)列為本文算法得到的背景成分。第(3)列為高低帽算法[3]的檢測結(jié)果，它是一種濾波算法，其中A、B、D行為采用TH(高帽)算法檢測暗背景下的高亮度目標(biāo)，利用開操作將暗背景下的小亮目標(biāo)侵蝕掉，再與原圖像相減，突出小目標(biāo)，其適用于平滑背景。對于噪聲大，局部對比不明顯的圖像，如D(3)的檢測效果不佳；C行為采用LH(低帽)算法檢測暗目標(biāo)，其利用閉操作將小的暗目標(biāo)填充掉，再與原圖像相減?？梢悦黠@看出檢測目標(biāo)的亮度發(fā)生了反轉(zhuǎn)，由暗目標(biāo)轉(zhuǎn)變?yōu)榱聊繕?biāo)。第(4)列為基于視覺顯著性[6]的檢測算法，其利用亮度信息和方向信息進行檢測，它可以同時檢測出亮度異常的點，在原始圖像中較亮和較暗的目標(biāo)在檢測結(jié)果中都呈現(xiàn)為亮點，如D(4)所示，因此容易造成誤檢；當(dāng)目標(biāo)過小時則容易造成漏檢，如A(4)所示。由圖6可以看出：不同的算法，檢測結(jié)果差異比較明顯。在本文的研究背景下，基于視覺顯著性的檢測算法的通用性較好的，但是容易造成誤檢和漏檢；高低帽方法在簡單背景下效果較好，在復(fù)雜背景下效果不佳；本文方法充分利用了原始圖像的低秩背景和目標(biāo)信息的先驗知識，有效降低了漏檢率，對紅外小目標(biāo)的檢測效果優(yōu)于其他算法，在通用性上則有一定不足，因為針對低溫和高溫泄漏需要設(shè)計不同的字典。

圖6 檢測結(jié)果

為了定量說明本文算法的優(yōu)點，利用受試者工作特征(Receiver Operating Characteristics，ROC)曲線對各個算法性能進行評價。ROC曲線以誤警率為橫軸，以正確檢出率為縱軸，曲線與X軸之間的面積越大，說明算法性能越好，正確檢出率和誤警率定義如下：

Pd=Nr/NT，

(24)

Pfa=Nf/NT，

(25)

其中：Pd為正確檢出率，Pfa為誤警率，Nr為正確檢測個數(shù)，Nf為錯誤檢測個數(shù)，NT為實際目標(biāo)總數(shù)。

在圖5所示的4組背景下，隨機添加模擬氣體泄漏點，其檢測結(jié)果如圖7所示。

(a)

(b)

(c)

(d)

6 結(jié) 論

本文根據(jù)現(xiàn)代石化工業(yè)對微小氣體泄漏檢測的要求，在理論上提出基于張量低秩分解和稀疏表示的方法，并給出了基于IALM的最優(yōu)化求解方法。然后針對實際應(yīng)用，給出了切實可行的實施步驟和方法。實驗結(jié)果證明：本文所提出的方法在檢測效果上優(yōu)于當(dāng)前流行的其他方法。由于在實際生產(chǎn)過程中，高溫氣體泄漏和低溫氣體泄漏都需要檢測，本文算法的通用性需要進一步驗證，同時本文算法的計算效率需要進一步提高。

[1] 何玉杰，李敏，張金利，等.基于低秩三分解的紅外圖像雜波抑制[J].光學(xué) 精密工程， 2015，23(7)： 2069-2078. HE Y J， LI M， ZHANG J L，etal.. Clutter suppression of infrared image based on there-component low-rank matrix decomposition [J].Opt.PrecisionEng.， 2015，23(7)： 2069-2078.(in Chinese)

[2] 趙愛罡，王宏力，楊小岡，等.基于結(jié)構(gòu)低秩編碼的復(fù)雜環(huán)境紅外弱小目標(biāo)檢測算法[J].中國慣性技術(shù)學(xué)報，2015，23(5)：662-669. ZHAO A G， WANG H L， YANG X G，etal.. Infrared dim small target detection algorithm based on structural low-rank coding under complex environment [J].JournalofChineseInertialTechnology，2015，23(5)：662-669.(in Chinese)

[3] TOM V T， PELI T， LEUNG M，etal.. Morphology-based algorithm for point target detection in infrared backgrounds [C].OpticalEngineeringandPhotonicsinAerospaceSensing.InternationalSocietyforOpticsandPhotonics， 1993： 2-11.

[4] DESHPANDE S D， MENG H E， VENKATESWARLU R，etal.. Max-mean and max-median filters for detection of small targets [C].SPIE’sInternationalSymposiumonOpticalScience，Engineering，andInstrumentation， 1999： 74-83.

[5] HAN J H， MA Y， ZHOU B，etal.. A robust infrared small target detection algorithm based on human visual system [J].GeoscienceandRemoteSensingLetters， 2014，11(12)： 2168-2172.

[6] QI S X， MA J， TAO C，etal.. A robust directional saliency-based method for infrared small target detection under various complex backgrounds [J].IEEEGeoscienceandRemoteSensingLetters， 2013，10(3)： 495-499.

[7] YANG C W， LIU H P，LIAO S Y，etal.. Small target detection in infrared video sequence using robust dictionary learning [J].InfraredPhysics&Technology，2015， 68： 1-9.

[8] ZHENG C Y， LI H. Small infrared target detection based on low-rank and sparse matrix decomposition [J].AppliedMechanicsandMaterials， 2013， 239：214-218.

[9] HE Y J， LI M， ZHANG J L，etal.. Small infrared target detection based on low-rank and sparse representation [J].InfraredPhysics&Technology， 2015， 68： 98-109.

[10] LIU X， ZHAO G Y. Background subtraction based on low-rank and structured sparse decomposition [J].IEEETransactionsonImageProcessingDRAFT，2015.

[11] YAO J， LIU X， QI C. Foreground detection using low rank and structured sparsity [C].Proc.IEEEInt.Conf.Multimed.Expo.， 2014： 1-6.

[12] BOUWMANS T， ZAHZAH E H. Robust PCA via principal component pursuit： A review for a comparative evaluation in video surveillance [J].Comput.Vis.ImageUnderst.， 2014， 122： 22-34.

隋中山(1985-)，男，山東平度人，博士研究生，2007年、2009年于第二炮兵工程學(xué)院分別獲得學(xué)士、碩士學(xué)位，主要從事圖像目標(biāo)識別方面的研究。E-mail：zclszs@163.com

李俊山(1956-)，男，陜西白水人，教授，博士生導(dǎo)師，1981于國防科技大學(xué)獲得學(xué)士學(xué)位，1988年于第二炮兵工程學(xué)院獲得碩士學(xué)位，2001年于西安微電子技術(shù)研究所獲得博士學(xué)位，主要從事智能圖像處理與目標(biāo)感知識別，電子對抗模擬與仿真等方面的研究。E-mail：lijunshan403@163.com

(版權(quán)所有 未經(jīng)許可 不得轉(zhuǎn)載)

Micro gas leakage detection based on tensor low rank decomposition and sparse representation from infrared images

SUI Zhong-shan1*， LI Jun-shan1， ZHANG jiao1， SUI Xiao-fei2

(1.Department of Information Engineering，RocketForceUniversityofEngineering，Xi′an710025，China；2.Unit96618，theChinesePeople’sLiberationArmy，Beijing100085，China)

To detect the micro gas leakage in petrochemical production， a single-frame small target detection method was proposed by using infrared images. The low-rank sparse decomposition theory and sparse representation theory were researched and an innovative method to detect a micro-target was proposed based on tensor low-rank decomposition and sparse representation. The tensor decomposition form was employed in exploiting the information contained in background matrices， The priori knowledge was used to construct a micro gas leakage target dictionary， meanwhile， the micro-gas leakage targets were decomposed by low-rank constraint in the background and sparse representation in the micro-target. Finally， the algorithm was solved optimally by using Inexact Augmented Lagrange Multiplier(IALM) method and its merits were compared with that of common methods. The results indicate that the proposed algorithm has better detection efficiency than that of common methods and it shows better ROC (Receiver Operating Characteristics)curves. It concludes that these results meet the requirements of micro gas leakage detection during industrial productions.

computer vision； infrared detection； leakage detection； tensor low-rank decomposition； sparse representation； infrared imaging

2016-07-11；

2016-08-12.

國家自然科學(xué)基金資助項目(No.61175120)

1004-924X(2016)11-2855-08

TP391.4；TN215

A

10.3788/OPE.20162411.2855

*Correspondingauthor，E-mail：zclszs@163.com