李賡飛，李桂菊，韓廣良，劉培勛，江 山

(中國科學院 長春光學精密機械與物理研究所，吉林 長春 130033)

?

紅外成像系統(tǒng)的非均勻性實時校正

李賡飛*，李桂菊，韓廣良，劉培勛，江 山

(中國科學院 長春光學精密機械與物理研究所，吉林 長春 130033)

針對紅外相機系統(tǒng)響應非均勻性所產(chǎn)生的豎條紋固定噪聲，提出了基于多線程優(yōu)化的單幀紅外圖像非均勻性實時校正算法。該算法參考全變分理論的非均勻性校正算法確定噪聲與圖像輸出之間的加性校正關(guān)系，建立了描述條紋噪聲的數(shù)學模型；通過濾波算法將圖像分為高頻分量和低頻分量，運用建立的條紋噪聲模型，基于最小二乘法從高頻部分中擬合出條紋噪聲。最后，使用全變分理論確定的加性關(guān)系對圖像進行校正。為了提高算法的計算速度，運用多線程技術(shù)對算法進行了優(yōu)化。對提出的算法與MIRE(Midway Infrared Equalization)算法及傳統(tǒng)雙邊濾波算法進行了圖像質(zhì)量評價和性能對比。結(jié)果顯示：本文算法使圖像非均勻性較原圖降低了3%；算法效率與MIRE算法無異，系統(tǒng)運行時間在320×256的14位紅外圖像上為1.5 ms/frame，達到了工程標準。

紅外圖像；成像系統(tǒng)；豎條紋噪聲；非均勻性校正；實時處理

1 引 言

受材料和工藝水平的限制，紅外焦平面陣列(Infrared Focal Plane Arrays，IRFPA)各探測器單元之間普遍存在著響應非均勻性，此外，圖像數(shù)據(jù)讀出電路中各通道之間的差異會使這種非均勻性更加明顯[1]。紅外圖像非均勻性主要分為時變噪聲及固定噪聲兩類[2]，通過觀察靶場實際圖像發(fā)現(xiàn)，多數(shù)紅外相機均表現(xiàn)出一定的非均勻性，這些非均勻性中對成像效果影響較大的是網(wǎng)狀及條紋狀噪聲。對于紅外圖像中的垂直或水平線，難以判斷其是場景目標自身的邊緣還是由錯誤的固定模式噪聲引起的，而傳統(tǒng)的去噪方法常常會引起圖像的過模糊化，在圖像細節(jié)很低的紅外圖像中尤為突出。因此，研究如何有效去除固定噪聲實現(xiàn)非均勻校正對于提高紅外圖像質(zhì)量具有重要的實用價值。

目前針對IRFPA響應的非均勻校正算法主要包括基于標準源定標和基于場景的校正方法兩類?；趨⒖荚炊祟惙椒ㄊ窃缙谔岢龅囊环N易于實現(xiàn)且較有效的方法。典型的如一點校正法，兩點校正法，多點校正法[3-4]等。這類方法一般通過更新非均勻性校正算法的參數(shù)實現(xiàn)，且需要一個額外的電子設(shè)備，進行多次重復定標。不但增加了系統(tǒng)復雜度，還會打斷圖像采集過程。在此基礎(chǔ)上，人們提出了基于場景的非均勻性校正算法[5-9]，如基于恒定統(tǒng)計算法[5]，神經(jīng)網(wǎng)絡(luò)算法[9]。這些基于場景的算法主要根據(jù)實際場景提取校正系數(shù)，可以自適應地校正探測器因溫度和時間漂移引起的變化。為了估計出可靠的參數(shù)，這類算法需要利用一組圖片進行提前訓練，一般收斂速度較慢，并且場景的幀間差異要足夠小。如果圖像序列幀間的差異過大，則在當前圖像上會有前一幀圖像遺留的殘影，而校正殘影的問題也是最近很多新算法必須研究的課題[7-8]。上述兩類非均勻校正算法具有普遍性，對于條紋噪聲，有時達不到滿意的效果。因此一些文獻又提出了針對性更強的算法[10-14]，例如基于小波的改進算法[10]，基于全變分理論的非均勻校正算法[11]，基于雙邊濾波的算法[12]，有的甚至直接從直方圖的特征來考慮[13-14]，目前效果比較好的是文獻[14]中提出的Midway InfRared Equalization(MIRE)算法，但以上算法在處理效果和運算速度的平衡上不能令人滿意，在應用到實際工程中需要一定的取舍，有時還會使架構(gòu)變的復雜。

