陳 貞，邢笑雪

（1.莆田學(xué)院信息工程學(xué)院，福建莆田351100；

2.長春大學(xué)電子信息工程學(xué)院，吉林長春130022）

基于非下采樣Shearlet變換與壓縮感知的圖像融合

陳 貞1?，邢笑雪2

（1.莆田學(xué)院信息工程學(xué)院，福建莆田351100；

2.長春大學(xué)電子信息工程學(xué)院，吉林長春130022）

針對非下采樣剪切波變換（NSST）分解后圖像的高頻系數(shù)數(shù)據(jù)量較大且具有較大稀疏性的問題，本文提出一種基于NSST和壓縮感知（CS）的圖像融合算法.算法首先采用NSST對源圖像進(jìn)行分解；其次利用CS算法將NSST分解后的圖像的高頻系數(shù)進(jìn)行壓縮、融合及重構(gòu)；然后利用“局部區(qū)域能量和局部區(qū)域方差”聯(lián)合指導(dǎo)待融合圖像的低頻系數(shù)的融合；最后利用NSST逆變換重構(gòu)融合圖像.由于只需要對高頻系數(shù)的壓縮值進(jìn)行融合，因此算法可以在不影響圖像融合效果的同時加快代碼的運(yùn)行速度.仿真實驗表明，該算法不需要原圖像的先驗知識就可以完成圖像的融合，當(dāng)圖像的尺寸較大時，該算法犧牲了微小的融合圖像質(zhì)量，但卻可以顯著提高算法的運(yùn)行速度，減小代碼的時間代價，降低對硬件系統(tǒng)的要求.該算法對于融合系統(tǒng)的實時性要求提供了一種思路，具有較大的應(yīng)用價值.

圖像融合；非下采樣剪切波變換；壓縮感知；局部區(qū)域能量；局部區(qū)域方差

1　引 言

多源圖像融合領(lǐng)域，融合算法的時間代價不僅取決于算法的好壞，很大部分還取決于待處理圖像的數(shù)據(jù)量的大小.如何高效地表示和分析圖像是圖像融合領(lǐng)域的一項核心問題.圖像表示方法的有效性直接決定了融合圖像的質(zhì)量.非下采樣剪切波變換（NSST）是目前最為先進(jìn)的一種多尺度多方向的分析方法［1-2］，它在各個方向和各個尺度上對圖像實現(xiàn)真正的稀疏表示.基于NSST的圖像融合算法能取得較好的融合效果，但是，圖像進(jìn)行多尺度分解后得到的高頻系數(shù)將隨著分解層數(shù)和方向數(shù)的增加而增加，對這些系數(shù)進(jìn)行處理時將大大加劇硬件和軟件的消耗.

壓縮感知（Compressed Sensing，CS），也稱壓縮傳感或壓縮采樣理論，2006年由Candès和Donoho提出［3-4］，其在信號采樣的同時完成信息的壓縮.將CS應(yīng)用到圖像融合可以減少圖像的采樣率、存儲量和計算復(fù)雜度，有效的降低對硬件的要求.因此，近年來，很多學(xué)者將目光轉(zhuǎn)向CS域圖像融合算法的研究［5］；2009年，Xiaoyan Luo等人在CS的基礎(chǔ)上提出了對測量值基于熵的加權(quán)平均的融合規(guī)則［6］；2010年，Juanjuan Han等對文獻(xiàn)［5］中的方法進(jìn)行了改進(jìn)，并與文獻(xiàn)［6］中的方法進(jìn)行比較，得到了更為優(yōu)越的結(jié)果［7］；2012年，Guidong Liu提出了基于CS的超聲波圖像融合算法［8］；2013年，Shutao Li提出了一種基于字典學(xué)習(xí)的稀疏表示的遙感圖像融合算法［9］.

如果能夠?qū)D像數(shù)據(jù)進(jìn)行壓縮感知并進(jìn)行無失真的重構(gòu)，即將CS應(yīng)用到圖像融合領(lǐng)域，既對融合的效果沒有影響，同時又可以加快算法的執(zhí)行效率.本文在綜合分析了CS的基礎(chǔ)上，提出一種基于NSST和CS的圖像融合算法.該算法具有以下特點:一方面，NSST可以完成圖像的多尺度分解和方向濾波，提取出圖像的低頻和高頻信息；另一方面，針對多個待處理的方向子帶，CS可以完成對其的壓縮、感知、融合及重構(gòu)，在較低采樣率下得到較好的融合效果.

