網(wǎng)格剖分及邊界對(duì)MT Occam反演的影響研究

鐘藝文，劉　羽

(桂林理工大學(xué) a.地球科學(xué)學(xué)院，b.機(jī)械與控制工程學(xué)院，桂林　541004)

?

網(wǎng)格剖分及邊界對(duì)MT Occam反演的影響研究

鐘藝文a，劉羽b*

(桂林理工大學(xué) a.地球科學(xué)學(xué)院，b.機(jī)械與控制工程學(xué)院，桂林541004)

摘要：在大地電磁有限元模擬計(jì)算中網(wǎng)格剖分及邊界放置是否得當(dāng)，影響著最終的反演結(jié)果。基于Occam反演法從網(wǎng)格剖分及邊界方面，設(shè)計(jì)不同模型，討論不同網(wǎng)格對(duì)反演精度的影響。研究表明，對(duì)于二維介質(zhì)體，只要將網(wǎng)格邊界置于適當(dāng)倍數(shù)趨膚深度的距離處，便可忽略截?cái)噙吔绲挠绊懀瑹o(wú)需放置在無(wú)窮遠(yuǎn)。在同一頻段及不均勻剖分的前提下，粗細(xì)網(wǎng)格均可獲得良好的反演結(jié)果。細(xì)網(wǎng)格的反演精度高于粗網(wǎng)格，但未能很好地圈定異常體的邊界。粗網(wǎng)格剖分下，MT Occam法對(duì)低阻異常體反演效果較好，而對(duì)高阻異常體效果較差。

關(guān)鍵詞：有限單元法；網(wǎng)格剖分；邊界；MT Occam反演

0引言

有限單元法因理論完善、能模擬復(fù)雜結(jié)構(gòu)等優(yōu)點(diǎn)被廣泛應(yīng)用于地球物理場(chǎng)的數(shù)值模擬中。Coggon[1]最先將有限單元法引入大地電磁領(lǐng)域，經(jīng)Wannamaker[2]、Rodi[3]電磁場(chǎng)的有限元模擬的發(fā)展與完善，有限單元法進(jìn)入實(shí)用階段。在國(guó)內(nèi)陳樂(lè)壽[4]首先介紹了有限單元法在大地電磁正演計(jì)算中的應(yīng)用，此后徐世浙[5]等對(duì)限單元法進(jìn)行了深入研究和發(fā)展。如今有限單元法已成為二、三維大地電磁數(shù)值模擬的重要方法，但隨著計(jì)算機(jī)的高速發(fā)展，人們對(duì)反演的精度要求越來(lái)越高。因此許多學(xué)者[6-8]從算法精度方面(如插值函數(shù)、線(xiàn)性方程組的求解)進(jìn)行改進(jìn)研究，而在區(qū)域剖分方面的研究相對(duì)較少[10-13]。區(qū)域剖分是插值函數(shù)選取和解線(xiàn)性方程組的基礎(chǔ)，剖分是否合理所產(chǎn)生的誤差影響遠(yuǎn)遠(yuǎn)大于算法方面計(jì)算精度的高低，即使有再高精度的求解系統(tǒng)，網(wǎng)格剖分得不合理也未必能得到合理的結(jié)果。區(qū)域剖分單元通常有四邊形和三角形兩種形式，前者易于剖分具有通用性，后者可以模擬復(fù)雜地形及不規(guī)則地質(zhì)體。這里基于Occam反演法，區(qū)域剖分引用的是Wannamaker等[2]將矩形單元二次剖分成4個(gè)三角形網(wǎng)格的方案。對(duì)于網(wǎng)格剖分中邊界的選取傳統(tǒng)上都選擇放置于無(wú)窮遠(yuǎn)處，但實(shí)際模擬區(qū)域并未能取到無(wú)窮遠(yuǎn)，所以在網(wǎng)格的邊界處勢(shì)必存在截?cái)噙吔绲挠绊憽?duì)于網(wǎng)格邊界的選取，在追求獲取高精度結(jié)果的同時(shí)還得兼顧計(jì)算時(shí)間和內(nèi)存需求。作者在前人[10-14]研究基礎(chǔ)上，同樣用趨膚深度δ作為網(wǎng)格邊界選取的依據(jù)，依次改變網(wǎng)格的邊界，對(duì)模型的正反演計(jì)算進(jìn)行比較，總結(jié)網(wǎng)格邊界置于何處方可忽略截?cái)噙吔绲挠绊憽Ｕ菔欠囱莸幕A(chǔ)，作者先研究網(wǎng)格剖分對(duì)正演的影響；在網(wǎng)格剖分疏密對(duì)正演影響不大的基礎(chǔ)下，繼而研究總結(jié)網(wǎng)格剖分對(duì)反演的影響，為后續(xù)解釋工作提供依據(jù)。