本文為解決實際工程中紅外相機普遍遇到的非均勻性問題，針對條紋樣式噪聲進行研究，參考全變分理論分析了圖像與噪聲的加性關(guān)系，建立了描述豎條紋噪聲的數(shù)學模型，提出了一種簡單的基于單幀的條紋狀噪聲校正算法，從根源上避免了殘影現(xiàn)象，并運用多線程技術(shù)分析優(yōu)化了算法，能夠高速實時地對圖像進行處理。

2 傳統(tǒng)去噪原理及校正模型

2.1 傳統(tǒng)圖像去噪原理

一般圖像大多采用濾波算法進行去噪，這類濾波器通過對鄰域像素的加權(quán)平均來獲取原始信號的可靠估計。其公式如下：

P(x)=k-1∫ΩI(y)fs(‖y-x‖)dy，

(1)

當圖像存在邊緣時，僅根據(jù)空間距離做處理，往往會使邊緣過度平滑，從而造成圖像模糊。而雙邊濾波器的主要思想是在空間位置上的加權(quán)操作中引入值域(像素值)加權(quán)，也就是說在分配圖像鄰域權(quán)值時，不僅考慮到了空間上是否相近，也考慮到了值域上是否相近，值域越相似，分配的權(quán)值也就越大，因此雙邊濾波也是一種非線性濾波器，加權(quán)組合后的雙邊濾波公式如下：

P(x)=k-1∫ΩI(y)fs(‖y-x‖)·

fr(I(y)-I(x))dy，

(2)

雙邊濾波器綜合考慮了值域(像素值之間的相似程度越大權(quán)值越大)與定義域(像素之間的距離越相近權(quán)值越大)，再分配權(quán)值，從而減少了去噪過程中圖像邊緣的過度平滑現(xiàn)象，在一定程度上保證了圖像邊緣。

2.2 非均勻校正模型

在紅外圖像中，由非均勻性造成的圖像退化是極其普遍的現(xiàn)象，而實際信號響應值和紅外線光子流間的關(guān)系正是非線性的。所以在多數(shù)文獻中，通用的非均勻性校正攝動模型為：

I(i，j)=g(S(i，j))+r(i，j)，

(3)

其中：(i，j)為像素坐標，I(i，j)為采集到的圖像，S(i，j)為真實紅外圖像，g(x)是關(guān)于傳感器的轉(zhuǎn)變函數(shù)，而r(i，j)則表示其他的隨機噪聲。

非均勻校正算法的主要目的就是找出針對每一個傳感器的g(x)函數(shù)，并對r(i，j)進行適當?shù)膬?yōu)化。

2.3 條紋校正模型

目前靶場上常見的紅外圖像中，以豎條紋樣噪聲對成像的影響最為顯著。為了簡化問題，將紅外相機的每一個傳感器在響應時均近似為線性模型，這樣它們的輸出形式為：

I(i，j)=g(i，j)*S(i，j)+o(i，j)+r(i，j)，

(4)

其中：(i，j)為像素坐標，I(i，j)為采集到的圖像，S(i，j)為真實紅外圖像，g(i，j)和o(i，j)兩個函數(shù)分別表示第(i，j)個傳感器的增益和偏移，而r(i，j) 表示其他的隨機噪聲。

3 實時紅外圖像的非均勻性校正算法

3.1 紅外噪聲模型的建立

全變分理論將圖像定義為一個二維空間上的有界變差函數(shù)，并假設(shè)非均勻條紋噪聲是一個固定的加性噪聲，因此采集到的圖像I(i，j)可分為兩部分：

I(i，j)=S(i，j)+N(i，j)，

(5)

其中：N(i，j)為條紋噪聲，S(i，j)為真實的紅外場景，I(i，j)為采集到的紅外圖像。

鑒于全變分理論的思想，可以將式(5)代入式(4)得到新的式(6)：

N(i，j)=(g(i，j)-1)*S(i，j)+o(i，j)+r(i，j).

(6)

