適用于聲波方程數(shù)值模擬的時(shí)間-空間域隱式有限差分算子優(yōu)化方法

陳東， 梁文全， 辛維*， 楊長(zhǎng)春

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 勝利油田分公司東勝精攻石油開(kāi)發(fā)集團(tuán)股份有限公司， 山東東營(yíng)　257000

適用于聲波方程數(shù)值模擬的時(shí)間-空間域隱式有限差分算子優(yōu)化方法

陳東1,2， 梁文全1， 辛維1*， 楊長(zhǎng)春1

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所， 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 2 勝利油田分公司東勝精攻石油開(kāi)發(fā)集團(tuán)股份有限公司， 山東東營(yíng)257000

摘要聲波方程數(shù)值模擬已廣泛應(yīng)用于理論地震計(jì)算，同時(shí)構(gòu)成了地震逆時(shí)偏移成像技術(shù)的基礎(chǔ).對(duì)于有限差分法而言，在滿足一定的穩(wěn)定性條件時(shí)，普遍存在著因網(wǎng)格化而形成的數(shù)值頻散效應(yīng).如何有效地緩解或壓制數(shù)值頻散是有限差分方法研究的關(guān)鍵所在.有限差分格式分為顯式有限差分和隱式有限差分.隱式有限差分能夠進(jìn)一步壓制數(shù)值頻散效應(yīng).因此本文提出了給定頻率范圍滿足時(shí)間-空間域隱式有限差分頻散關(guān)系的方法，并根據(jù)震源頻率、波速和網(wǎng)格間距確定波數(shù)范圍，在此基礎(chǔ)上建立方程確定了相應(yīng)的隱式有限差分系數(shù)，使得差分系數(shù)能在更大頻率范圍符合波場(chǎng)傳播規(guī)律.通過(guò)頻散分析和正演模擬，驗(yàn)證了本文方法的有效性.

關(guān)鍵詞聲波正演； 時(shí)間-空間域； 隱格式有限差分； 頻散關(guān)系

1引言

因?yàn)槁暡ǚ匠逃邢薏罘址椒ㄓ?jì)算效率高、所需內(nèi)存相對(duì)較小、實(shí)現(xiàn)簡(jiǎn)單，而廣泛應(yīng)用于地震正演研究(Alford, et al.，1974；Kelly, et al.，1976；Basabe and Sen，2007；馮英杰等，2007；王珺等，2007；韓令賀等，2011；蘭海強(qiáng)等，2011；Chu and Stoffa, 2012a, 2012b; Wang, 2014)，同時(shí)也是地震逆時(shí)偏移成像技術(shù)得以快速發(fā)展的基礎(chǔ)(李博等，2012;劉紅偉等，2011,2012).

網(wǎng)格頻散是有限差分中至關(guān)重要的問(wèn)題，直接影響著有限差分法在波動(dòng)方程數(shù)值解過(guò)程中的應(yīng)用效果.網(wǎng)格頻散是由對(duì)時(shí)間和空間偏導(dǎo)數(shù)網(wǎng)格化造成的，相速度變成了網(wǎng)格間距的函數(shù)(Dablain 1986)，這會(huì)導(dǎo)致地震波的數(shù)值相速度不等于地球介質(zhì)的真實(shí)相速度，使得波場(chǎng)模擬的精度降低.一般來(lái)說(shuō)，如果存在時(shí)間頻散，則高頻的相速度增大；如果存在空間頻散，則高頻的相速度減小(Dablain 1986).

