周亞同,劉志峰,張志偉
(1.河北工業(yè)大學電子信息工程學院,天津300401;2.中國科學院地質與地球物理研究所,北京100029)
形態(tài)分量分析框架下基于DCT與曲波字典組合的地震信號重建
周亞同1,2,劉志峰1,張志偉1
(1.河北工業(yè)大學電子信息工程學院,天津300401;2.中國科學院地質與地球物理研究所,北京100029)
針對單一型數(shù)學變換或字典不能有效刻畫地震信號的形態(tài)特征多樣性這一問題,在形態(tài)分量分析(MCA)框架下,提出了一種基于離散余弦變換(DCT)與曲波字典組合的地震信號重建方法。該方法首先建立MCA框架下的地震信號重建模型,依托模型將信號分解成局部奇異形態(tài)分量以及平滑與線狀形態(tài)分量。然后采用DCT字典表示局部奇異分量,采用曲波字典表示平滑與線狀分量。再以迭代求解方式逐一重建各分量,最后將重建后的分量合并。人工合成地震信號和二維疊前及疊后實際地震信號重建實驗結果表明,該方法能很好完成信號重建,重建精度不僅要高于非抽樣小波變換(UDWT)與曲波字典組合、曲波與曲波字典組合、余弦調制濾波器組與曲波字典組合,而且更高于DCT,UDWT,或曲波等單一型字典。
形態(tài)分量分析;地震信號重建;離散余弦變換;曲波變換;字典
地震信號重建在地震資料處理中具有重要意義[1-5]?,F(xiàn)有地震信號重建方法大致分為以下3類:基于波場算子、基于預測濾波、基于數(shù)學變換的方法。在基于數(shù)學變換這類重建方法中,早年采用的是傅里葉變換[1]和Randon變換[6]等單一型經(jīng)典變換。近年隨著信號的多尺度幾何分析理論的發(fā)展,涌現(xiàn)出一些新的變換,如曲波變換[7]、Seislet變換[8]和Dreamlet變換[9]等并被用于地震信號重建。另外隨著稀疏信號表示及壓縮感知理論的發(fā)展,又有學者嘗試用超完備學習型字典來重建地震信號[10-11]。字典作為由若干個基函數(shù)所組成的超完備冗余框架,與單一型變換相比擁有更強的表達能力。但無論是疊前還是疊后地震剖面,通常都是由多種形態(tài)特征明顯的元素組合而成的。以疊后剖面為例,既包含斷點或尖滅點等局部奇異元素,也包含弱反射區(qū)或連續(xù)同相軸等平滑與線狀元素。上述不同的元素對應著地震信號的不同內(nèi)蘊特征,因此僅靠單一型變換或單一字典未必能有效刻畫這種多樣性的特征。
李海山等曾提出一種基于形態(tài)分量分析(morphological component analysis,MCA)的地震信號重建方法[12]:首先將地震信號分解成多個形態(tài)分量,然后從字典組合中挑選合適的字典分別重建各分量,最后將重建結果合并。這種分而治之的策略受數(shù)字圖像修復啟發(fā)[13],用字典組合來表達地震信號中的不同內(nèi)蘊特征,探索了地震信號重建的新途徑。MCA框架下的信號分解與重建,面臨的一個關鍵問題是字典組合的選擇。字典組合選擇不僅要匹配被表達信號的不同內(nèi)蘊特征,而且要保證信號分解和重建可以數(shù)值化實現(xiàn)。李海山等在重建地震信號時選用的是非抽樣小波變換(UDWT)與曲波字典組合,即用UDWT字典表示局部奇異分量,用曲波字典表示平滑與線狀分量,取得了良好的重建效果。
本文在李海山等工作的基礎上,嘗試在MCA框架下用多種不同的字典組合重建地震信號,并對重建結果進行定量評價,以此為依據(jù)挑選出最優(yōu)的字典組合。合成地震信號、疊前及疊后實際地震信號重建結果表明,本文選用的離散余弦變換(DCT)與曲波字典組合具有最優(yōu)的重建精度。這可能是因為DCT字典依靠局域離散余弦變換,由多個固定大小的子塊構成冗余字典。字典中稀疏的高頻部分包含的是奇異性較強的局部細節(jié)特征信息,適合表示斷點或尖滅點等局部奇異分量。而曲波字典具有很強的方向性,適合處理形態(tài)各異的同相軸紋理,即能對剖面的平滑與線狀分量提供穩(wěn)定、高效和近似最優(yōu)表示。
1.1 地震信號重建問題
二維地震信號重建問題可通過(1)式描述:
Y=MF
(1)
式中:F代表完整的二維原始地震采樣信號;Y代表二維缺失地震采樣信號;M為由元素0和1組成的采樣矩陣,且元素滿足:
(2)
式中:Ω代表地震信號中的缺損待重建區(qū)域。地震信號重建就是由缺失信號Y和采樣矩陣M恢復出完整信號F的過程。
1.2 形態(tài)分量分析(MCA)框架基本原理
MCA是一種基于多種基函數(shù)的稀疏信號分解框架[14],可以看作是經(jīng)典基追蹤算法(BP)[15]和匹配追蹤算法(MP)[16]的結合。MCA利用信號組成分量的形態(tài)多樣性,將信號分解成多個分量的線性組合,不同分量用不同字典稀疏表示。每個分量對應的字典僅能稀疏表示該分量,而對其它分量不能稀疏表示?,F(xiàn)假設二維信號F可表示成K個形態(tài)分量Fk之和:
(3)
式中:αk為分解系數(shù);Φk為過完備字典??梢酝ㄟ^求解(4)式約束規(guī)劃得到系數(shù)αk:
(4)
式中:符號“‖‖1”代表L1范數(shù)。(4)式也可改寫成以下形式:
(5)
式中:λ為給定閾值。上述求解系數(shù){α1,…,αK}也可以轉化為求解信號分量{F1,…,FK}問題:
(6)
1.3 MCA框架下的地震信號重建模型
根據(jù)MCA的基本原理,若在該框架下由缺失地震信號Y和采樣矩陣M重建完整信號F,考慮到(1)式和(3)式,可以通過改寫(5)式來實現(xiàn)重建:
(7)
MCA框架下的地震信號重建可用圖1直觀闡釋:采用分而治之策略,首先將地震信號分解成多個形態(tài)分量Fk(k=1,2,…,K),然后用不同的字典Φk分別表示各分量,最后將重建結果合并。圖1中字典Φk對重建質量有重要影響。理論上可以根據(jù)一些預設準則自適應構造字典。但這樣的構造方法通常比較復雜,在工程中往往根據(jù)經(jīng)驗選用一些現(xiàn)有變換基作為字典。如曲波字典、局部脊波字典、DCT字典、小波字典、小波包字典、Gabor字典等。
圖1 MCA框架下的地震信號重建模型
MCA框架下的地震信號分解與重建,面臨的一個關鍵問題是字典組合{Φ1,…,Φk,…,ΦK}的選擇。本文嘗試使用多種不同的字典組合進行地震信號重建,并對重建結果進行定量評價,以此為依據(jù)挑選出最優(yōu)的字典組合。為計算方便,假設(7)式中形態(tài)分量數(shù)目K=2,所選用的4種字典組合分別為:DCT與曲波字典組合、UDWT與曲波字典組合、曲波與曲波字典組合、余弦調制濾波器組(CMFB)與曲波字典組合。其中,CMFB通過對線性相位低通原型濾波器進行余弦調制來實現(xiàn)。另外為便于對比,重建中還采用了單一型字典如DCT字典、UDWT字典和曲波字典等,此時形態(tài)分量數(shù)目K=1。
根據(jù)MCA框架下的地震信號重建模型,在選定字典組合以后,首先需要構造出各字典,然后基于(7)式計算各字典對應的分解系數(shù),再基于(3)式獲得重建信號。目前大致存在兩種構造過完備字典的方式:第一種是直接基于現(xiàn)有線性變換構造;第二種是基于學習方式,在給定樣本集上學習出能夠稀疏表示樣本的字典[11]。由于第一種方式能借助快速變換和反變換實施數(shù)值計算,而且產(chǎn)生的字典往往高度結構化,因此本文采用第一種方式構造字典。
2.1 采用DCT字典表示局部奇異分量
DCT是數(shù)字信號處理領域常用的一種正交變換,有很強的能量集中特性,算法復雜度也比較低。設某信號塊F(x,y)經(jīng)離散化后共有M道,每道有N個采樣點,對應的二維離散余弦變換定義如下:
(8)
其對應的反變換為:
(9)
DCT字典基于DCT變換構造。對于DCT變換后所獲得的完備字典,采用分數(shù)頻率法將其擴展成為過完備字典,具體做法是將得到的字典在其頻率上做更精細的遍歷和抽樣,從而獲得一個新的過完備字典??紤]到DCT字典自身特點,本文用其表示地震信號的局部奇異分量。
2.2 采用曲波字典表示平滑與線狀分量
曲波變換是一種建立在脊波變換基礎上的廣泛使用的多尺度幾何分析手段[17]。該變換將任意平方可積空間L2(R)中的函數(shù)映射為系數(shù)序列,實質上是通過把曲線無限分割,近似的把每段看作是直線段,然后對每個直線段進行脊波變換。曲波變換具有多尺度、多方向性、各向異性等特點,能有效表示和檢測到信號中包含的多方向的線狀變化特征[18-19],因此十分適合用于表示地震信號中的線狀同相軸。由于曲波變換在時間域內(nèi)對信號的定義域進行正方形塊分割,會造成邊緣處出現(xiàn)不連續(xù)性,因此,Candès等[20]又提出了第二代曲波變換,把頻率域劃分為一組具有不同尺度的楔形區(qū)域,然后在每個區(qū)域上定義一個基函數(shù)。
曲波變換采用基函數(shù)和信號內(nèi)積的形式實現(xiàn)信號稀疏表示,曲波變換可表示如下:
(10)
式中:符號“〈·〉”代表內(nèi)積;cj,l,k為系數(shù)序列;φj,l,k為曲波函數(shù),包含尺度j,方向l和位置k等參數(shù)。曲波變換緊框架通過標準正交基分解,可以表示一系列曲波函數(shù)的加權疊加,在L2范數(shù)意義下,對應的重建公式可以表示為:
(11)
(12)
曲波超完備字典基于曲波緊框架構造。首先對該緊框架中的尺度參數(shù)j,方向參數(shù)l以及位置參數(shù)k進行離散化,然后再平移滑動窗對最小尺度下不同方向的曲波原子進行截取,最終獲得曲波超完備字典。由于曲波緊框架對信號的平滑與線狀部分提供穩(wěn)定、高效和近似最優(yōu)表示,因此基于該框架得到的超完備字典,能夠很好地匹配地震信號中的平滑與線狀分量,能處理形態(tài)各異的同相軸紋理。鑒于曲波超完備字典的上述特性,本文用其表示二維地震信號中的平滑與線狀分量。
2.3 MCA迭代重建的數(shù)值實現(xiàn)
在構造出各個字典后,需要基于(7)式計算各系數(shù)αk(k=1,2,…,K),進而基于(3)式重建出完整的地震信號。直接求解(7)式十分困難,下面以形態(tài)分量數(shù)目K=2為例,分步闡述以迭代求解方式實現(xiàn)重建的全過程。
1) 根據(jù)地震信號缺失情況,確定(1)式中的采樣矩陣M,同時根據(jù)地震信號的形態(tài)特征選擇合適的字典組合{Φ1,Φ2}。
3) 初始化總殘差R(0)=Y,其中Y為(1)式中的已知缺失地震信號。
4) 由硬閾值法得到起始閾值λ0,同時設置迭代停止閾值λstop。
5) 設置最大迭代次數(shù)Imax,在本文中最大迭代次數(shù)一般設置為50次。
8) 計算重建地震信號和原始地震信號間的各種誤差及相似度,以此評價重建質量。
為了檢驗重建方法的效果,特設計人工合成地震信號和二維疊前、疊后實際地震信號重建3個實驗。實驗在CPU 2.53GHz,內(nèi)存2GB,預裝Windows 7旗艦版32位操作系統(tǒng)的個人計算機上進行,軟件編程環(huán)境為MATLAB 7.10.0(R2010a)。實驗中采用如下5個指標定量評價地震信號的重建精度:均方誤差(MSE)、均方根誤差(RMSE)、平均絕對誤差(MAE)、峰值信噪比(PSNR)和結構相似度(SSIM)。
(13)
(14)
(15)
(16)
另外為度量地震信號重建前、后在整體結構形態(tài)方面的差異,定義二者間的結構相似度[21]為:
(17)
3.1 人工合成地震信號的重建過程展示
為清晰展示重建過程,選擇如圖2a所示某簡單楔形合成地震信號,共有96道,單道160個采樣點,道間距為10m,時間采樣率為0.001s。圖2b為隨機缺失48道后的信號(缺失率達50%)?,F(xiàn)采用本文算法對之進行重建,閾值類型設置為硬閾值,形態(tài)分量數(shù)目設為2,迭代起始閾值λ0設為1.3178,停止閾值λstop為10-6。圖2c和圖2d是在MCA框架下,分別采用UDWT與曲波字典組合、DCT與曲波字典組合分別進行重建的結果。從圖2c和圖2d可以看出,兩種字典組合均能很好地完成重建。表1給出了采用兩種字典組合重建時,合成地震信號的6個重建指標。從表1可以看出,DCT與曲波字典組合比UDWT與曲波字典組合的重建性能更好,重建時間也短一些(因為DCT的計算復雜度要比UDWT的計算復雜度低)。
圖2 合成地震信號以及MCA框架下兩種字典組合方式的重建結果
圖3給出了在MCA框架下采用各種字典組合重建時,峰值信噪比(PSNR)和均方根誤差(RMSE)隨迭代次數(shù)的變化曲線。從圖3中可以看出:①在MCA框架下采用字典組合方式,其重建效果比采用單一字典要好;②在4種字典組合方式中,DCT與曲波字典組合的重建效果最好,其次是UDWT與曲波字典組合,再次是曲波與曲波字典組合,CMFB與曲波字典組合的重建效果最差;③MCA框架下的3種單一字典重建方式的效果大致相當。
表1 MCA框架下采用兩種字典組合時的合成地震信號重建質量評價
圖3 各種字典組合重建的峰值信噪比(a)和均方根誤差(b)隨迭代次數(shù)的變化曲線
圖4a給出了缺失地震信號的采樣矩陣,即(1)式中的矩陣M。圖4b給出了采用DCT與曲波字典組合進行迭代重建時,閾值λi隨迭代次數(shù)的變化曲線。圖4c為重建前、后地震信號的采樣值統(tǒng)計直方圖,其中左圖是重建前原始合成地震信號(圖2a)的采樣值統(tǒng)計直方圖,右圖是由DCT與曲波字典組合重建后的地震信號(圖2d)的統(tǒng)計直方圖。對比發(fā)現(xiàn)二圖非常接近,說明DCT與曲波字典組合確實能很好地完成重建。
圖4 缺失地震信號的采樣矩陣(50%道缺失)(a)、迭代重建閾值曲線(b)以及重建前、后地震信號采樣值統(tǒng)計直方圖(c)
3.2 二維疊前實際海洋地震信號重建結果比較
本實驗選取某二維疊前實際海洋地震信號進行重建。圖5a為原始地震信號,共有64道,每道含256個采樣點,道間距為10m,時間采樣率為0.001s。圖5b為隨機缺失19道后的信號。在MCA框架下,采用實驗一中提及的各種字典組合方式對其進行重建,重建結果分別如圖5c至圖5i所示。對于DCT與曲波字典組合,所采用的重建參數(shù)為:閾值類型為硬閾值,形態(tài)分量數(shù)目為2,迭代起始閾值λ0為824.9917,停止閾值λstop為10-6。從圖5 中可以看出,采用各種字典組合都能恢復缺失道,而且能保持與相鄰道的振幅關系。但仔細對比各圖發(fā)現(xiàn),每種字典組合的重建結果均不相同。
圖5 二維疊前實際海洋地震信號以及MCA框架下各種字典組合方式的重建結果
為了定量比較評價各種字典組合的重建效果,逐一計算6個重建指標并列于表2中。從表2可以看出:①在MCA框架下采用字典組合,其重建效果比采用單一字典要好;②在所有單一字典或字典組合中,DCT與曲波字典組合具有最優(yōu)重建精度;③在相同的計算條件下,字典組合所需重建時間要比采用單一字典長。在實際計算中可通過減少迭代次數(shù)、選擇合適的停止閾值等方式來減少計算時間。
3.3 二維疊后實際地震信號重建結果比較
疊后地震信號與疊前地震信號相比,通常同相軸的形態(tài)更加多樣化。而且部分同相軸可能存在彎曲并伴隨斷點等,因此包含的元素類型更多,重建難度更大。選取如圖6a所示某疊后實際地震剖面片段,該剖面片段共128道,單道含256個采樣點,道間距為10m,時間采樣率為0.001s。圖6b顯示的是缺道疊后剖面,在MCA框架下采用不同的字典組合對其進行重建。圖6c和圖6d 分別給出了UDWT與曲波字典組合、DCT與曲波字典組合的重建結果。從圖6可以看出,兩種字典組合均能完成重建,但重建結果仍有差異。
表2 MCA框架下采用不同字典組合時的疊前實際海洋地震信號重建質量評價
表3列出了采用不同字典組合時,疊后實際地震剖面的6個重建指標。無論是從MSE,RMSE或MAE,還是從PSNR或SSIM等指標來看,DCT與曲波字典組合均具有最優(yōu)的重建精度。另外采用字典組合方式的重建效果要比采用單一字典好,但所需重建時間更長。上述結論與疊前地震信號重建一致。
圖6 疊后實際地震信號以及MCA框架下兩種字典組合方式的重建結果
表3 MCA框架下采用不同字典組合時的疊后實際地震剖面重建質量評價
字典組合類型MSERMSEMAEPSNRSSIM重建時間/sDCT與曲波字典組合3.42251.84980.801526.86500.9531104.7915UDWT與曲波字典組合5.47182.33921.119924.82610.9252124.5500曲波與曲波字典組合6.15642.48121.290624.31420.9159131.3001CMBF與曲波字典組合5.67312.49631.256124.58680.9209221.1349單一DCT字典6.67942.58451.331823.96010.909810.9905單一UDWT字典6.29752.50951.415824.21580.908426.6257單一曲波字典6.21202.48231.313024.29260.913524.3287
在MCA框架下,提出了一種基于DCT與曲波字典組合的地震信號重建方法。3個重建實驗從不同角度驗證了該方法的有效性,實驗結果表明,DCT與曲波字典組合的重建精度不僅要高于UDWT與曲波字典組合、曲波與曲波字典組合、CMBF與曲波字典組合,而且要高于DCT,UDWT,或曲波等單一型字典。本文提出的重建方法在借助MCA框架發(fā)掘地震信號的內(nèi)蘊特征、進而對信號進行更精細、靈活刻畫方面提供了啟示。下一步研究工作包括MCA框架下的字典自適應選擇、地震信號分量分解質量的精細評價、迭代重建數(shù)值計算效率提升等。另外地震信號的反假頻重建在實際資料處理中非常重要,下一步研究將考慮在本文的重建方法中引入反假頻環(huán)節(jié)。
[1] Duijndam A J W,Schonewille M A,Hindriks C O H.Reconstruction of band-limited signals,irregularly sampled along one spatial direction[J].Geophysics,1999,64(2):524-538
[2] Zhang Y,Zhang G,Bleistein N.True amplitude wave equation migration arising from true amplitude one-way wave equations[J].Inverse Problems,2003,19(5):1113-1138
[3] Hindriks K,Duijndam A J W.Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates[J].Geophysics,1999,64(1):253-263
[4] Zwartjes P,Gisolf A.Fourier reconstruction of marine-streamer data in four spatial coordinates[J].Geophysics,2006,71(6):171-186
[5] Schonewille M A,Romijn R,Duijndam A J W,et al.A general reconstruction scheme for dominant azimuth 3D seismic data[J].Geophysics,2003,68(6):2092-2105
[6] Kabir M M N,Verschuur D J.Restoration of missing offsets by parabolic Radon transform[J].Geophysical Prospecting,1995,43(3):347-368
[7] 張華,陳小宏.基于jitter采樣和曲波變換的三維地震數(shù)據(jù)重建[J].地球物理學報,2013,56(5):1637-1649 Zhang H,Chen X H.Seismic data reconstruction based on jittered sampling and curvelet transform[J].Chinese Journal of Geophysics,2013,56(5):1637-1649
[8] Fomel S,Liu Y.Seislet transform and seislet frame[J].Geophysics,2010,75(3):25-38
[9] Wu R S,Geng Y,Ye L L.Preliminary study on dreamlet based compressive sensing data recovery[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,3585-3590
[10] 唐剛.基于壓縮感知和稀疏表示的地震數(shù)據(jù)重建與去噪[D].北京:清華大學,2010 Tang G.Seismic data reconstruction and denoising based on compressive sensing and sparse representation[D].Beijing:Tsinghua University,2010
[11] 周亞同,王麗莉,蒲青山.壓縮感知框架下基于K-奇異值分解字典學習的地震數(shù)據(jù)重建[J].石油地球物理勘探,2014,49(4):652-660 Zhou Y T,Wang L L,Pu Q S.Seismic data reconstruction based on K-SVD dictionary leaning under compressive sensing framework[J].Oil Geophysical Prospecting,2014,49(4):652-660
[12] 李海山,吳國忱,印興耀.形態(tài)分量分析在地震數(shù)據(jù)重建中的應用[J].石油地球物理勘探,2012,47(2):236-243 Li H S,Wu G C,Yin X Y.Morphological component analysis in seismic data reconstruction[J].Oil Geophysical Prospecting,2012,47(2):236-243
[13] 張濤,洪文學.基于自適應字典選擇的MCA圖像修復方法[J].光學技術,2010,36(5):672-676 Zhang T,Hong W X.Image inpainting based on MCA featured adaptive dictionary selection[J]. Op-tical Technique,2010,36(5):672-676
[14] Starck J L,Moudden Y,Bobin J,et al.Morphological component analysis[J].Proceedings of SPIE,2005(5914):1-15
[15] Chen S S.Basis pursuit[D].California:Stanford University,1995
[16] Mallat S,Zhang Z.Matching pursuit with time frequency dictionaries[J].IEEE Transactions on Signal Processing,1993,41(12):3397-3415
[17] Starck J L,Candes E J,Donoho D L.The curvelet transform for image denoising[J].IEEE Transactions on Image Processing,2002,11(6):670-684
[18] Candès E J,Donoho D L.Curvelets a surprisingly effective nonadaptive representation for objects with edges[M].Vanderbilt:Vanderbilt University Press,1999:71-89
[19] Candès E J,Demanet L,Donoho D L,et al.Fast discrete curvelet transforms[J].Multiscale Modeling and Simulation,2006,5(3):861-899
[20] Candès E J,Donoho D L.New tight frames of curvelets and optimal representations of objects with smooth singularities[R].California:Department of Statistics,Stanford University,2002
[21] Wang Z,Bovik A C,Sheikh H R,et a1.Image quality assessment:from error visibility to structural similarity[J].IEEE Transactions on Image Processing,2004,13(4):600-612
(編輯:陳 杰)
Seismic signal reconstruction under the morphological component analysis framework combined with discrete cosine transform (DCT) and curvelet dictionary
Zhou Yatong1,2,Liu Zhifeng1,Zhang Zhiwei1
(1.SchoolofElectronicandInformationEngineering,HebeiUniversityofTechnology,Tianjin300401,China;2.InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China)
Aiming at the problem that the mathematic transforms and dictionaries can not effectively depict the morphological features diversity of seismic signals,we propose a seismic signal reconstruction method under the morphological component analysis (MCA) framework combined with discrete cosine transform (DCT) and curvelet dictionary.Firstly the seismic signal reconstruction model is built under the MCA framework.Then the signal is decomposed into local singular and smooth linear component based on the model.Following that,the local singular component is represented by DCT dictionary,and the smooth linear component is represented by curvelet dictionary.We combine two kinds of components together after iterative reconstruction.The experiments on synthetic and real seismic signals illustrates that the proposed method can be used to reconstruct signals very well.The reconstruction precision of the method is not only higher than some dictionary combinations such as UDWT & curvelet dictionary combination,curvelet & curvelet dictionary combination,CMBF & curvelet dictionary combination,but also higher than some single dictionaries such as DCT,UDWT or curvelet etc.
morphological component analysis,seismic signal reconstruction,discrete cosine transform,curvelet transform,dictionary
2014-12-02;改回日期:2015-05-30。
周亞同(1973—),男,博士,教授,主要研究方向為地震信號處理、模式識別與機器學習等。
國家自然科學基金項目(60972106,51475136)、河北省自然科學基金項目(F2013202254)和中國博士后科學基金項目(2014M56053)共同資助。
P631
A
1000-1441(2015)05-0560-09
10.3969/j.issn.1000-1441.2015.05.009