2　非下采樣剪切波變換

當(dāng)維數(shù)n=2時，具有合成膨脹的仿射系統(tǒng)AAB（ψ）形式如下:

其中:y∈L2（R2），A，B是2×2可逆矩陣，|det B|=1.如果A AB（ψ）具有如下形式的Parseval框架（也稱緊框架），則這個系統(tǒng)的元素稱為合成小波［10，11］，即對任意的f∈L2（R2），有

其中，矩陣Aj關(guān)聯(lián)尺度變換，Bl關(guān)聯(lián)幾何變換，例如旋轉(zhuǎn)和剪切變換操作等.

為剪切波系統(tǒng)，ψast（x）即為剪切波.

NSST的實現(xiàn)主要分為兩個步驟［12，13］:

（1）非下采樣的多尺度剖分.非下采樣金字塔分解（Non-subsampled Pyramid，NSP）采用二通道非下采樣濾波器組使NSST具備了多尺度性，源圖像經(jīng)一層NSP后將會得到圖像的低頻系數(shù)fa1和高頻系數(shù)fd1，之后每層的NSP分解都在上層分解得到的低頻分量上進(jìn)行迭代以獲取圖像的奇異點.由于在NSST過程中不存在下采樣，因此，圖像f經(jīng)層NSP分解后，最終可得k+1個與f大小相同的子帶圖像.這k+1個圖像包括1個低通圖像和k個帶通子帶圖像.

（2）方向局部化.NSST利用剪切波濾波器實現(xiàn)高頻圖像的方向局部化.具體實現(xiàn)過程為:

1）將偽極化坐標(biāo)映射到笛卡爾坐標(biāo)；

2）利用”Meyer”小波構(gòu)造窗函數(shù)，生成剪切波濾波器；

3）將k個帶通子帶圖像和”Meyer”窗函數(shù)進(jìn)行卷積操作，然后獲得方向子帶圖像.

綜上所述，在固定分辨率水平j(luò)的剪切波變換的框圖如圖1所示.

圖1　NSST的實現(xiàn)流程圖Fig.1　Flow chart of NSST

3　壓縮感知的基本理論

CS是一種全新的、能夠在采樣的同時實現(xiàn)壓縮的理論模型，CS的壓縮、采樣過程見圖2.壓縮感知理論框架主要包括三個重要的內(nèi)容:（1）信號的稀疏表示；（2）觀測矩陣的設(shè)計；（3）信號的重構(gòu).

圖2　CS模型Fig.2　CS model

實現(xiàn)步驟如下:

（1）首先，如果信號X∈RN在某個正交基或緊框架Ψ上是稀疏的或者是可壓縮的，求出變換系數(shù)Θ=ΨTX，Θ是X的等價或逼近的稀疏表示；

假設(shè)一維離散信號X的列向量為X=［X1，X2，…，XN］T，則任意的一組N維向量可以用一個N×1維的向量基｛Ψi｝Ni=1表示，假設(shè)這個向量基正交，則通過矩陣Ψ=［Ψ1，Ψ2，…，ΨN］，信號X可以表示成:

其中，θi=＜ΨiT，X＞，Θ和X是N×1維矩陣，Ψ為N×N維矩陣.

（2）其次，設(shè)計一個平穩(wěn)的、與變換基Ψ不相關(guān)的M×N維的觀測矩陣Φ；

用一個大小為M×N（M＜N）的觀測矩陣Φ=｛φ1，φ2，…，φM｝對稀疏信號Θ進(jìn)行線性投影，計算｛φj｝jM=1與各向量Θ之間的內(nèi)積，經(jīng)過M次觀測，得到線性觀測值yj=＜Θ，φjT＞（j=1，2，…，M），記M×1（M≤N）維的觀測集合為Y=（y1，y2，…，yM）.這個過程也可以看做是信號在觀測矩陣上從高維向低維投影的過程，則Y滿足:

其中:Φ為M×N的矩陣，Y為M×1的矩陣，即矩陣X的觀測值.

（3）最后，利用0-范數(shù)意義下的優(yōu)化問題求解X的精確或近似逼近.

4　基于NSST與CS的圖像融合算法

4.1融合框架