由于豎條紋噪聲屬于固定噪聲，所以公式(6)中的增益g(i，j)和偏移函數(shù)o(i，j)在時間上的變化是十分緩慢的，可以認為g(i，j)和o(i，j)為常數(shù)。針對條紋噪聲，假設(shè)隨機噪聲為零，可得條紋噪聲公式為：

N(i，j)=a1*I(i，j)+a0.

(7)

雖然公式(7)為線性模型。但根據(jù)多項式的形式，不難將其擴展為高階的非線性模型，即：

(8)

其中：N(i，j)為豎條紋噪聲，I(i，j)為采集到的圖像，ak(k=1，2，3，…)為誤差參數(shù)，k為多項式階數(shù)。

若對圖像中的像素逐一進行處理，則對每個像素都需要一組參數(shù)ak(k=1，2，3，…)來進行校正，導致算法復雜度極高，因此本文針對豎條紋固定噪聲，采用以列為整體的處理策略，對公式(8)增加簡單約束：

(9)

其中：j表示圖像的任意列，這樣針對每一列，求出一組參數(shù)ak即可。

3.2 單一圖像的非均勻性校正算法

綜上所述，本文采用的算法流程如圖1所示。由式(2)分離出高頻分量和低頻分量，具體步驟為：高頻分量包含了圖像部分邊緣細節(jié)以及豎條紋固定模式噪聲，以原始圖像作為輸入，通過式(9)與高頻分量擬合，得出純條紋噪聲圖像，最后將條紋噪聲帶入式(5)得到校正后的圖像。

圖1 NUC算法流程

算法擬合過程中，以列為單位通過尋找最優(yōu)解的方法，使得條紋噪聲與高頻分量差異最小。從而將問題轉(zhuǎn)化為最小二乘優(yōu)化問題，得到成本函數(shù)E(A)：

Ej(A)=(Nj(i)-Hj(i))2，

(10)

其中：Ej(A) 為第j列成本函數(shù)，Nj(i)為條紋噪聲，以采集到的紅外圖像作為輸入，根據(jù)式(9)得出，k為多項式階數(shù)，Hj(i)為雙邊濾波分離出的高頻分量。

使成本函數(shù)Ej(A)最小，令其為零，再將式(9)帶入可得：

(11)

將式(11)中的I(i)k寫成系數(shù)矩陣P：

(12)

其中：R為圖像行數(shù)，則系數(shù)向量A=(a0，a1，a2，…，ak)的求解公式為：

(13)

其中：H為經(jīng)雙邊濾波處理后得到的高頻分量。Pj為j列的Ik～I以及常量1組成的(k+1)×R矩陣，R為圖像總行數(shù)。Hj為圖像高頻部分的j列。

4 實驗結(jié)果與分析

4.1 實驗設(shè)備

本文算法在Matrix MCX-6300嵌入式工控機(CPU：8核酷睿i7，主頻：2.30 G，內(nèi)存4 G)上運行，使用了Matrox solios eCL/XCL-B圖像采集卡和Matrox vio圖像處理卡作為圖像輸入輸出卡。以Visual Studio 2010為編譯環(huán)境，運用MFC搭建了系統(tǒng)框架，算法代碼調(diào)用了OpenCV 2.4.9庫。紅外熱像儀采用分辨率為320×256的14 bit圖像，輸出幀頻有三檔，分別為50，100，200 Hz 可變。

4.2 參數(shù)設(shè)置

本文可調(diào)節(jié)的參數(shù)主要有3個，其中兩個為雙邊濾波中空域及值域的相似參數(shù)，經(jīng)試驗測試知，空域的相似參數(shù)對算法影響較小，而值域的相似參數(shù)對條紋大小較為敏感，本文試驗取空域相似參數(shù)為10，值域參數(shù)為3。

k值的設(shè)置方法如下，在k=0時，公式(9)變?yōu)镹j(i)=aj0，可以看到誤差與輸入圖像邊緣無關(guān)，因此校正并不能使圖像在有目標前景時得到滿意的結(jié)果。而當k>0時，公式(9)中N為I的函數(shù)，誤差隨輸入圖像的不同而變化，使得在不同場景下具有不同的校正值。條紋噪聲的計算過程如圖2所示，實驗采用320×256的14位紅外圖像，從圖中可以看出，當k取較大值時，條紋噪聲變化不明顯，而處理時間卻有所增加。綜合考慮運行效率及校正效果，本文采用k=1。

