劉 偉,曹思遠(yuǎn),崔 震

(1.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室,北京昌平102249;2.中國石油大學(xué)CNPC物探重點實驗室,北京昌平102249)

基于壓縮感知和TV準(zhǔn)則約束的地震資料去噪

劉 偉1,2,曹思遠(yuǎn)1,2,崔 震1,2

(1.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室,北京昌平102249;2.中國石油大學(xué)CNPC物探重點實驗室,北京昌平102249)

考慮到傳統(tǒng)曲波變換閾值法難以獲得理想的去噪效果,提出一種基于壓縮感知和全變差(Total Variation,TV)準(zhǔn)則約束的曲波變換地震資料去噪方法。該方法將傳統(tǒng)的去噪問題轉(zhuǎn)化成TV準(zhǔn)則約束下的最優(yōu)化問題,通過求取最優(yōu)解重構(gòu)原始信號,從而實現(xiàn)對地震資料的去噪處理。理論模型數(shù)據(jù)測試結(jié)果表明,該方法在地震資料去噪過程中可以有效壓制同相軸邊緣附近產(chǎn)生的不光滑畸變,相比傳統(tǒng)的曲波變換閾值法能夠獲得更好的去噪效果;與F-X反褶積相比,該方法不僅能夠有效壓制地震資料中的隨機噪聲,提高地震資料的信噪比,而且能夠較好地保護(hù)有效信號。某工區(qū)實際地震資料去噪結(jié)果進(jìn)一步驗證了上述方法的有效性。

壓縮感知;TV準(zhǔn)則;曲波變換;隨機噪聲;噪聲衰減

地震信號在傳播過程中不可避免地受到噪聲的影響,使地震資料的品質(zhì)降低,給后續(xù)解釋工作帶來困難；因此,提高地震資料信噪比是地震資料處理的主要任務(wù)之一。在實際生產(chǎn)中,常用的壓制噪聲的方法主要有F-X反褶積[1]、多項式擬合[2]、F-K濾波[3]、奇異值分解[4]和小波變換[5]等。然而這些方法在去噪處理過程中都不同程度地存在局限性,如作為全局變換的傅里葉變換對奇異點或奇異曲線的局部特征刻畫不夠準(zhǔn)確,F-X反褶積使地震數(shù)據(jù)中的所有相干信息都加強,小波變換對于高維信號的邊緣信息表達(dá)不夠準(zhǔn)確等,所有這些都極大地影響了地震資料的去噪效果。

壓縮感知是近幾年發(fā)展起來的一種新的圖像處理方法[6-7],目前已經(jīng)廣泛應(yīng)用于天文學(xué)、核磁共振和模式識別等領(lǐng)域,且得到了良好的發(fā)展。在醫(yī)學(xué)信號檢測、特征提取、故障診斷和圖像融合等方面,壓縮感知也是眾多學(xué)者研究的熱點[8]。該方法突破了Nyquist采樣定理的限制,在一定條件下可以利用欠采樣數(shù)據(jù)重構(gòu)原始數(shù)據(jù);利用地震數(shù)據(jù)在變換域中的稀疏特性,通過設(shè)定稀疏基矩陣和測量矩陣,可以將地震數(shù)據(jù)去噪問題轉(zhuǎn)化成求解最優(yōu)化問題。

用于地震數(shù)據(jù)稀疏表示的方法有很多,常用的有傅里葉變換、離散余弦變換、短時傅里葉變換、小波變換和曲波變換等[9]。傅里葉變換[10]和離散余弦變換[11]都是一種全局變換,不能有效地刻畫信號的局部特征;短時傅里葉變換[12]是傅里葉變換的進(jìn)一步發(fā)展,使用窗函數(shù)作用于信號后再進(jìn)行傅里葉變換,本質(zhì)上是一種單一分辨率分析方法;小波變換[13]繼承并發(fā)展了短時傅里葉變換的局部化思想,可以根據(jù)待分析信號的特征自動調(diào)節(jié)窗口大小,從而實現(xiàn)對信號的多分辨率分析,但不能很好地用于具有曲線狀特征的地震數(shù)據(jù)稀疏表示。近些年興起的曲波變換[14]由不同尺度、不同方向的曲線元素構(gòu)成變換基,在保持小波變換多分辨率特性的同時增加了小波變換的多方向性和各向異性,更加適用于刻畫地震數(shù)據(jù)的特征,是對地震數(shù)據(jù)進(jìn)行稀疏表示的理想工具。因此,曲波變換目前已經(jīng)被廣泛應(yīng)用于多次波衰減[15]、地震數(shù)據(jù)去噪[16-17]、面波壓制[18]和圖像處理[19]等領(lǐng)域。