利用NSST對待融合圖像進(jìn)行分解后，由于低頻系數(shù)不具有稀疏性，因此基于NSST和CS的圖像融合算法只對高頻系數(shù)進(jìn)行CS處理，將低頻系數(shù)和高頻系數(shù)分開進(jìn)行處理可以得到更好的融合效果.圖3為NSST和CS相結(jié)合的圖像融合方法示意圖.

圖3　基于NSST和CS的圖像融合框架Fig.3　Image fusion frame based on NSST and CS

設(shè)待融合圖像為A和B，融合后圖像為F.基于NSST和CS的圖像融合的具體操作步驟如下:

Step 1:將待融合的圖像A和B進(jìn)行配準(zhǔn)，本文算法中A和B均已完成配準(zhǔn)；

Step 2:對已配準(zhǔn)的圖像A和B進(jìn)行NSST分解，得到其子帶系數(shù)｛AL，AHj，k｝和｛BL，BHj，k｝.其中，AL、BL表示圖像的低頻系數(shù)，AHj，k、BHj，k表示圖像在第j層第k個高頻子帶系數(shù).分解層數(shù)設(shè)為3層，即j∈｛1，2，3｝，其中，j=1，k=4；j=2，k=4；j=3，k=8；

Step 3:采用“局部區(qū)域方差和局部區(qū)域能量”聯(lián)合指導(dǎo)的策略對Step 2中AL和BL加權(quán)處理，進(jìn)而得到FL；

Step 7:對｛FL，F(xiàn)Hj，k｝進(jìn)行NSST逆變換，得到最終的融合圖像F.

4.2低頻系數(shù)的融合規(guī)則

低頻系數(shù)代表輪廓信息，圖像的低頻分量集中了圖像的主要能量，反映圖像的近似和平均特性，是源圖像的逼近，不具有稀疏特點，因此不能采用CS對其進(jìn)行處理.局部區(qū)域方差表明了該區(qū)域灰度變化的劇烈程度，在一定程度上也反映了該區(qū)域圖像的清晰程度.為此，本文將局部區(qū)域能量和局部區(qū)域方差結(jié)合對待融合圖像的低頻系數(shù)進(jìn)行融合.

圖像A和圖像B的低頻系數(shù)AL（r，c）和BL（r，c）的局部區(qū)域方差為VA（r，c）和VB（r，c），區(qū)域能量EA（r，c）和EB（r，c），計算公式如式（7）～式（9）所示.

其中:Ω（r，c）為以像素點（r，c）為中心的m×n大小的鄰域窗口，一般為3×3或5×5，本文中選取3×3的鄰域窗口.fPL（m，n）為圖像PL在（m，n）處的像素值，w（m，n）為（m，n）處對應(yīng)的窗函數(shù)，EPL（r，c）為像素（r，c）的局部能量.

融合圖像F的低頻系數(shù)FL（r，c）采用加權(quán)平均算子來進(jìn)行融合，計算公式如式（10）所示.

其中:ωA（r，c）和ωB（r，c）滿足式（11）.

4.3高頻系數(shù)的融合規(guī)則

由于圖像的高頻系數(shù)反映的是圖像的細(xì)節(jié)信息，如直線、曲線、邊緣、紋理等信息，系數(shù)的絕對值越大，表示該處亮度變化越劇烈，因此在進(jìn)行圖像融合時，應(yīng)盡量選擇源圖像中清晰的部分來獲得融合系數(shù).NSST分解得到的高頻波段的系數(shù)大部分在零值左右波動，只有小部分的高頻系數(shù)具有較大的幅值，所以依據(jù)CS理論，可以認(rèn)為高頻系數(shù)具有較大的稀疏性.為此，本文選用CS算法對稀疏度較大的高頻系數(shù)進(jìn)行壓縮，這樣可以加快算法運(yùn)行速度，降低時間損耗.

其中，R是n×m的一個矩形區(qū)域，nR=n×m，一般為3×3或5×5，本文中選取3×3的鄰域窗口.j=1，2，3，…代表分解尺度，k=｛k1，k2，…，kj｝表示在每一個分解層上對應(yīng)的方向個數(shù).

5　仿真實驗和結(jié)果分析