1Occam反演法理論

地球物理中觀測(cè)數(shù)據(jù)是有限的、離散的，同時(shí)存在觀測(cè)誤差，因此造成地球物理反演具有不適定性。在眾多的反演方法中，Occam反演是能夠克服此不足的方法之一。Constable等[15-17]提出了Occam法反演理論，Occam反演法尋求的是最光滑反演模型，因此引入了粗糙度(R)項(xiàng)以壓制非數(shù)據(jù)產(chǎn)生的冗余構(gòu)造。Occam反演算法在計(jì)算反演模型m的過(guò)程中，不僅要求正演響應(yīng)數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)d的擬合差在設(shè)定精度內(nèi)，同時(shí)要求粗糙度盡可能小。因此，使用拉格朗日乘子μ來(lái)平衡模型光滑和數(shù)據(jù)擬合程度，最終所構(gòu)造的目標(biāo)泛函如式(1)所示。

(1)

(2)

(3)

其中：Nm是模型參數(shù)個(gè)數(shù)；Npt是約束項(xiàng)項(xiàng)數(shù)。

將式(1)中的問(wèn)題線(xiàn)性化，引用公式F[mk+△]=F[mk]+JK△,△=(mk+1-mk)表示的是模型修改量，mk是第k迭代的解，Jk為雅克比矩陣，即F[mk]相對(duì)于mk的偏導(dǎo)數(shù)。因此可以得到模型的迭代公式(式(4))。

mk+1(μ)=[μ(?T?+TT)+(WJk)TWJK]-1

(4)

2模型計(jì)算

2.1網(wǎng)格邊界的影響

在二維結(jié)構(gòu)中，在地面使用右旋制，令x軸為走向方向向北，y軸垂直x軸水平向東，z軸垂直向下。如圖1所示，假設(shè)u為所要求的電磁場(chǎng)，在不同模式下，令上邊界AB處的u都為常數(shù)“1”，值得一提的是TM模式下AB邊直接取在地面，而TE模式下，為了消除地下電磁場(chǎng)的變化對(duì)地面電磁場(chǎng)的影響，AB邊必須距離地面一定的距離。左右邊界AC、BD應(yīng)取在距橫向不均勻體足夠遠(yuǎn)的地方，使得在該處電磁場(chǎng)沿深度的分布可以被看成是與地下一維介質(zhì)或均勻半空間介質(zhì)時(shí)的情況相類(lèi)似。對(duì)于下邊界，要求異常場(chǎng)在CD邊上的場(chǎng)值為“0”，CD邊以下可看作是均勻半空間，所以CD邊滿(mǎn)足第三類(lèi)邊界條件。綜合以上所討論u滿(mǎn)足的邊界條件如式(5)所示。

(5)

圖1　研究區(qū)域及坐標(biāo)Fig.1　Study area and coordinates