傳統(tǒng)的曲波變換閾值法利用有效信號和隨機噪聲在曲波系數(shù)上的分布差異,選擇合適的閾值函數(shù)對曲波系數(shù)進(jìn)行閾值處理,通過曲波反變換得到去噪后的結(jié)果。這種方法在去除噪聲的同時使信號邊緣附近產(chǎn)生了不光滑畸變現(xiàn)象,且對有效信號局部造成一定損傷。本文基于壓縮感知理論,在全變差(TV)準(zhǔn)則約束下,將傳統(tǒng)去噪問題轉(zhuǎn)化成求解最優(yōu)化問題,通過最優(yōu)解重構(gòu)原始信號,從而實現(xiàn)對地震資料的去噪處理。該方法能夠有效壓制同相軸邊緣附近產(chǎn)生的不光滑畸變,較好地衰減地震數(shù)據(jù)中的隨機噪聲,最大限度地保留有效信號,從而為地震數(shù)據(jù)的高保真處理奠定基礎(chǔ)。

1 基本原理

1.1 壓縮感知理論

信號的壓縮重建問題可以描述為由一組不完整的數(shù)據(jù)通過算子作用恢復(fù)出完整的數(shù)據(jù)，即：

(1)

其中,y∈RM為測得的不完整數(shù)據(jù),f∈RN為期望恢復(fù)的數(shù)據(jù),ψ∈RM×N(N>M)為測量矩陣,n∈RM為噪聲。

由于N>M,即方程的個數(shù)少于未知數(shù)的個數(shù),因此,這是一個欠定問題,也就是說方程(1)的解不唯一。壓縮感知理論為求取f提供了可能,但要求滿足3個條件：①f是稀疏的,或者在某個變換域內(nèi)是稀疏的;②采用隨機采樣,使假頻轉(zhuǎn)換成容易壓制的噪聲而將其濾除;③使用有效的求解策略尋找最優(yōu)解。

如果f在某個變換域中的系數(shù)x是稀疏的,即x=Φf(Φ為稀疏變換基矩陣),則(1)式可以寫成如下形式：

(2)

在壓縮感知理論框架下,信號的恢復(fù)問題可以表示為基于L0范數(shù)最小化的稀疏優(yōu)化問題，即：

(3)

(3)式本質(zhì)上是一個NP難問題[20],Donoho[21-22]采用最小L1范數(shù)代替最小L0范數(shù)的方法來求解,即：

(4)

信號f可以通過(5)式來重建：

(5)

在信號重建過程中,首先要保證測量矩陣ψ和稀疏變換基矩陣Φ非相干,其次,測量矩陣ψ要滿足有限等距特性[21]。因此,文中稀疏變換基矩陣選用曲波變換基矩陣,測量矩陣選用高斯隨機矩陣[22]。

1.2 稀疏變換

稀疏變換是壓縮感知的重要組成部分,變換域中的數(shù)據(jù)越稀疏,重構(gòu)效果就越好。本文使用對于地震數(shù)據(jù)的表示具有獨特優(yōu)勢的曲波變換作為稀疏變換。曲波變換正變換可以簡單地表示成信號f與曲波φj,l,k的內(nèi)積[19]：

(6)

由(6)式可知,給定一個曲波函數(shù),經(jīng)過伸縮、平移和旋轉(zhuǎn),可以生成平方可積函數(shù)空間的緊標(biāo)架,這就意味著它具有如下重構(gòu)公式：

(7)

1.3 基于壓縮感知和TV準(zhǔn)則約束的地震數(shù)據(jù)去噪原理

令含有噪聲的地震數(shù)據(jù)表示為：

(8)

式中：d為含噪數(shù)據(jù);s為有效信號;n為隨機噪聲。

首先通過一個M×N維(M?N)的測量矩陣ψ作用于含噪數(shù)據(jù)d,使其從N維(高維)信號投影到M維(低維)信號y:

(9)

然后對含噪數(shù)據(jù)d進(jìn)行稀疏變換：

(10)

將(9)式和(10)式整理得到：

(11)

對于含噪數(shù)據(jù)d而言,有效信號s在稀疏變換域內(nèi)是稀疏的,但是噪聲n是不稀疏的,根據(jù)壓縮感知理論,通過估計有效信號的稀疏表示來重建原始數(shù)據(jù),可以達(dá)到去噪的目的。上述過程等價于求解如下最優(yōu)化約束問題：

(12)

(13)

式中：n1和n2為x的維數(shù)。

這樣,根據(jù)拉格朗日法則,公式(12)可轉(zhuǎn)化成如下無約束的正則化問題[23]：

(14)

通過基于壓縮感知的TV重構(gòu)TVAL3算法[23]求解方程(14),地震數(shù)據(jù)去噪后的結(jié)果可以表示為：

(15)

2 數(shù)值模型試算

2.1 合成記錄A

首先建立如圖1a所示的合成地震記錄A,該記錄共50道,60個采樣點,采樣率為1ms,加入隨機噪聲之后的地震記錄如圖1b所示。圖2a為傳統(tǒng)的曲波變換閾值法去噪結(jié)果,圖2b為本文基于壓縮感知和TV準(zhǔn)則約束的曲波變換(以下簡稱為TV準(zhǔn)則約束的曲波變換)去噪結(jié)果。對比圖2a和圖2b可以發(fā)現(xiàn),兩種方法都能夠?qū)Φ卣鹩涗浿械碾S機噪聲進(jìn)行有效壓制,但是,傳統(tǒng)的曲波變換閾值法在去除噪聲的同時,使同相軸邊緣附近產(chǎn)生了不光滑畸變現(xiàn)象(橢圓區(qū)域所示),嚴(yán)重影響了去噪效果;而基于TV準(zhǔn)則約束的曲波變換能夠有效地壓制同相軸邊緣附近的不光滑干擾,去噪效果更好。此外,傳統(tǒng)的曲波變換閾值法使非有效信號區(qū)域變得模糊,而TV準(zhǔn)則約束的曲波變換在這些區(qū)域較好地保持了原始數(shù)據(jù)的特征。

圖1 合成記錄A(a)及加噪后記錄A(b)

圖2 傳統(tǒng)的曲波變換閾值法去噪(a)和TV準(zhǔn)則稀疏約束的曲波變換去噪(b)效果對比

2.2 合成記錄B

為了說明TV準(zhǔn)則約束的曲波變換具有更好的去噪效果和保幅能力,建立如圖3a所示的合成地震記錄B,該記錄共50道,512個采樣點,采樣率為2ms,圖3b為該記錄加入隨機噪聲后的結(jié)果。圖3c和圖3d分別為F-X反褶積和TV準(zhǔn)則約束的曲波變換去噪結(jié)果。對比圖3c和圖3d可見,前者雖然可以壓制絕大多數(shù)的隨機噪聲,但是去噪后的結(jié)果信噪比降低,去噪效果不夠理想;而TV準(zhǔn)則約束的曲波變換去噪結(jié)果不僅有效壓制了隨機噪聲,而且信噪比較高,去噪效果更好。圖3e 和圖3f為兩種方法去噪前、后的殘差剖面,可見F-X反褶積的殘差剖面上有效同相軸隱約可見,尤其是彎曲同相軸和相交同相軸(圖3e橢圓區(qū)所示),說明F-X反褶積對有效信號的損傷較大;相比而言,TV準(zhǔn)則約束的曲波變換去除的噪聲里幾乎沒有有效同相軸的痕跡,說明該方法在去噪過程中對有效信號的損傷相對較小。

為了更好地評價F-X反褶積與TV準(zhǔn)則約束的曲波變換去噪效果,定義如下參數(shù)：

1) 信噪比

(16)

2) 峰值信噪比

(17)

(16)和(17)式中：s0為無噪聲數(shù)據(jù);s為含噪聲數(shù)據(jù);M為數(shù)據(jù)矩陣的行數(shù);N為數(shù)據(jù)矩陣的列數(shù)。

3) 均方差

(18)

4) 平均絕對誤差

(19)