為了驗證本文提出的基于NSST和CS的融合算法的有效性和正確性，本文選用紅外與可見光圖像圖4所示、多聚焦圖像圖6所示和醫(yī)學(xué)圖像圖8所示來進(jìn)行仿真實驗.實驗平臺為Matlab2011（a），計算機(jī)主機(jī)配置為:Intel（R）處理器，CPU主頻3.4 GHz，內(nèi)存8 GB.實驗對比了本文提出的算法和其他4種算法在相同的融合規(guī)則下的時間代價及融合效果.其他4種算法分別為:算法（1）:基于NSST的圖像融合算法；算法（2）:基于CS的圖像融合算法；算法（3）:文獻(xiàn)［14］中的基于DWT和CS的圖像融合算法［14］；算法（4）:文獻(xiàn)［15］中的基于NSCT和CS的圖像融合算法［15］.

由于不同的采樣率M/N下圖像的重構(gòu)效果不同，采樣率M/N越高，圖像的重構(gòu)效果越好，但是需要的重構(gòu)時間越長.本文對不同采樣率M/N下（M/N設(shè)置為0.1～1，M以5遞增）的6幅實驗圖像的PSNR、MI、CORR和重構(gòu)時間T進(jìn)行了對比，當(dāng)采樣率M/N為0.3時，在不太影響圖像的重構(gòu)效果的同時所需的重構(gòu)時間也較小.因此，在本文中，CS算法中的采樣率均選擇為0.3.

圖4　“Bristol Queen's Road”源圖像Fig.4　Source images of Bristol Queen's Road

參數(shù)設(shè)置方面，算法（1）中，NSST設(shè)置為3層尺度分解，剪切方向數(shù)為｛4，8，8｝；算法（2）中，CS采用小波稀疏基進(jìn)行稀疏，高斯隨機(jī)矩陣進(jìn)行觀測，OMP算法進(jìn)行重構(gòu)，M/N設(shè)置為0.3；算法（3）中，小波分解層設(shè)定為3層，采用“sym8”小波函數(shù)，CS的參數(shù)設(shè)定與本文算法（2）一致；算法（4）中，NSCT分解層數(shù)設(shè)為3層，尺度分解濾波器采用“maxflat”濾波器，方向分解濾波器選用“dmaxflat”濾波器，CS的參數(shù)設(shè)定與算法（2）一致；算法（5）即本文算法中，NSST設(shè)置與算法（1）一致，CS參數(shù)設(shè)定與算法（2）一致.

對于所有的算法，融合規(guī)則均采用本文提出的融合方法.實驗選用熵值（E）、平均梯度（AG）、標(biāo)準(zhǔn)差（SD）、邊緣信息保留量（QAB/F）、空間頻率（SF）和算法耗時（T）6個參數(shù)來綜合評價各種方法融合圖像的質(zhì)量.

（1）紅外和可見光源圖像

圖5　“Bristol Queen's Road”源圖像融合結(jié)果Fig.5　Fusion results of source images of Bristol Queen's Road

從圖5可以看出，5種算法均可以實現(xiàn)“Bristol Queen's Road”圖像的融合.從視覺效果來看，NSST算法可以得到最好的視覺效果，道路的邊緣、廣告牌上的細(xì)節(jié)信息均得到了很好的體現(xiàn)；CS算法得到的視覺效果是最差的，分析其原因主要是因為實際圖像的稀疏性較差導(dǎo)致重構(gòu)效果較差；NSST+CS算法中首先對圖像進(jìn)行NSST分解，只對其稀疏的高頻系數(shù)進(jìn)行觀測重構(gòu)，其融合效果僅次于NSST算法且優(yōu)于DWT+CS和NSCT+CS.從客觀參數(shù)上來看，NSST+CS算法的耗時T是所有算法中耗時最低的，其SD、E、AG、QAB/F和SF僅小于NSST，但高于CS算法、DWT+CS算法和NSCT+CS算法.

（2）多聚焦源圖像

對圖6和圖7的融合結(jié)果進(jìn)行分析，可以發(fā)現(xiàn)NSST+CS算法可以將多聚焦圖像進(jìn)行較好的融合，源圖像的邊緣和細(xì)節(jié)信息可以被較好的傳遞到融合圖像中去.對表1分析表明，NSST+ CS算法得到的客觀參數(shù)值較大，且算法的時間消耗最小.主觀效果評價和客觀參數(shù)度量均論證了NSST+CS算法是一種切實可行的融合算法.

圖6　“Pepsi”多聚焦源圖像Fig.6　Multi-focus source images of Pepsi