為了緩解或者壓制網(wǎng)格數(shù)值頻散效應(yīng)，一般采用兩種措施：其一是采用低階差分格式，要求時(shí)間步長(zhǎng)和空間步長(zhǎng)非常小，這會(huì)造成所需內(nèi)存和計(jì)算量過(guò)大而無(wú)法應(yīng)用于三維的情況；其二是采用高階差分格式，一般情況下，時(shí)間方向采用二階差分格式，空間方向采用高階差分格式或偽譜法.高階差分格式主要是利用等波紋準(zhǔn)則或最小誤差準(zhǔn)則優(yōu)化設(shè)計(jì)差分系數(shù)實(shí)現(xiàn)對(duì)空間偏導(dǎo)數(shù)的近似.偽譜法則是通過(guò)正、反傅里葉變換來(lái)實(shí)現(xiàn)空間偏導(dǎo)數(shù)的精確求解.兩種方法相比，高階空間差分格式具有精度適中、計(jì)算效率高的特征而得到廣泛應(yīng)用，特別是用于地震逆時(shí)偏移成像.偽譜法計(jì)算精度高，但因計(jì)算量過(guò)大，一般常用于地震正演模擬.此外，Dablain(Dablain 1986)和Chen(Chen 2007, 2011)提出了時(shí)間方向?qū)?shù)的四階有限差分格式，提高了波場(chǎng)模擬的精度.

在空間方向偏導(dǎo)數(shù)差分格式近似計(jì)算方面已開(kāi)展了大量研究工作，積累了大量的研究成果.考慮到在地震波場(chǎng)數(shù)值模擬中，需在時(shí)間方向和空間方向同時(shí)網(wǎng)格化，兩者共同決定著數(shù)值頻散特征.因此僅依靠空間方向偏導(dǎo)數(shù)計(jì)算精度的提高，并不能壓制時(shí)間頻散，甚至加強(qiáng)了時(shí)間頻散.鑒于上述原因，F(xiàn)inkelstein提出了在時(shí)間-空間域內(nèi)設(shè)計(jì)與確定有限差分格式的新方法，以期提高正演模擬計(jì)算精度.其方法原理是空間域?qū)Σ〝?shù)進(jìn)行泰勒展開(kāi)、特定頻率點(diǎn)滿足相速度關(guān)系、群速度關(guān)系的組合，用以控制精度的階數(shù)和頻散(Finkelstein and Kastner, 2006; Finkelstein and Kastner, 2008).

偏微分方程的有限差分格式分為顯式的和隱式的.顯式有限差分格式中，某點(diǎn)導(dǎo)數(shù)的計(jì)算需要該點(diǎn)和它周?chē)c(diǎn)的值.隱式有限差分格式中，某點(diǎn)導(dǎo)數(shù)的計(jì)算除了需要該點(diǎn)和它周?chē)c(diǎn)的值之外，還需要它周?chē)c(diǎn)的導(dǎo)數(shù)值(Liu and Sen，2009).使用隱式有限差分格式需要求解多對(duì)角陣，因而計(jì)算量較大.但是為了達(dá)到相同的精度，使用隱式有限差分算子比使用顯式有限差分算子一般需要更少的網(wǎng)格點(diǎn).Liu和 Sen(2009)的研究指出，對(duì)于二階導(dǎo)數(shù)來(lái)說(shuō)，隱式差分格式(2N+2)階有限差分算子的精度相當(dāng)于(4N+2)階顯式有限差分算子的精度.但是傳統(tǒng)的隱式有限差分有兩個(gè)缺點(diǎn)：第一是在空間域進(jìn)行泰勒展開(kāi)求取有限差分系數(shù)，僅能在低頻或者低波數(shù)較好的保持頻散關(guān)系；第二是在空間域求取隱式有限差分系數(shù)，不能減小時(shí)間頻散.因此，本文提出了使用線性化的方法在時(shí)間-空間域確定隱式有限差分算子，使得算子在更大的波數(shù)范圍保持時(shí)間-空間域頻散關(guān)系，從而提高正演模擬的精度和效率.

2方法

本文以一維聲波方程為例分析差分格式的優(yōu)化問(wèn)題：

(1)

其中，p(x,t)是波場(chǎng)，v是縱波速度.

二階導(dǎo)數(shù)的顯式有限差分算子公式：

(2)

函數(shù)p(x,t)二階中心差分：

(3)

二階導(dǎo)數(shù)的隱式有限差分算子公式為(Liu and Sen, 2009)

(4)

其中，cm(m=0,1,2,…,M)和b是隱格式有限差分系數(shù).當(dāng)b=0時(shí)，得到的是顯式的二階導(dǎo)數(shù)有限差分算子.

(5)其中，r=vτ/h， α=kh/2， ω=kv.

=(1-4bsin2α)α2,

(6)