在做邊界截取研究時(shí)，網(wǎng)格尺寸對(duì)計(jì)算結(jié)果有很大影響，所以在不改變網(wǎng)格尺寸的情況下，只對(duì)外邊界做調(diào)整。基于Wannamaker等[2]在PW2D的做法，TE模式下，將研究區(qū)域總橫向距離的一半作為AB邊的距離。MT Occam反演不依賴(lài)初始模型，通常取均勻半空間，電阻率取100 Ω·m。正演網(wǎng)格數(shù)為59(橫向)×17(縱向)，Occam要求反演網(wǎng)格是正演網(wǎng)格的子集，所以正反演網(wǎng)格必須是對(duì)齊的。反演網(wǎng)格設(shè)計(jì)：橫向按等間距1 000 m剖分并取兩個(gè)間距作為網(wǎng)格單元的寬，在此基礎(chǔ)上首尾兩網(wǎng)格單元的寬分別由±1 km、±3 km、±9 km、±27 km、±81 km、±243 km、±729 km七個(gè)間距組成；縱向網(wǎng)格由地面往下除前三個(gè)是按等間距剖分外，以下的都是按常數(shù)比例增長(zhǎng)剖分，縱向最后一層的厚度則是由多個(gè)縱向剖分間距組成。頻率分別為：0.937 5 Hz，0.625 Hz，0.468 8 Hz，0.312 5 Hz，0.234 4 Hz，0.156 3 Hz，0.117 2 Hz，0.078 13 Hz，0.058 59 Hz，0.039 06 Hz，0.029 3 Hz，0.019 53 Hz，0.014 65 Hz，0.009 77 Hz，0.007 33 Hz，0.004 88 Hz，0.003 66 Hz，0.002 44 Hz，0.001 83 Hz，0.001 22 Hz。

作者采用二維異常嵌入體模型做模擬計(jì)算，模型如圖 2所示。若嵌入體是低阻體，等于1 Ω·m。頻率取0.390 6 Hz，取背景電阻率為100 Ω·m，依據(jù)公式可知，趨膚深度大約是8 km。按上述的網(wǎng)格剖分左右邊界已置于足夠遠(yuǎn)。根據(jù)吳娟[10]、湯井田等[9]研究，只要下邊界置于1 δ距離以上，就不會(huì)影

響計(jì)算精度。因縱向剖分的層數(shù)足夠多，下邊界所在的地方也可以認(rèn)為是足夠遠(yuǎn)。二維介質(zhì)沒(méi)有解析解，所以認(rèn)為在此邊界條件下計(jì)算得的是相當(dāng)精確的數(shù)值解，所以把計(jì)算得到的這組數(shù)據(jù)作為參考解。

圖2　二維地質(zhì)體模型Fig.2　Two-dimensional geological model

首先在以上的網(wǎng)格剖分為基礎(chǔ)，將正演網(wǎng)格的左右邊界分別置于120 km、40 km、24 km處。在測(cè)點(diǎn)O處觀測(cè)得到的視電阻率如表1所示。

表1　兩側(cè)邊界變化下視電阻率觀測(cè)結(jié)果Tab.1　Apparent resistivity result when changing the both sides of mesh boundary

注：R120、R40、R24分別是左右邊界置于120 km、40 km、24 km處觀測(cè)的視電阻率；ε為誤差百分比。

從表1可知，TE模式下，當(dāng)網(wǎng)格邊界從1 000 km外不斷縮減到120 km處時(shí)，所觀測(cè)到的視電阻率沒(méi)有變化。當(dāng)縮減至40 km(5δ)、24 km時(shí)，相對(duì)高頻所觀測(cè)到的視電阻率的相對(duì)誤差在1 %左右。隨著頻率繼續(xù)減小，相對(duì)誤差增加到了5%左右。對(duì)TM模式作了相關(guān)研究，發(fā)現(xiàn)不管邊界取在120 km還是40 km或者24 km處，所觀測(cè)到的視電阻率相對(duì)誤差都維持在小于1 %的數(shù)量級(jí)。以上實(shí)驗(yàn)表明，只要將兩側(cè)邊界置于5 δ以上的距離處，便可以當(dāng)作不受截?cái)噙吔缬绊憽ｋm然在相對(duì)低頻處會(huì)存在一些小的偏差，但整體上可以反映淺部及深部的電性分布信息。

其次，在正演網(wǎng)格數(shù)59(橫向)×17(縱向)基礎(chǔ)上，反演網(wǎng)格數(shù)為24×9+13×3，兩側(cè)邊界分別置于120 km、40 km、24 km處，研究邊界的選取對(duì)反演結(jié)果的影響。

Occam反演所尋求的目標(biāo)擬合差是“1”，反演規(guī)定的最大迭代次數(shù)是20次。由表2的數(shù)據(jù)可知，反演的最終擬合差精度都相當(dāng)高，但是Occam反演結(jié)果尋求的是最光滑模型，所以還得考慮粗糙度這方面。隨著左右邊界的位置不斷拉近，反演模型的粗糙度會(huì)逐漸增加，反演不能正常收斂，以至于反演的迭代次數(shù)都達(dá)到限定次數(shù)時(shí)反演才結(jié)束，因此就會(huì)增加反演時(shí)間。綜合上述討論，截?cái)噙吔鐚?duì)正反演結(jié)果都有一定的影響，左右邊界都應(yīng)取在大于5δ距離處，但也無(wú)需置于無(wú)窮遠(yuǎn)。