圖7　“Pepsi”多聚焦源圖像融合結(jié)果Fig.7　Fusion results of Multi-focus source images of Pepsi

圖8　“腦血管”源圖像Fig.8　Source images of Cerebrovascular

（3）醫(yī)學(xué)源圖像

圖8和圖9的融合結(jié)果表明，NSST+CS的融合效果弱于NSST，主要是因為CS對高頻系數(shù)進(jìn)行了70%的壓縮，這就導(dǎo)致在重構(gòu)時圖像的邊緣和細(xì)節(jié)信息沒有NSST算法那么豐富；同時，由于NSST+CS算法采用了NSST多尺度分解工具，圖像的細(xì)節(jié)和邊緣可以被很好的提取，因此其融合效果要優(yōu)于DWT+CS和NSCT+CS的算法.從表1可以看出，由于醫(yī)學(xué)圖像的稀疏性較好，所以采用NSST+CS進(jìn)行圖像融合時，不僅可以得到較大的客觀參數(shù)值，并且算法運(yùn)行的時間較快，可以較好的滿足硬件系統(tǒng)的實時性的要求.

圖9　“腦血管”源圖像融合結(jié)果Fig.9　Fusion results of source images of Cerebrovascular

6　結(jié) 論

針對NSST分解后圖像的高頻系數(shù)數(shù)據(jù)量較大且具有較大稀疏性的問題，本文提出一種基于NSST和CS的圖像融合算法.算法的核心思想是首先采用NSST對源圖像進(jìn)行分解；其次利用CS算法將NSST分解后的圖像的高頻子帶系數(shù)進(jìn)行壓縮、融合及重構(gòu)，其中，高頻觀測值的融合采用“局部區(qū)域能量和全局梯度”聯(lián)合指導(dǎo)的策略；然后利用“局部區(qū)域能量和局部區(qū)域方差”聯(lián)合指導(dǎo)低頻分量的融合；最后利用NSST逆變換重構(gòu)融合圖像.

由于只需要對高頻系數(shù)的壓縮值進(jìn)行融合，因此算法可以在不影響圖像的融合效果時又可以加快代碼的運(yùn)行速度.本文分別使用紅外及可見光圖像、醫(yī)學(xué)圖像和多聚焦圖像對算法進(jìn)行了驗證，仿真實驗表明，該算法不需要原圖像的先驗知識就可以完成圖像的融合，當(dāng)圖像的尺寸較大時，該算法雖然犧牲了融合圖像質(zhì)量，但是卻可以顯著提高算法的運(yùn)行速度，減小代碼的時間代價，降低對硬件系統(tǒng)的要求.該算法對于融合系統(tǒng)的實時性要求提供了一種思路，具有較大的應(yīng)用價值.

［1］ Guo K，Labate D，Lim W Q.et al.Wavelets with composite dilations and their MRA properties［J］.Applied and Computational Harmonic Analysis，2006，20（2）:202-236.

［2］ Lim W Q.The discrets shearlet transform:A new directional image representation and compactly supported shearlet frames［J］.IEEE Transactions on Image Processing，2010，19（5）:1166-1180.

［3］ Candès E J.Compressive Sampling［C］.Proceedings of the International Congress of Mathematicians，Madrid，Spain:European Mathematical Society，2006.

［4］ Donoho D L.Compressed sensing［J］.IEEE Transactions on Information Theory，2006，52（4）:1280-1306.

［5］ Wan T，Canagarajah N，Achim A.Compressive image fusion［C］.15th IEEE International Conference on Image Processing，San Diego，CA:IEEE，2008:1308-1311.

［6］ Luo X Y，Zhang J，Yang J Y，et al.Image fusion in compressed sensing［C］.2009 16th IEEE International Conference on Image Processing，Cairo:IEEE，2009:2205-2208.

［7］ Han J J，Loffeld O，Hartmann K，et al.Multi image fusion based on compressive sensing［C］.2010 International Conference on Audio Language and Image Processing，Shanghai:IEEE，2010:1463-1469.

［8］ Liu G D，Shen Y.Ultrasonic image fusion using compressed sensing［J］.Electronics Letters，2012，48（19）:1182.

［9］ Li S T，Yin H T，F(xiàn)ang L Y.Remote sensing image fusion via sparse representations over learned dictionaries［J］. IEEE Transactions on Geosciences and Remote Sensing，2013，51（9）:4779-4789.