則時(shí)間-空間域隱式有限差分格式退化為空間域隱式有限差分格式.

根據(jù)震源最高頻率，波速和網(wǎng)格間距確定需要優(yōu)化的波數(shù)上限(梁文全等, 2013)：

(7)

其中，f是最高震源頻率.

在確定的波數(shù)范圍內(nèi)，選擇M+1個(gè)均勻分布的波數(shù)點(diǎn)滿足時(shí)間-空間域隱式差分格式有限差分算子頻散關(guān)系(公式(5))，則得到下面的線性方程組：

(8)

解線性方程組，可以求得時(shí)間-空間域隱式有限差分算子的系數(shù)，其中

(9)

由公式(8)可見(jiàn)，時(shí)間-空間域隱式有限差分算子的系數(shù)與波速、所取波數(shù)范圍、時(shí)間步長(zhǎng)和空間步長(zhǎng)有關(guān).也就是說(shuō)，時(shí)間-空間域中，對(duì)于不同的波速，隱格式有限差分的系數(shù)b是不同的.但是隱格式有限差分中，空間偏導(dǎo)數(shù)的求取需要解一個(gè)三對(duì)角陣，這個(gè)三對(duì)角陣主要有b構(gòu)成，一般希望b值是固定的.同時(shí)我們知道，低波速需要的算子長(zhǎng)度最長(zhǎng).因此根據(jù)最低波速確定b值，然后把b作為已知，求取系數(shù)cm(m=0,1,2,…,M).空間域隱格式有限差分算子，空間域交錯(cuò)網(wǎng)格隱格式有限差分算子，時(shí)間-空間域隱格式交錯(cuò)網(wǎng)格有限差分算子的推導(dǎo)與上述相同，不再贅述.

3頻散分析

使用下面的公式衡量時(shí)間-空間域頻散(Liu and Sen, 2009)：

(10)

δ的值越接近1，數(shù)值頻散越小；否則數(shù)值頻散越大.

由圖1可以看到，隨著有限差分算子長(zhǎng)度的增加，傳統(tǒng)隱式有限差分算子和時(shí)間-空間域隱式有限差分算子的數(shù)值頻散都減小.時(shí)間-空間域新方法的數(shù)值頻散在整個(gè)kh范圍內(nèi)分布比較均勻，在更廣的kh范圍內(nèi)保持較小誤差，且可以根據(jù)震源頻率、空間步長(zhǎng)調(diào)整需要滿足頻散關(guān)系的波數(shù)上限使得kh比較小的情況下數(shù)值頻散也很小.

由圖2和3可以看出，隨速度和時(shí)間步長(zhǎng)的變化，傳統(tǒng)隱式有限差分算子的數(shù)值頻散變化比較大，而時(shí)間-空間域隱式有限差分算子頻散誤差變化比較小，因此新方法比傳統(tǒng)的方法穩(wěn)定.

4正演例子

4.1均勻介質(zhì)的二維聲波正演

首先使用均勻介質(zhì)模型，模型大小為220×200.對(duì)于所有的隱式有限差分算子M=8,v=2000 m·s-1,

圖1　兩種不同方法隨M增加的數(shù)值頻散圖(a) 頻散對(duì)比； (b) a圖的局部放大, 其中，v=2000 m·s-1時(shí)間步長(zhǎng)h=10 m,空間步長(zhǎng)τ=0.001 s.Fig.1　Plot of dispersion curves for different M (a) the comparison of FD dispersion; (b) zoom of (a)

圖2　兩種不同方法隨M增加的數(shù)值頻散圖(a) 頻散對(duì)比； (b) a圖的局部放大； 其中，v=2000 m·s-1時(shí)間步長(zhǎng)h=10 m,空間步長(zhǎng)τ=0.002 sFig.2　Plot of dispersion curves for different M when τ=0.002 s (a) the comparison of FD dispersion; (b) zoom of (a)

圖3　兩種不同方法隨M增加的數(shù)值頻散圖(a) 頻散對(duì)比； (b) a圖的局部放大；其中，v=4500 m·s-1時(shí)間步長(zhǎng)h=10 m,空間步長(zhǎng)τ=0.001 sFig.3　Plot of dispersion curves for different M when v=4500 m·s-1 (a) the comparison of FD dispersion; (b) zoom of (a)