k=0，t=17.51 ms

k=1，t=18.39 ms

k=2，t=20.21 ms

k=3，t=21.70 ms

4.3 算法對比及質(zhì)量評價

圖像質(zhì)量的評價實驗在分辨率為320×256的14 bit 圖像上進行，本文選取MIRE算法，雙邊濾波算法作為對比算法，從參照處理時間(t)、信噪比(SNR)、平方差誤差(MSE)、圖像粗糙度(ρ)及國家標準的非均勻性響應率(UR)幾方面進行分析。

表1 幾種算法性能對比

算法質(zhì)量評價結(jié)果如表1所示，可見，本文算法的運算時間遠小于MIRE算法，而信噪比略高于MIRE算法；通過比較MSE可知，本文算法對原圖的改動較大，僅次于雙邊濾波方法，但是雙邊濾波方法對原圖像的模糊作用在圖像粗糙程度這一指標中有所體現(xiàn)；ρ表示圖像的粗糙程度，即高頻信息。分析得到以下結(jié)論：ρ越大，算法保留的邊緣越好，相較其他算法，本文算法的邊緣保留效果更好。本文的非均勻響應率較原圖降低了3%。

幾種算法對3種圖像的處理效果如圖3～5所示。從圖3可以看出，本文算法的視覺效果與MIRE相當，而運算時間遠低于MIRE算法。對于圖4的無條紋噪聲圖像，本文算法對原圖的失真沒有過增強，對于圖5的存在條紋邊緣干擾的圖像，本文算法并沒有出現(xiàn)明顯的邊緣去除現(xiàn)象。

(a)原始圖像

(b)MIRE算法結(jié)果

(c)雙邊濾波算法結(jié)果

(d)本文算法結(jié)果

Fig.3 Performance comparison of different algorithms in sea with complex background

(a)原始圖像

(b)MIRE算法結(jié)果

(c)雙邊濾波算法結(jié)果

(d)本文算法結(jié)果

(a)原始圖像

(b)MIRE算法結(jié)果

(c)雙邊濾波算法結(jié)果

(d)本文算法結(jié)果

4.4 多線程優(yōu)化策略

多線程技術(shù)是使計算機可以在同一時間處理不同任務，可以提高算法的運算速度。鑒于實際應用需求，本程序啟用了4個線程：圖像輸入、圖像輸出、豎條紋噪聲估計和圖像校正。

圖像輸入使用雙內(nèi)存的乒乓操作來捕獲圖像，可與其他線程并行運行；條紋噪聲的估計以一個輔助線程做運算。由于豎條紋噪聲的固定特點，可以近似認為其在短時間內(nèi)變化極小，則對于圖片序列中鄰近幀的校正可以用同一條紋噪聲圖片重復進行處理。該線程循環(huán)地復制當前最新圖片，并作為輸入，計算條紋噪聲后再存儲圖片；圖像校正線程時，則使用已估算出的條紋噪聲圖片，根據(jù)公式(5)進行校正操作。最后，圖像輸出線程經(jīng)過自適應直方圖拉伸算法，將校正后的14位圖像實時轉(zhuǎn)化成8 bit圖像并顯示。這種策略的原理性缺陷在于近似認為條紋噪聲在短時間內(nèi)是不變的。若輔助線程估計噪聲耗時過多，則會導致條紋噪聲更新過慢，校正不準。

以幀頻為200 Hz的相機為例，假設(shè)算法噪聲估計時間為20 ms，輔助線程每4幀估計出一張噪聲圖像，隨后根據(jù)噪聲圖像和當前采集到的圖像進行校正，并即時對最新的圖像進行新的噪聲估計，以更新噪聲圖像。線程工作安排如圖6所示，其中Fx表示第x幀圖像，M1，M2表示乒乓操作的緩沖區(qū)，橫線“—”表示處理中，Nx表示圖像根據(jù)第x幀估計出的噪聲，F(xiàn)xNy表示第x幀圖像根據(jù)第y幀估計噪聲進行校正的結(jié)果，Dx表示轉(zhuǎn)換后的顯示結(jié)果。圖6中噪聲估計線程在復制輸入信息后便在后臺單獨運算，實測系統(tǒng)從幀信號開始到顯示圖像完成的周期為1.5 ms 左右，其中：采集圖像與噪聲估計都是異步處理，在時間上與圖像顯示是并行的，采集時間和估計時間均不包含在內(nèi)，而有效占用時間包括校正圖像耗時(≈1.15 ms)和輸出轉(zhuǎn)換圖像耗時(≈0.35 ms)兩部分。