［10］ Easley G，Labate D，Lim W Q.Sparse directional image representations using the discrete shearlet transform［J］. Applied and Computational Harmonic Analysis，2008，25（1）:25-46.

［11］ Guo K H，Labate D.Optimally sparse multidimensional representation using shearlets［J］.SIAM J.Math. Anal.，2007，39（1）:298-318.

［12］ 王雷，李彬，田聯(lián)房.基于平移不變剪切波變換的醫(yī)學(xué)圖像融合［J］.華南理工大學(xué)學(xué)報（自然科學(xué)版），2011，39（12）:13-19.

Wang L，Li B，Tian L F.Medical image fusion based on shift-invariant shearlet transformation［J］.Journal of South China University of Technology（Natural Science Edition），2011，39（12）:13-19.（in Chinese）

［13］ 高印寒，陳廣秋，劉妍妍.基于圖像質(zhì)量評價參數(shù)的非下采樣剪切波域自適應(yīng)圖像融合［J］.吉林大學(xué)學(xué)報（工學(xué)版），2014，44（1）:225-234.

Gao Y H，Chen G Q，Liu Y Y.Adaptive image fusion based on image quality assessment parameter in NSST system［J］.Journal of Jilin University（Engineering and Technology Edition），2014，44（1）:225-234.（in Chinese）

［14］ 楊粵濤.基于非采樣Contourlet變換的圖像融合［D］.長春:中國科學(xué)院長春光學(xué)精密機(jī)械與物理研究所，2012.

Yang Y T.Research on image fusion based on nonsubsamepled Contourlet transform［D］.Changchun:Changchun Institute of Optic，F(xiàn)ine Mechanics and Physics of Chinese Academy of Science，2012.（in Chinese）

［15］ 沈瑜.一種改進(jìn)小波變換的圖像融合方法［J］.自動化與儀器儀表，2013（5）:202-204.

Shen Y.An improved wavelet transform image fusion algorthm［J］.Automation&Instrumentation，2013（5）: 202-204.（in Chinese）

Image fusion algorithm based on non-subsampled Shearlet transform and compressed sensing

CHEN Zhen1?，XING Xiao-xue2

（1.Information Engineering Department，Putian University，Putian 351100，China；
2.College of Information Engineering，Changchun University，Changchun 130022，China）

After the image decomposition with NSST，the high-frequency coefficients have a large amount of data and greater sparsity.In order to obtain fusion results rapidly，an image fusion algorithm based on Non-subsampled Shearlet Transform（NSST）combined with Compressed Sensing（CS）is presented.Firstly，the source images are decomposed with NSST；secondly，the high-frequency sub-band coefficients of the decomposed images are compressed，fused and reconstructed by CS；then，based on local area variance and local area energy，the low-frequency coefficients was fused；finally，the inverse NSST is used to get the final fused image.Because only the compressed values of the high frequency coefficients are fused，the image fusion effects can't be affected，and the running time of the algorithm can be reduced.In this paper，the multi-focus image，medical image and infrared andvisible images are used to verify the effectiveness of the algorithm.The simulation results indicate that this algorithm can achieve the fusion of the image without prior knowledge of the original image. When the image size is larger，although the fusion image quality is sacrificed，it can significantly improve the speed to reduce the time cost and hardware requirements.The algorithm provides an idea on how to satisfy the real time requirements in the fusion system，which has a great practical value.

image fusion；non-subsampled Shearlet transform；compressed sensing；local area variance；local area energy

TP391.41

A doi:10.3788/YJYXS20153006.1024

1007-2780（2015）06-1024-08

陳貞（1977-），女，福建莆田人，碩士，講師，主要從事圖像處理、粗糙集理論方面的研究.E-mail:chenfjpt@ yeah.net

邢笑雪（1981-），女，山西霍州人，博士，講師，研究方向為數(shù)字圖像處理、模式識別等.

2015-01-05；

2015-02-05.

國家教育部“春暉計劃”項目（No.Z2014138）；吉林省教育廳基金項目（No.20140529）；莆田市科技計劃項目（No.2014G20）

Supported by Chunhui Project of the National Education Department，P.R.（No.Z2014138）；Education Department Foundation of Jilin Province（No.20140529）；Science and Technology Plan of Putian City（No.2014G20）

?通信聯(lián)系人，E-mail:chenfjpt@yeah.net