h=10 m,時(shí)間步長(zhǎng) dt=1 ms.震源放在中心位置. 震源子波定義為

(11)

其中，f0=150 Hz，f0/π是震源主頻,t0=4/f0，c是常數(shù).

儒家思想之內(nèi)涵、處事原則及方式，無(wú)不體現(xiàn)著對(duì)經(jīng)與權(quán)的闡釋。權(quán)，從字面意思來(lái)理解是變化，似乎與原則相悖，但卻是對(duì)經(jīng)的完美補(bǔ)充，于是便有君子之中庸，有經(jīng)有權(quán)，小人反中庸。權(quán)，這一尺度生于人心權(quán)衡機(jī)能，生成固定之經(jīng)后，又可以反過(guò)來(lái)幫助人在具體情境中進(jìn)行尺度的判斷和運(yùn)用。

傳統(tǒng)隱式有限差分格式和時(shí)間-空間域隱式有限差分格式在350 ms時(shí)的波場(chǎng)快照如圖4所示.圖4的切片(z/dz=110)如圖5所示，我們以時(shí)間步長(zhǎng)0.5 ms的偽譜法為解析解.由圖4可見(jiàn)，偽譜法(τ=1 ms)的時(shí)間頻散嚴(yán)重，傳統(tǒng)的隱式差分格式時(shí)間頻散和空間頻散都嚴(yán)重，時(shí)間-空間域隱式有限差分格式同時(shí)有效壓制了時(shí)間頻散和空間頻散.4.2Marmousi模型

圖6顯示了Marmousi模型，其速度在1500 m·s-1和4700 m·s-1之間.震源子波如公式 (11)所示，其中f0=120 Hz.對(duì)于不同的有限差分算子h=10 m, dt=1 ms,M=8.

圖7a是空間域泰勒展開(kāi)法隱式有限差分格式的地震記錄，頻散明顯；圖7b是時(shí)間-空間域隱式有限差分算子的地震記錄, 時(shí)間和空間頻散都得到壓制.同時(shí)我們注意到，使用空間域隱式有限差分格式進(jìn)行正演模擬的時(shí)候，在1000 ms時(shí)已經(jīng)出現(xiàn)不穩(wěn)定現(xiàn)象，數(shù)值模擬崩潰.而使用時(shí)間-空間域隱式有限差分格式則可以保持穩(wěn)定.兩種方法在x/dx=250處的地震記錄如圖8所示，新方法明顯同時(shí)減小了時(shí)間頻散和空間頻散.

圖4　傳統(tǒng)隱式有限差分算子(a)和新時(shí)間-空間域有限差分算子(b)的波場(chǎng)快照Fig.4　Comparison of snapshots for the conventional implicit method and the new time-space domain implicit method

圖5　圖4的切片圖Fig.5　Slice of snapshots in Fig.4

圖6　Marmousi模型Fig.6　Marmousi velocity model

5結(jié)論

本文提出了時(shí)間-空間域確定聲波方程隱式有限差分算子的準(zhǔn)則，即給定波數(shù)范圍滿足頻散關(guān)系.方法本身無(wú)需在波數(shù)方向上進(jìn)行泰勒展開(kāi)，而是通過(guò)震源頻率、波速和網(wǎng)格間距確定波數(shù)上限，這使得新的差分格式在更大范圍符合波場(chǎng)傳播規(guī)律.此外，與傳統(tǒng)的隱式有限差分格式相比，新的隱式有限差分格式穩(wěn)定性更好.

圖7　Marmousi模型地震記錄(a) 傳統(tǒng)方法隱式有限差分算子; (b) 時(shí)間-空間域隱式有限差分算子.Fig.7　Seismic records from the Marmousi model(a) The conventional implicit method; (b) The new time-space domain implicit method.

圖8　Marmousi模型地震道對(duì)比紅線表示傳統(tǒng)隱式方法計(jì)算的地震道，綠線表示時(shí)空域隱式方法計(jì)算的地震道.Fig.8　Comparison of seismic seismograms from Marmousi model (the red line is the traditional implicit method, whereas the green line is the new time-space domain method).