(18)和(19)式中：f(m,n)為去噪后數(shù)據(jù);f0(m,n)為無噪聲數(shù)據(jù);M為數(shù)據(jù)矩陣的行數(shù);N為數(shù)據(jù)矩陣的列數(shù)。

5) 去噪能力

(20)

其中,RSN為去噪后數(shù)據(jù)的信噪比，RSN0為原含噪數(shù)據(jù)的信噪比，EN0為原含噪數(shù)據(jù)中噪聲的能量，EN為去噪后數(shù)據(jù)中噪聲的能量。

表1是合成記錄B分別使用F-X反褶積和TV準(zhǔn)則約束的曲波變換去噪處理后計算得到的信噪比RSN,峰值信噪比RPSN,均方差SMS,平均絕對誤差SMA和去噪能力AD。由表1可知,TV準(zhǔn)則約束的曲波變換方法提高信噪比的能力更強,在信號保真方面也比F-X反褶積有優(yōu)勢。總體而言,TV準(zhǔn)則約束的曲波變換去噪能力(0.0661)明顯優(yōu)于F-X反褶積(0.0342),而且去噪過程中對有效信號的損傷也相對較小。

圖3 TV準(zhǔn)則約束的曲波變換與F-X反褶積去噪效果分析

表1 不同方法去噪結(jié)果評價參數(shù)對比

評價參數(shù)F?X反褶積TV準(zhǔn)則約束的曲波變換RSN/dB7．33478．9837RPSN/dB20．042521．8914SMS0．00440．0017SMA0．03850．0163AD0．03420．0661

3 實際資料處理效果分析

為了進(jìn)一步測試本文方法在實際地震資料處理中的應(yīng)用效果,選取某工區(qū)實際含噪聲地震資料進(jìn)行TV準(zhǔn)則約束的曲波變換,并與傳統(tǒng)的F-X反褶積進(jìn)行對比。圖4為研究區(qū)某實際單炮記錄,共150道,每道750個采樣點,采樣率為2ms。可以看出,隨機噪聲的存在使得該炮集數(shù)據(jù)的信噪比較低,同相軸連續(xù)性較差,尤其在800～1200ms,一些有效同相軸被完全淹沒在噪聲之中。圖5a 是基于TV準(zhǔn)則約束的曲波變換去噪結(jié)果,可見隨機噪聲得到了有效壓制,同相軸清晰、連續(xù)(圖5a中矩形框所示),位于800～1200ms的有效同相軸也得到恢復(fù),炮集數(shù)據(jù)整體信噪比得到明顯提高。圖5b為F-X反褶積去噪結(jié)果,可見大部分隨機噪聲也得到壓制,但是仔細(xì)對比圖5a與圖5b不難發(fā)現(xiàn)TV準(zhǔn)則約束的曲波變換去噪效果更好。圖6是TV準(zhǔn)則約束的曲波變換與F-X反褶積壓制的噪聲剖面,可見F-X反褶積壓制的噪聲剖面上彎曲同相軸和傾斜同相軸隱約可見,說明該方法對有效信號造成了一定程度損傷;而TV準(zhǔn)則約束的曲波變換壓制的噪聲剖面基本沒有有效同相軸的痕跡,說明該方法對有效信號的損傷較小。

圖4 某工區(qū)實際單炮記錄(含噪)

圖5 某工區(qū)含噪炮集TV準(zhǔn)則約束的曲波變換(a)與F-X反褶積(b)去噪結(jié)果對比

圖6 某工區(qū)含噪炮集TV準(zhǔn)則約束的曲波變換(a)與F-X反褶積(b)去除的噪聲對比

圖7為研究區(qū)去噪前疊加剖面,共1000道,時間為1000～1325ms,采樣率為2ms。圖8a為TV準(zhǔn)則約束的曲波變換去噪結(jié)果,可以看出1150～1250ms及1260～1310ms處(圖8a中矩形框所示)斷面變得更加清晰,一些淹沒在噪聲之中的同相軸也顯現(xiàn)出來,而且同相軸更加連續(xù)(圖8a 中橢圓框所示)。作為對比,應(yīng)用F-X反褶積對該數(shù)據(jù)進(jìn)行了處理(圖8b),兩種方法去掉的噪聲分別如圖9a和圖9b所示。由圖7至圖9可見,與疊前炮集去噪類似,TV準(zhǔn)則約束的曲波變換和F-X反褶積都可以壓制隨機噪聲,但是前者的去噪效果更好,對有效同相軸的損傷更小,從而為地震資料后續(xù)處理奠定了良好的基礎(chǔ)。