表2　反演取不同網(wǎng)格邊界的反演結(jié)果Tab.2　Inversion results when inversion take different boundary of grid

2.2粗細(xì)網(wǎng)格對(duì)反演結(jié)果的影響

反演最終結(jié)果是每個(gè)單元上的電阻率值，求解方程組后將節(jié)點(diǎn)場(chǎng)值作垂向偏導(dǎo)，采用差分公式計(jì)算偏導(dǎo)數(shù)，由此可知網(wǎng)格剖分粗細(xì)與偏導(dǎo)數(shù)計(jì)算存在某種關(guān)系。以下將從網(wǎng)格的粗細(xì)做研究，依然采用如圖2所示模型，均勻半空間中有一個(gè)矩形異常體，可取1 Ω·m或者1 000 Ω·m。網(wǎng)格剖分成兩種形式：粗網(wǎng)格，正演網(wǎng)格數(shù)都是59(橫向)×35(縱向)，反演網(wǎng)格數(shù)為303，測(cè)點(diǎn)數(shù)為21，頻點(diǎn)數(shù)是20；細(xì)網(wǎng)格，正演網(wǎng)格數(shù)是103(橫向)×35(縱向)，反演網(wǎng)格數(shù)為889，測(cè)點(diǎn)數(shù)為21，頻點(diǎn)數(shù)是20。

反演是一個(gè)不斷正演的過(guò)程，為了量化影響反演的因素，首先討論了不同網(wǎng)格剖分對(duì)正演的影響。由圖3可知，低阻模型的低頻段粗細(xì)網(wǎng)格正演的數(shù)據(jù)有些許的偏差，整體上粗網(wǎng)格與細(xì)網(wǎng)格都獲得基本一致的正演的結(jié)果。由此說(shuō)明，只要網(wǎng)格橫縱向剖分的尺寸在某個(gè)范圍內(nèi)網(wǎng)格剖分的疏密對(duì)正演計(jì)算影響不大。

整個(gè)研究過(guò)程使用的是MT Occam 3.0版本程序，Occam 反演需要輸入四個(gè)文件，即DATA、MESH、MODEL、startup，DATA文件是經(jīng)正演而產(chǎn)生的，Occam反演的最終結(jié)果是反演網(wǎng)格數(shù)個(gè)電阻率，保存在迭代文件ITERXX中，這些電阻率是地下介質(zhì)的綜合反映，即某種平均下的電阻率，所以會(huì)有一定偏差。使用畫(huà)圖程序?qū)⒆罱K迭代文件中電阻率進(jìn)行繪圖便可得到圖4、圖5中的圖形。由圖4可知，不管是粗網(wǎng)格還是細(xì)網(wǎng)格在埋深5 km處都出現(xiàn)低阻異常，異常體中心電阻率值與實(shí)驗(yàn)值相差甚微，異常體中心位置都能與初始模型基本吻合，但細(xì)網(wǎng)格的模擬精度更優(yōu)于粗網(wǎng)格。而左右邊界方面，粗細(xì)網(wǎng)格異常體都未能很好的重現(xiàn)初始模型的邊界。細(xì)網(wǎng)格上邊界出現(xiàn)向上擴(kuò)撒現(xiàn)象，而粗網(wǎng)格沒(méi)有。粗細(xì)網(wǎng)格的下邊界都出現(xiàn)異常體向下擴(kuò)散延伸現(xiàn)象，細(xì)網(wǎng)格更是明顯。圖5是高阻模型下的反演結(jié)果，粗細(xì)網(wǎng)格反演的整個(gè)模擬區(qū)域都未出現(xiàn)明顯的高阻異常現(xiàn)象，很難判斷異常體的中心位置。出現(xiàn)此情況主要是由于在正演時(shí)產(chǎn)生的數(shù)據(jù)偏差較大加上繪圖的色彩較容易混淆。結(jié)合反演數(shù)據(jù)分析可知，在原設(shè)計(jì)異常體位置處反演得出的電阻率略微的高于周邊的電阻率，因此從數(shù)據(jù)方面是可以圈定異常體。