通過(guò)頻散分析，以及二維的地震波數(shù)值模擬，驗(yàn)證了本文提出的新方法的有效性.因此，時(shí)間-空間域隱式有限差分算子方法可以代替空間域泰勒展開(kāi)法隱式有限差分算子用于地震的正演、逆時(shí)偏移中，以提高地震波數(shù)值模擬的精度和效率.

References

Alford R M, Kelly K R, Boore D M. 1974. Accuracy of finite-difference modeling of the acoustic wave equation.Geophysics, 39(6): 834-842.

Basabe J D D, Sen M K. 2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations.Geophysics, 72: T81-T95. Chen J B. 2007. High-order time discretizations in seismic modeling.Geophysics, 72(5): SM115-SM122.Chen J B. 2011. A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation.Geophysics, 76(2): T37-T42.

Chu C C, Stoffa P L. 2012a. Implicit finite-difference simulations of seismic wave propagation.Geophysics, 77(2): T57-T67.

Chu C C, Stoffa P L. 2012b. Determination of finite-difference weights using scaled binomial windows.Geophysics, 77(3): W17-W26.

Dablain M A 1986. The application of high-order differencing to the scalar wave equation.Geophysics, 51(1): 54-66.

Feng Y J, Yang C C, Wu P. 2007. The review of the finite-difference elastic wave motion modeling.ProgressinGeophysics(in Chinese), 22(2): 487-491.

Finkelstein B, Kastner R. 2008. A comprehensive new methodology for formulating FDTD schemes with controlled order of accuracy and dispersion.IEEETransactionsonAntennasandPropagation, 56(11): 3516-3525.Finkelstein B, Kastner R. 2006. Finite difference time domain dispersion reduction schemes.JournalofComputationalPhysics, 221(1): 422-438. Han L H, He B S, Zhang H X. 2011. Prestack reverse-time depth migration of quasi-P wave equations in VTI media.ActaSeismologicaSinica(in Chinese), 33(2): 209-218.Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms: A finite-difference approach.Geophysics, 41(1): 2-27.

Lan H Q, Liu J, Bai Z M. 2011. Wave-field simulation in VTI media with irregular free surface.ChineseJournalofGeophysics(in Chinese), 54(8): 2072-2084, doi: 10.3969/j.issn.0001-5733.2011.08.014.Li B, Li M, Liu H W, et al. 2012. Stability of reverse time migration in TTI media.ChineseJ.Geophys. (in Chinese), 55(4): 1366-1375, doi: 10.6038/j.issn.0001-5733.2012.04.032.

Liang W Q, Yang C C, Wang Y F, et al. 2013. Acoustic wave equation modeling with new time-space domain finite difference operators.ChineseJ.Geophys. (in Chinese), 56(10): 3497-3506, doi: 10.6038/cjg20131024.

Liu H W, Liu H, Zou Z, et al. 2010. The problems of denoise and storage in seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 53(9): 2171-2180, doi: 10.3969/j.issn.0001-5733.2010.09.017.

Liu H W, Liu H, Li B, et al. 2011. Pre-stack reverse time migration for rugged topography and GPU acceleration technology.ChineseJ.Geophys. (in Chinese), 54(7): 1883-1892, doi: 10.3969/j.issn.0001-5733.2011.07.022.Liu Y, Sen M K. 2009. A practical implicit finite-difference method: examples from seismic modelling.JournalofGeophysicsandEngineering, 6(3): 231-249.Wang J, Yang C C, Feng Y J. 2007. Utilizing optimal flux-corrected transport technology to suppress dispersion simulation.ProgressinExplorationGeophsics(in Chinese), 30(4): 252-256.Wang Y F, Liang W Q, Nashed Z, et al. 2014. Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method.Geophysics, 79(5): T277-T285.

附中文參考文獻(xiàn)

馮英杰, 楊長(zhǎng)春, 吳萍. 2007. 地震波有限差分模擬綜述. 地球物理學(xué)進(jìn)展, 22(2): 487-491.