圖7 研究區(qū)含噪地震剖面

圖8 研究區(qū)含噪地震剖面TV準(zhǔn)則約束的曲波變換(a)與F-X反褶積(b)去噪結(jié)果對比

圖9 研究區(qū)含噪地震剖面TV準(zhǔn)則約束的曲波變換(a)與F-X反褶積(b)去除的噪聲對比

4 結(jié)論與認(rèn)識

本文將壓縮感知和TV準(zhǔn)則與曲波變換結(jié)合起來,提出了一種基于壓縮感知和TV準(zhǔn)則約束的曲波變換地震資料去噪方法,通過理論模型測試和實際資料處理,可以得到以下結(jié)論與認(rèn)識。

1) 與傳統(tǒng)的曲波變換閾值去噪方法相比,基于壓縮感知和TV準(zhǔn)則約束的曲波變換去噪方法在去除噪聲的同時可以很好地壓制同相軸邊緣附近產(chǎn)生的不光滑畸變,有效改善去噪效果。

2) 與F-X反褶積相比,基于壓縮感知和TV準(zhǔn)則約束的曲波變換去噪方法不僅能夠很好地壓制地震資料中的隨機噪聲,而且能夠減少對有效信號的損傷。

3) 無論對于炮集數(shù)據(jù)還是對于疊加剖面,基于壓縮感知和TV準(zhǔn)則約束的曲波變換都是一種有效的去噪方法。

基于曲波變換求解最優(yōu)化解是本文方法的核心,對計算機硬件有較高的要求,如何實現(xiàn)并行算法以提高計算效率,是我們需要進(jìn)一步研究的內(nèi)容。

[1] 康治,于承業(yè),賈臥,等.f-x域去噪方法研究[J].石油地球物理勘探,2003,38(2)：136-138 Kang Z,Yu C Y,Jia W,et al.A study on noise suppression method in f-x domain [J].Oil Geophysical Prospecting,2003,38(2):136-138

[2] Yu S,Cai X,Su Y.Seismic signal enhancement by polynomial fitting [J].Applied Geophysics,1989,1(1):57-65

[3] 張軍華,呂寧,田連玉,等.地震資料去噪方法綜合評述[J].石油地球物理勘探,2005,40(增刊)：121-127 Zhang J H,Lv N,Tian L Y,et al.Comprehensive review of seismic data denoising[J].Oil Geophysical Prospecting,2005,40(S1):121-127

[4] 劉偉,楊凱,王征,等.基于SVD的正交多項式變換及其在地震資料處理中的應(yīng)用[J].工程地球物理學(xué)報,2010,7(4):433-437 Liu W,Yang K,Wang Z,et al.Orthogonal polynomials transform based on SVD and its application in seismic data processing[J].Chinese Journal of Engineering Geophysics,2010,7(4):433-437

[5] 王勇,酈軍,張奎鳳.基于小波變換的地震信號降噪處理[J].石油物探,1998,37(3)：72-76 Wang Y,Li J,Zhang K F.The noise elimination processing of seismic signal based on wavelet transform [J].Geophysical Prospecting for Petroleum,1998,37(3):72-76

[6] Canales E J,Wakin M B.An introduction to compressive sampling [J].IEEE Signal Processing Magazine,2008,25(2):21-30

[7] Donoho D L.Compressed sensing [J].IEEE Transactions on Information Theory,2006,52(4):1289-1306

[8] 尹宏鵬,劉兆棟,柴毅,等.壓縮感知綜述[J].控制與決策,2013,28(10)：1441-1445 Yin H P,Liu Z D,Chai Y,et al.Survey of compressed sensing [J].Control and Decision,2013,28(10):1441-1445

[9] Shan H,Ma J W,Yang H Z.Comparisons of wavelets,contourlets and curvelets in seismic denoising[J].Journal of Applied Geophysics,2009,69(2):103-115

[10] Bochner S,Chandrasekharan K.Fourier Transforms [M].New Jersey:Princeton University Press,1949:45-52