本文算法對實驗用條紋噪聲圖像的估計處理時間小于20 ms/frame，在200FPS的成像系統(tǒng)下的更新延遲僅為4 frame。而MIRE算法是以列為單位對圖像進行中值直方圖均衡化處理的，故算法運算量大，耗時達1 000 ms，導致其更新的條紋噪聲與當前幀差異過大，容易引起校正錯誤。

圖6 多線程工作安排

圖7 算法多線程處理效果

實際間隔幀數(shù)以輔助線程的計算效率為準，結(jié)果表明圖像顯示周期約為1.5 ms，經(jīng)多線程優(yōu)化后，算法的處理時間較單線程提升了近10倍，而且處理結(jié)果與其并無明顯差異。與經(jīng)緯儀對接的實際處理效果如圖7所示。

5 結(jié) 論

針對紅外相機的非均勻性對成像影響較大這一問題，本文研究了單幀的實時校正算法，從根源上避免了基于場景算法的殘影現(xiàn)象，并結(jié)合噪聲特點對算法進行優(yōu)化。最后在分辨率為320×256的14位紅外圖像上進行了性能驗證實驗，并將MIRE、雙邊濾波算法作為對比算法，進行處理效果及時間對比實驗。實驗結(jié)果證明：算法在分辨率為320×256的14位紅外圖像上的運算時間為1.5 ms/frame，在200 FPS的成像系統(tǒng)上能夠達到實時處理。較好地保持了圖像邊緣，非均勻響應率較原圖像降低了3%。

[1] 劉永進，朱紅，趙亦工. 基于紅外焦平面陣列讀出結(jié)構(gòu)的非均勻校正算法[J]. 光學 精密工程，2008，16(1)：128-133. LIU Y J， ZHU H， ZHAO Y G. Nonuniformity correction algorithm based on infrared focal plane array readout architecture[J].Opt.PrecisionEng.， 2008， 16(1)： 128-133.(in Chinese)

[2] 曹揚，金偉其，劉崇亮，等. 紅外焦平面陣列的自適應非均勻性校正及硬件實現(xiàn)[J]. 光學 精密工程，2011，19(12)：2985-2991. CAO Y， JIN W Q， LIU CH L，etal.. Adaptive nonuniformity correction and hardware implementation of IRFPA [J].Opt.PrecisionEng.， 2011， 19(12)： 2985-2991.(in Chinese)

[3] 代少升，張新立. 利用單片可編程系統(tǒng)對紅外焦平面陣列進行實時非均勻性校正[J]. 光學 精密工程，2009，17(4)：854-858. DAI SH SH， ZHANG X L. Real-time nonuniformity correction of infrared focal plane arrays by system on programmable chip [J].Opt.PrecisionEng.， 2009，17(4)：854-858.(in Chinese)

[4] 張紅輝，羅海波，余新榮，等. 多點標定的自適應非均勻性校正方法[J]. 紅外與激光工程，2014，43(11)：3651-3654. ZHANG H H，LUO H B，YU X R，etal.. Adaptive non-uniformity correction algorithm based on multi-point correction [J].InfraredandLaserEngineering， 2014，43(11)：3651-3654.(in Chinese)

[5] 代少升，張?zhí)祢U. 基于非線性響應的紅外焦平面陣列非均勻性自適應校正算法[J]. 光學 精密工程，2008，16(4)：714-719. DAI SH SH， ZHANG T Q. Adaptive nonuniformity correction algorithm for IRFPA with nonlinear response [J].Opt.PrecisionEng.， 2008，16(4)：714-719.(in Chinese)

[6] NING L，JUN X.Inter frame phase-correlated registration scene-based nonuniformity correction technology [J].InfraredPhysics&Technology， 2015，69： 198-205.

[7] PABLO M， GUILLWEMO M， SERGIO T，etal.. Simultaneous digital super-resolution and nonuniformity correction for infrared imaging system [J].AppliedOptics， 2015，54(21)： 6508-6515.