圖3　不同網(wǎng)格下測(cè)點(diǎn)O的視電阻率曲線(xiàn)對(duì)比圖Fig.3　Comparison of apparent resistivity curve at station O of different grid(a)低阻模型；(b)高阻模型

圖4　低阻模型下不同網(wǎng)格的MT Occam反演Fig.4　The MT Occam inversion of different grid in low resistance model(a)粗網(wǎng)格；(b)細(xì)網(wǎng)格

圖5　高阻模型下不同網(wǎng)格的MT Occam反演Fig.5　The MT Occam inversion of different grid in high resistance model(a)粗網(wǎng)格；(b)細(xì)網(wǎng)格

3結(jié)論

針對(duì)Occam框架下的MT反演問(wèn)題，通過(guò)將正反演網(wǎng)格的邊界置于不同位置，用不同反演網(wǎng)格作對(duì)比實(shí)驗(yàn)得出的結(jié)果可知：

1)用趨膚深度作為網(wǎng)格左右邊界劃分的度量，正反演網(wǎng)格的左右邊界置于5倍趨膚深度以上但無(wú)需置于無(wú)窮遠(yuǎn)，便可忽略截?cái)噙吔绲挠绊憽?/p>

2)在同一頻段內(nèi)，細(xì)網(wǎng)格的模擬精度略高于粗網(wǎng)格，而細(xì)網(wǎng)格在圈定異常體范圍方面稍有欠缺。粗網(wǎng)格剖分下Occam反演法對(duì)低阻異常相對(duì)靈敏，反演結(jié)果精度相當(dāng)高，而反演高阻異常體得的結(jié)果較差，與真實(shí)相差甚遠(yuǎn)。

3)基于淺部剖分都相對(duì)比較細(xì)，而深部則是相對(duì)較粗的不均勻剖分粗細(xì)網(wǎng)格兩種形式都可獲得良好的反演效果。

參考文獻(xiàn)：

[1]COGGON J H.Electromagnetic and electrical modeling by the finite element method[J].Geophysics,1971,36(1):132-155.

[2]WANNAMAKER P.E，STODT J.A.A stable finite element solution for two-dimensional magnetotelluric modelling[J].Geophysical Journal of the Royal Astronomical Society，1987：277-296.

[3]RODI W L.A technique for improving the accuracy of finite element solution for magnctotelluric data [J].Geophysics，1976,44(2)：483-506.

[4]陳樂(lè)壽.有限元法在大地電磁場(chǎng)正演計(jì)算中的應(yīng)用于改進(jìn)[J].石油物探，1981,20(3):84-103.

CHEN L S.Application and improvement of finite-element method in forward calculation of geo-electromagnetic field[J].geophysical prospecting for petrole,1981,20(3):84-103.(In Chinese)

[5]徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社，1994.

XU S Z.The finite element method in geophysics [M].Beijing:Science Press,1994.(In Chinese)

[6]陳小斌,胡文寶．有限元直接迭代算法在MT二維正演計(jì)算中的應(yīng)用[J]．石油地球物理勘探，2000,35(4):487-496.

CHEN X B,HU W B.Application of finite-element direct iteration algorithm to mt 2-d forward computation[J].oil geophysical prospecting，2000，35(4):487-496.(In Chinese)

[7]馬為，陳小斌.大地電磁測(cè)深二維正演中輔助場(chǎng)的新算法[J].地震地質(zhì)，2008,30(2):525-533.

MA W,CHEN X B.A new algorithm for the calculation of auxiliary field in mt 2-d forward modeling [J].Seismology and Geology,2008,30(2):525-533.(In Chinese)

[8]柳建新，郭榮文，童孝忠，等.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用[J].中南大學(xué)學(xué)報(bào)，2009,40(2):484-491.

LIU J X，GUO Z W,TONG X Z,et al.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling [J].Journal of central university，2009,40(2):484-491.(In Chinese)

[9]湯井田,薛帥.MT有限元模擬中截?cái)噙吔绲挠绊慬J].吉林大學(xué)學(xué)報(bào)，2013,43(1):267-274.