韓令賀, 何兵壽, 張會(huì)星. 2011. VTI介質(zhì)中準(zhǔn)P波方程疊前逆時(shí)深度偏移. 地震學(xué)報(bào), 33(2): 209-218.

蘭海強(qiáng), 劉佳, 白志明. 2011. VTI介質(zhì)起伏地表地震波場(chǎng)模擬. 地球物理學(xué)報(bào), 54(8): 2072-2084, doi: 10.3969/j.issn.0001-5733.2011.08.014.

李博, 李敏, 劉紅偉等. 2012. TTI介質(zhì)有限差分逆時(shí)偏移的穩(wěn)定性探討. 地球物理學(xué)報(bào), 55(4): 1366-1375, doi: 10.6038/j.issn.0001-5733.2012.04.032.

梁文全, 楊長(zhǎng)春, 王彥飛等. 2013. 用于聲波方程數(shù)值模擬的時(shí)間-空間域有限差分系數(shù)確定新方法. 地球物理學(xué)報(bào), 56(10): 3497-3506, doi: 10.6038/cjg20131024.

劉紅偉, 劉洪, 鄒振等. 2010. 地震疊前逆時(shí)偏移中的去噪與存儲(chǔ). 地球物理學(xué)報(bào), 53(9): 2171-2180, doi: 10.3969/j.issn.0001-5733.2010.09.017.

劉紅偉, 劉洪, 李博等. 2011. 起伏地表疊前逆時(shí)偏移理論及GPU加速技術(shù). 地球物理學(xué)報(bào), 54(7): 1883-1892, doi: 10.3969/j.issn.0001-5733.2011.07.022.

王珺, 楊長(zhǎng)春, 馮英杰. 2007. 用優(yōu)化通量校正傳輸技術(shù)壓制數(shù)值模擬的頻散. 勘探地球物展, 30(4): 252-256.

(本文編輯汪海英)

Acoustic wave equation modeling based on implicit finite difference operators in the time-space domain

CHEN Dong1,2, LIANG Wen-Quan1, XIN Wei1*, YANG Chang-Chun1

1KeyLaboratoryofPetroleumResourcesResearch,InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China2DongshengJinggongPetroleumDevelopmentCo.Ltd.,ShengliOilfieldCompany,SINOPEC,ShandongDongying257000,China

AbstractNumerical simulation of acoustic wave equation is widely used to synthesize seismograms theoretically, and is also the basis of reverse time migration. With some stability conditions, dispersion often appears when time and the spatial derivatives in the wave equation are gridded. How to suppress such grid dispersion is therefore a key problem for finite difference approaches. Finite-difference formulations of partial differential equations are either explicit or implicit. Implicit formulations can suppress the grid dispersion further and the stability condition is better. Therefore, we propose to satisfy the time-space domain dispersion relationship in a proper range of frequencies. Dispersion analysis and seismic numerical simulation demonstrate the effectiveness of the proposed method.

KeywordsAcoustic wave equation modeling; Time-space domain; Implicit finite difference formulation; Dispersion relationship.

基金項(xiàng)目國(guó)家自然科學(xué)基金項(xiàng)目(41504142)資助. 國(guó)家高技術(shù)研究發(fā)展計(jì)劃(2011AA060301)，國(guó)家自然科學(xué)基金(41204086，41374122)項(xiàng)目聯(lián)合資助.

作者簡(jiǎn)介李振春，男，1963年生，教授，博士生導(dǎo)師，主要從事地震波正演及偏移成像方面的研究.E-mail：leonli@upc.edu.cn *通訊作者楊富森，男，碩士研究生.E-mail：sende12345@163.com 陳東，男，1966年生，高級(jí)工程師， 2001年于美國(guó)科羅拉多礦業(yè)學(xué)院獲石油工程碩士學(xué)位，目前在中石化勝利油田東勝公司從事開(kāi)發(fā)研究與管理工作.E-mail：chendong359.slyt@sinopec.com *通訊作者辛維，中國(guó)科學(xué)院地質(zhì)與地球物理研究所博士后.E-mail：xinwei@mail.iggcas.ac.cn

doi:10.6038/cjg20160429 中圖分類(lèi)號(hào)P631

收稿日期2015-08-27，2016-01-13收修定稿