[8] FAN F， MA Y， HUANG J，etal.. A combined temporal and spatial deghosting technique in scene based nonuniformity correction [J].InfraredPhysics&Technology， 2015，71： 408-415.

[9] 段程鵬，劉偉，謝慶勝，等.改進的神經(jīng)網(wǎng)絡(luò)非均勻校正算法的研究與實現(xiàn)[J].光學學報，2015，35(7)：0710001. DUAN CH P， LIU W， XIE Q SH，etal.. Research and implementation of improved neural network non-uniformity correction algorithm [J].ActaOpticaSinica， 2015， 35(7)：0710001.(in Chinese)

[10] TORRES J，INFANTE S O.Wavelet analysis for the elimination of striping noise in satellite images[J].Opt.Eng.，2001，40(7)：1309-1314.

[11] 潘科辰，顧國華，陳錢，等.一種改進的全變分模型校正紅外焦平面陣列條紋非均勻性[J].紅外技術(shù)，2016，38(2)：138-143. PAN K CH，GU G H，CHEN Q，etal.. An improved total variation model for correcting the strip nonuniformity in IRFPA image [J].InfraredTechnology， 2016， 38(2)：138-143.(in Chinese)

[12] 王書朋，高騰.基于雙邊濾波器的紅外圖像條紋噪聲消除算法[J].紅外技術(shù)，2014，36(9)：728-731. WANG SH P，GAO T. Destriping method for infrared image based on bilateral filter [J].InfraredTechnology， 2014， 36(9)：728-731.(in Chinese)

[13] 張磊，王岳環(huán)，宋瓊.紅外焦平面條帶狀非均勻校正硬件實現(xiàn)方法[J].紅外與激光工程，2015，44(10)：2886-2890. ZHANG L， WANG Y H， SONG Q. Hardware implementation of stripe non-uniformity correction method of infrared focal plane arrays [J].InfraredandLaserEngineering， 2015， 44(10)：2886-2890.(in Chinese)

[14] TENDERO Y，LANDEAU S，GILLES J.Non-uniformity correction of infrared images by midway equalization [J].IPOL， 2012：doi：10.5201/ipol.2012.

李賡飛(1989-)，男，吉林通化人，博士研究生，2012年于吉林大學獲得學士學位，主要研究方向為圖像和視頻信息處理。E-mail：killcolours@126.com.

導師簡介：

李桂菊(1964-)，女，吉林吉林人，研究員，博士生導師，1988年于中國科學院長春光學精密機械與物理研究所獲得碩士學位，研究方向為圖像和視頻信息處理、DSP開發(fā)與應用。E-mail： lgjciom666@126.com

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

Real-time non-uniformity correction of infrared imaging system

LI Geng-fei*， LI Gui-ju， HAN Guang-liang， LIU Pei-xun，JIANG Shan

(Changchun Institute of Optics， Fine Mechanics and Physics，ChineseAcademyofSciences，Changchun130033，China)

A real-time Non-Uniformity Correction(NUC) algorithm for a single-frame infrared image was proposed based on multi-thread optimization to correct the strip-type Fixed Pattern Noise(FPN) generated by non-uniformity response from an infrared camera. Firstly， the plus-relationship of the output value and the noise was determined by NUC algorithm based on total variation theory and a mathematical model was established to describe the strip noise. The image was split into a high-frequency component and a low-frequency component by filters， and the strip noise was fitted from the high-frequency component by the mathematical model. Finally， the real infrared image was restored by the plus-relationship obtained from the total variation theory. To improve the computing speed， the proposed algorithm was optimized with the multi-thread. The image quality evaluation of the proposed algorithm was carried out and its performance was compared with the Midway Infrared Equalization(MIRE) and the bilateral filter methods. Experimental results indicate that the proposed algorithm has reduced the response non-uniformity by 3% as compaired with original image according to national standard GB/T 17444-1998， and its efficiency is no difference with that of the MIRE. Moreover， system running time for an 14 bit infrared image with a resolution of 320×256 is 1.5 ms/frame， which meets the engineering requirements.

infrared image； imaging system； strip-type noise； non-uniformity correction； real-time processing

2016-05-13；

2016-06-20.

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

1004-924X(2016)11-2841-07

TN216

A

10.3788/OPE.20162411.2841

*Correspondingauthor，E-mail：killcolours@126.com