TANG J T,XUE S.Influence of Truncated Boundary in FEM Numerical Simulation of MT [J].Journal of Jilin university,2013,43(1):267-274.(In Chinese)

[10]吳娟，席振銖，王鶴.網(wǎng)格尺寸及邊界對(duì)大地電磁有限元正演精度的影響[J].物探化探計(jì)算技術(shù)，2012,34(1):27-32.

WU J,XI Z Z，WANG H.Effect of grid size and boundary on MT forward modeling using finite element method [J].Computing Technology of Geophysical and Geochemical Exploration,2012,34(1):27-32.(In Chinese)

[11]朱崇利，董云,王延平，等.網(wǎng)格剖分對(duì)大地電磁測(cè)深反演精度的影響[J].物探與化探，2014,38(4):737-741.

ZHU C L,DONG Y，WANG Y P，et al.The effect of grid subdivision on the accuracy of MT inversion [J].Geophysical and Geochemical Exploration,2014,38(4):737-741.(In Chinese)

[12]朱崇利.網(wǎng)格剖分對(duì)反演的影響[J].地球物理學(xué)進(jìn)展，2014,29(2):889-893.

ZHU C L.The influence of grid subdivision on inversion [J].Progress in Geophysics,2014,29(2):889-893.(In Chinese)

[13]歐東新.計(jì)算機(jī)精度和網(wǎng)格大小對(duì)大地電磁有限單元法正演的影響[J].桂林工學(xué)院學(xué)報(bào)，2007,27(3):329-332.

OU D X.Effect of computer precision and grid length on MT simulating using FEM [J].Journal of Guilin Institute of Technology，2007,27(3):329-332.(In Chinese)

[14]SHARMA P S,KAIKKONEN P.An automated finite element mesh generation and element coding in 2-D electromagnetic inversion[J].Geophysics，1998,34(3):93-114.

[15]CONSTABLE S C,PARKER R L,CONSTABLE C G.Occam’s inversion:a practical algorithm for generating smooth models from electromagetic sounding Data [J] Geophysics ,1987,52(3):289-300.

[16]DEGROOT-HEDLIN C,CONSTABLE S C.Occam’s inversion to generate smooth two dimensional models from magnetotelluric data[J].Geophysics ,1990,55(12):1613-1624.

[17]DEGROOT-HEDLIN C,STEVEN CONSTABLE .Occam’s inversion and the north American central plains electrical anomaly[J].Geomag.eoelectr,1993,45:985-99.

The study of effect of mesh generation and boundary on MT Occam inversion

ZHONG Yi-wena,LIU Yub*

(Guilin University of Technology a.College of Earth Sciences,b.College of Mechanical and Control Engineering,Guilin541004,China)

Abstract:The grid subdivision is reasonable or not and the boundary select appropriately of MT modeling using finite element method have influence on the inversion results.Based on Occam inversion method,using different models from grid subdivision and boundary settings,discussed their effect to the inversion results.The result show that mesh boundary are not necessary to place at infinity only in the appropriate multiple skin depth distance can ignore the influence of the truncation boundary on two-dimensional media.At the premise of same frequency band and non-uniform subdivision,good inversion results can be obtained from both fine mesh and coarse mesh.The inversion accuracy of fine grid better than that of the coarse grid,but not good at delineating the boundary of the abnormal body.The MT Occam results of coarse grid for the abnormal body of low resistance is better,however for abnormal body of high resistance,the results is worse.

Key words:finite element method;mesh generation;boundary;MT Occam inversion

收稿日期：2015-04-22改回日期：2015-06-09

基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目(41264005)

作者簡(jiǎn)介：鐘藝文(1990-)，男，碩士，研究方向?yàn)殡姶艛?shù)值模擬，E-mail：zhongyw321 @163.com。*通信作者：劉羽(1961-)，男，教授，博士，主要研究方向?yàn)閿?shù)據(jù)處理及并行計(jì)算，E-mail：lewis_5709 @163.com。

文章編號(hào)：1001-1749(2016)02-0145-06

中圖分類(lèi)號(hào)：P 631.3

文獻(xiàn)標(biāo)志碼：A

DOI：10.3969/j.issn.1001-1749.2016.02.01