[11] Feig E,Winograd S.Fast algorithms for the discrete cosine transform[J].IEEE Transactions on Signal Processing,1992,40 (9):2174-2193

[12] Chikkerur S,Cartwright A N,Govindaraju V.Fingerprint enhancement using STFT analysis[J].Pattern Recognition,2007,36(2):303-312

[13] Rioul O,Vetterli1 M. Wavelet and signal processing [J].IEEE Signal Processing,1991,8 (4):14-37

[14] Candès E J,Demannet S L,Donoho D L,et al.Fast discrete curvelet transforms[J].SIAM Multiscale Modeling and Simulation,2006,5(3):861-899

[15] Herrmann F J,Verschuur E.Curvelet-domain multiple elimination with sparseness constraints [J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1333-1336

[16] 張恒磊,張云翠,宋雙,等.基于Curvelet 域的疊前地震資料去噪方法[J].石油地球物理勘探,2008,43(5)：508-513 Zhang H L,Zhang Y C,Song S,et al.Curvelet domain-based prestack seismic data denoise method[J].Oil Geophysical Prospecting,2008,43(5):508-513

[17] 仝中飛,王德利,劉冰.基于Curvelet 變換閾值法的地震數(shù)據(jù)去噪方法[J].吉林大學(xué)學(xué)報(地球科學(xué)版),2008,38(增刊)：48-52 Tong Z F,Wang D L,Liu B.Seismic data denoise based on curvelet transform with the threshold method [J].Journal of Jilin University (Earth Science Edition),2008,38(S1):48-52

[18] 董烈乾,李振春,王德營,等.第二代Curvelet變換壓制面波方法[J].石油地球物理勘探,2012,46(6)：897-904 Dong L Q,Li Z C,Wang D Y,et al.Ground-roll suppression based on the second generation Curvelet transform [J].Oil Geophysical Prospecting,2012,46(6):897-904

[19] Ying L,Demanet L,Candes E.3D discrete curvelet transform [J].International Society for Optics and Photonics,2005,SPIE591413

[20] Amaldi E,Kann V.On the approximability of minimizing nonzero variables or unsatisfied relations in linear systems [J].Theoretical Computer Science,1998,209(1):237-260

[21] Donoho D.For most large underdetermined systems of linear equations,the minimal l1-norm near-solution approximates the sparsest near-solution [J].Communications on Pure and Applied Mathematics,2006,59(6):797-829

[22] Donoho D.Compressed sensing [J].IEEE Transactions on Information Theory,2006,51(4):1289-1306

[23] Li C B.An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing [D].USA：Rice University,2009

(編輯：戴春秋)

Random noise attenuation based on compressive sensing and TV rule

Liu Wei1,2,Cao Siyuan1,2,Cui Zhen1,2

(1.StateKeyLaboratoryofPetroleumResourcesandProspecting,ChinaUniversityofPetroleum,Beijing102249,China; 2.CNPCKeyLabofGeophysicalProspecting,ChinaUniversityofPetroleum,Beijing102249,China)

Owing to the undesired noise attenuation effect of conventional threshold method based on Curvelet transform,a seismic data denoising approach based on compressive sensing and Total Variation (TV ) rule is proposed in this paper.It changes noise attenuation into the optimization problem under the constraint of TV rule,the original signal is reconstructed by optimum solution and thereby the denosing results of seismic data are obtained.Theoretical numerical model test shows the proposed method can not only effectively suppress the distortion near the edges of seismic events but also achieve better results compared with conventional Curvelet transform threshold method.Moreover,compared with F-X deconvolution,this approach can attenuate the random noise effectively and preserve effective signal with an improvement of signal to noise ratio in seismic data.Finally,the denoising results of actual seismic data prove the effectiveness of the proposed method further.

compressive sensing,TV norm,Curvelet transform,random noise,noise attenuation

2014-09-28;改回日期：2014-12-18。

劉偉(1982—),男,博士在讀,主要從事地震數(shù)據(jù)處理新方法研究和儲層預(yù)測工作。

國家科技重大專項(2011ZX05024-001-01)和國家自然科學(xué)基金項目(41140033)共同資助。

P631

A

1000-1441(2015)02-0180-08

10.3969/j.issn.1000-1441.2015.02.009