郭楚楓,張世暉,劉天佑

(中國地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

有限單元法是地球物理數(shù)值正演模擬方法的一種，以變分問題為基礎(chǔ)，根據(jù)地球物理中的偏微分方程和邊界條件，用數(shù)值方法近似計算地球物理場值，適用于復(fù)雜物性分布和復(fù)雜邊界形狀的地球物理計算[1]。自1971年Coggon首次系統(tǒng)地闡述了有限單元法在電、電磁法中的應(yīng)用以來[2]，有限單元法已經(jīng)廣泛應(yīng)用于復(fù)雜形態(tài)、復(fù)雜地質(zhì)條件的地球物理正演中。

在電磁法領(lǐng)域，基于非結(jié)構(gòu)化網(wǎng)格的矢量有限單元法，具有很強(qiáng)的適應(yīng)性，是目前國內(nèi)外研究的熱點[3-5]。隨著研究的深入，各向同性假設(shè)模型已難以滿足解釋的需要，基于有限元的各向異性研究受到了廣泛關(guān)注[6-8]。在重力勘探中，有限單元法常用來進(jìn)行梯度張量及重力矢量的正演模擬計算[9-12]，Cai和Wang[11]實現(xiàn)了利用有限元法計算非均勻復(fù)雜形體重力異常的快速算法，朱自強(qiáng)等[13]則利用有限單元法進(jìn)行重力異常的地形校正。在磁場正演模擬中，有限單元法的優(yōu)勢在于不需要引入均勻磁化設(shè)定，能同時考慮磁各向異性、退磁效應(yīng)、剩磁三種因素的影響[14]，所以在充分考慮退磁作用、進(jìn)行高精度磁測資料反演解釋時，常用有限單元法進(jìn)行正演模擬[15-16]。

上述的有限元計算存在相同的問題，即采用傳統(tǒng)截斷邊界，在一個相對較大的區(qū)域內(nèi)，將無窮的邊界問題近似為有限區(qū)域進(jìn)行[17]。而有限單元法是一種區(qū)域型方法，必須在全區(qū)域進(jìn)行剖分，對于地球物理正演中的無界區(qū)域問題，在邊界處產(chǎn)生畸變會影響對異常的解釋。一般的解決方法是進(jìn)行擴(kuò)邊，即設(shè)定更大的求解區(qū)域以減小傳統(tǒng)截斷邊界的影響。如蔣甫玉等[9]進(jìn)行三度體重力矢量有限元正演計算時提出，設(shè)置求解域邊界長度應(yīng)大于場源體長度的7.5倍以達(dá)到理想效果；Ren和Kalscheuer等[5]進(jìn)行三維有限元大地電磁正演時，設(shè)置了邊長為20 km的求解域，以保證長度100 m的板狀體模型正演能達(dá)到模擬精度。擴(kuò)邊往往造成有限元計算區(qū)域大，節(jié)點數(shù)多，進(jìn)一步降低了有限元正演的計算效率[17]，較合理的策略是將計算區(qū)域限定在較小的目標(biāo)區(qū)域范圍內(nèi)，減少單元與節(jié)點數(shù)量，以節(jié)約計算資源，提高正演效率。利用無限單元取代傳統(tǒng)的邊界條件，通過坐標(biāo)映射將“無限單元”沿某一方向無限擴(kuò)展，最后結(jié)合有限單元法對內(nèi)部區(qū)域進(jìn)行正演模擬，提高正演效率。

無限元的概念最早由Ungless在1973年給出，并實現(xiàn)了衰減型無限單元[18]。無限元大體可以分為3類：Bettess映射無限元[19-20]、Astley波包映射無限元[21-22](mapped wave envelope elements)和Burnett多級展開無限元[23-24]，其中Astley波包映射無限元應(yīng)用最廣。無限元的基本思路是通過坐標(biāo)變換，采用和位場不同的插值函數(shù)，通過坐標(biāo)映射實現(xiàn)無限域積分，并在映射后的局部坐標(biāo)里建立位場插值形函數(shù)，使位場在有限元網(wǎng)格和無限元網(wǎng)格結(jié)合處保持一致性，而在無限遠(yuǎn)處衰減為0[25]。目前，無限元廣泛應(yīng)用于聲學(xué)、巖土力學(xué)等方面[26-27]，在地球物理學(xué)中則主要應(yīng)用于地震、測井和電法領(lǐng)域：Fu和Wu[28]應(yīng)用無限元處理了彈性波的吸收邊界條件；朱軍等[29]則利用無限元模擬測井中面臨的無窮邊界問題；湯井田等[30]應(yīng)用無限元進(jìn)行了三維直流電阻率的數(shù)值模擬；張林成等[17]應(yīng)用無限元進(jìn)行了基于二次場的可控源電磁法的三維數(shù)值模擬。目前，混合有限元—無限元在磁場正演模擬中應(yīng)用較少。為提高有限元三維磁場正演模擬的效率，本文基于COMSOL Multiphysics軟件，在求解域外部邊界上設(shè)置無限元，綜合考慮剩磁、退磁及地表起伏，通過一系列的數(shù)值模擬，與解析算法和傳統(tǒng)有限元方法進(jìn)行對比，一方面驗證了方法的有效性，另一方面證明本算法在保證精度的同時能夠大大提高正演效率。

1 方法理論

1.1 有限單元法原理

對于地球物理中磁場的正演問題，其邊界條件是依據(jù)麥克斯韋方程組給出的：

(1)

即在無電流源的磁場空間中磁場強(qiáng)度H和磁感應(yīng)強(qiáng)度B，應(yīng)該滿足：

H和B的關(guān)系為：

B=μ0(H+M)，

其中：μ0為真空中的磁導(dǎo)率，M為磁化強(qiáng)度，由感應(yīng)磁化強(qiáng)度Mi和剩余磁化強(qiáng)度Mr組成

M=Mi+Mr，

而感應(yīng)磁化強(qiáng)度Mi與磁場強(qiáng)度H成正比

Mi=χH，

其中：χ為磁化率，由以上關(guān)系式可以推出

(2)

圖1 非均勻介質(zhì)分布(改自徐世浙，1994)Fig.1 Distribution of inhomogeneous medium(modified from Xu Shizhe，1994)

帶入式(2)中，得

(3)

這是包括感磁和剩磁的磁標(biāo)位基本方程。

圖1為磁場正演中非均勻介質(zhì)分布的簡單模型，Γ0為兩介質(zhì)共有邊界，Γ1、Γ2為兩介質(zhì)的外部邊界，Γ∞為區(qū)域外部邊界，Vm1、Vm2、μr1、μr2、Mr1、Mr2分別代表界面兩側(cè)介質(zhì)中的磁位、相對磁導(dǎo)率和剩余磁化強(qiáng)度。徐世浙[1]已經(jīng)證明，磁位邊值問題中的內(nèi)部邊界條件屬于自然邊界條件，使用有限單元法解偏微分方程(3)時，無需討論磁位的內(nèi)部邊界條件。因此有限單元法具有解決非均勻介質(zhì)的磁場正演能力。

給定Vm的外邊界條件，將計算區(qū)域Ω取足夠大，使得區(qū)域的邊界Γ∞與磁性體足夠遠(yuǎn)，這時磁性體引起的磁位在邊界上近似為零。但在磁法勘探中存在著正常地磁場T，?！奚系拇盼皇怯蒚決定的，所以?！奚嫌械谝活愡吔鐥l件：

Vm=-T·r， ∈Γ∞；

(4)

式中：r是坐標(biāo)原點至邊界點的矢徑。轉(zhuǎn)為第二類邊界條件

(5)

針對基本方程(3)及邊界條件(5)，可以構(gòu)造泛函

(6)

由于內(nèi)部邊界條件為自然邊界條件，考慮?！薜倪吔鐥l件式(5)，以及自由空間中μr=1，Mr=0，式(6)的變分為

因此，無電流靜磁場的磁位的邊值問題與下列的變分問題相應(yīng)：

δF(Vm)=0。

(7)

1.2 有限元—無限元耦合算法

有限元—無限元耦合算法將整個求解區(qū)域分為有限元區(qū)域和無限元區(qū)域，用無限元區(qū)域代替?zhèn)鹘y(tǒng)的外邊界，在兩個區(qū)域分別采取有限元分析和無限元分析，通過總體剛度矩陣組裝將兩者結(jié)合起來從而進(jìn)行數(shù)值求解。

圖2a中劃分了有限元和無限元計算區(qū)域，圖中有限元區(qū)域為目標(biāo)區(qū)域，包含場源、目標(biāo)體及測點等；無限元區(qū)域為有限元區(qū)域邊界至無窮遠(yuǎn)，為邊界計算區(qū)域。無限元分析就是在該區(qū)域的某一方向上，采用無限元映射和形函數(shù)，將整體坐標(biāo)映射到局部坐標(biāo)上來，其原理與有限元分析相同。

使用有限單元法計算時，采用四面體單元網(wǎng)格剖分，線性插值單元中的插值函數(shù)為：

(8)

式中：Ni為插值形函數(shù)，與四面體自然坐標(biāo)形同；Vmi為節(jié)點上的磁標(biāo)位值。

圖2 三維無限元映射(改自湯井田等，2010)Fig.2 3-D infinite element mapping(modified from Tang et al.，2010)

無限元單元分析時，采用六節(jié)點，圖2a給出了三維無限元映射示意，圖中節(jié)點1、2、3、4、5、6為無限元基本要素，P點為坐標(biāo)原點，節(jié)點1、2、3為有限元邊界上某單元的三節(jié)點，P點到節(jié)點4、5、6的距離分別為P點到節(jié)點1、2、3距離的兩倍。無限單元最外圍的3個節(jié)點為無窮遠(yuǎn)。無限元分析就是通過無限元映射將無窮坐標(biāo)映射到圖中的局部坐標(biāo)上。

圖2b中ζ方向表示無限元映射的方向，在ξ-η平面內(nèi)，無限單元與有限單元采用相同的映射形式和形函數(shù)。無限元的坐標(biāo)映射為

式中：Li為三角形面積坐標(biāo)。結(jié)合圖2有：

(9)

結(jié)合ζ方向的坐標(biāo)映射關(guān)系式，可以得到六節(jié)點無限單元的結(jié)點坐標(biāo)映射函數(shù)：

(10)

其中：

經(jīng)過上述映射之后，磁位借助無限單元延伸至無限遠(yuǎn)處，并且在無限遠(yuǎn)處衰減為零。對于無限元形函數(shù)，其具體的表達(dá)式如下：

(11)

該形函數(shù)是Astley映射無限元理論[21]中采用的形函數(shù)中的系數(shù)(1-ζ)/2與二階Lagrange插值多項式的乘積。

對有限元—無限元進(jìn)行單元分析后，對方程組進(jìn)行總體合成，再解相應(yīng)的大型線性方程組即可得到各個節(jié)點的磁位。得到各節(jié)點的磁位后，通過磁場強(qiáng)度與磁位的關(guān)系，求得磁場的各個分量：

(12)

2 模型試算

為說明有限元—無限元算法的優(yōu)勢，綜合考慮剩磁和退磁的影響，設(shè)計了兩種模型：球體模型和組合體模型。使用COMSOL Multiphysics軟件進(jìn)行建模，以下數(shù)值模擬平臺均為Intel(R) Xeon(R) CPU3.91G，64GB RAM，16CPUs。

2.1 球體模型

在磁場正演模擬中，忽略退磁和剩磁會對正演結(jié)果產(chǎn)生較大的影響，忽略剩磁的相對誤差與Q(科尼斯布格比)相關(guān)，而忽略退磁的相對誤差與磁化率相關(guān)，若Q=0.5，忽略剩磁的相對誤差可達(dá)到33%以上[31]。相比于一般的正演方法，有限單元法的優(yōu)勢在于適用于復(fù)雜條件、復(fù)雜形狀的地球物理計算。為突出有限單元法的優(yōu)勢，驗證基于有限元—無限元耦合的磁場正演算法的準(zhǔn)確性，設(shè)計了如圖3所示的模型。由于球體的退磁系數(shù)為常數(shù)，可以利用退磁改正公式計算磁異常，將傳統(tǒng)有限元方法，有限元—無限元耦合方法和解析解進(jìn)行對比。

圖3 球體模型及網(wǎng)格剖分示意(藍(lán)色為球體模型，紅色直線代表測線，黃色為觀測平面范圍)Fig.3 Diagram of the Sphere model and mesh generation(the blue region is the sphere, the red line is observation line, the yellow plane is the range of observation plane)

在相同的求解域設(shè)置和相同網(wǎng)格剖分條件下，分別使用有限元—無限元耦合算法和傳統(tǒng)有限元方法進(jìn)行數(shù)值模擬。圖4a為經(jīng)過退磁改正后的球體ΔT磁異常特征圖，圖4b為有限元—無限元耦合算法的數(shù)值模擬結(jié)果，圖4c為傳統(tǒng)有限單元法數(shù)值模擬結(jié)果。圖4d、圖4e為平面誤差分布。圖4b中基于有限元—無限元耦合算法的結(jié)果與圖4a中解析解的異常形態(tài)特征及幅值基本一致，圖4d中，有限元—無限元耦合算法在邊界處幾乎沒有產(chǎn)生誤差。傳統(tǒng)有限元數(shù)值模擬結(jié)果與解析解差別較大，圖4c和圖4e中，求解域邊長為200 m的傳統(tǒng)有限元數(shù)值模擬得到的ΔT磁異常幅值與解析解相比，最大偏差達(dá)到837 nT，而在邊界處異常形態(tài)差別較大。此例說明，傳統(tǒng)有限單元法進(jìn)行數(shù)值模擬時在邊界處會產(chǎn)生較大的誤差，嚴(yán)重影響了正演模擬的精度。而在相同的計算范圍內(nèi)，基于有限元—無限元耦合的算法的數(shù)值模擬結(jié)果非常優(yōu)秀，特別在邊界附近達(dá)到了很高的精度，明顯優(yōu)于傳統(tǒng)有限元算法的結(jié)果。

擴(kuò)大求解域的范圍，利用傳統(tǒng)有限單元法進(jìn)行數(shù)值模擬，再與有限元—無限元耦合算法進(jìn)行對比。求解域一的長寬高分別為300、300、300 m，求解域二為400、400、400 m。圖5為相同磁場參數(shù)的球體模型，在相同網(wǎng)格剖分條件下，使用不同方法計算得到的ΔT磁異常曲線。曲線1為解析解，曲線2為基于有限元—無限元的計算結(jié)果，曲線3為求解域邊界長度為200 m的傳統(tǒng)有限元計算結(jié)果，曲線4為求解域邊界長度為300 m的傳統(tǒng)有限元計算結(jié)果，曲線5為求解域邊界長度為400 m的傳統(tǒng)有限元計算結(jié)果。

圖4 不同方法正演的球體平面ΔT磁異常及其平面誤差分布Fig.4 The plane total-field anomaly of the sphere using different methods

圖5 不同方法正演的球體磁異常曲線Fig.5 The total-field anomaly curve of the sphere using different methods

在進(jìn)行三維有限元數(shù)值模擬時，由于受到6個方向的截斷邊界影響，使用傳統(tǒng)有限單元法進(jìn)行三維正演時，若要提高精度，需要在空間的6個方向?qū)δＰ瓦M(jìn)行延展。圖5和圖6中，對比解析解(曲線1)和傳統(tǒng)有限元(曲線3、曲線4和曲線5)的結(jié)果，隨著求解區(qū)域的逐漸增大，傳統(tǒng)有限單元法在觀測平面邊界處產(chǎn)生的誤差逐漸減小，正演精度也在逐漸提高。當(dāng)求解域邊界長度達(dá)到400 m時，即對應(yīng)曲線5的結(jié)果，傳統(tǒng)有限元的數(shù)值模擬結(jié)果和解析計算結(jié)果已經(jīng)十分接近，但是仍然存在約85 nT的平均絕對誤差，且隨著求解區(qū)域的擴(kuò)大，傳統(tǒng)有限單元法的精度增加有限，圖5中曲線5和曲線1也不重合。圖5中，解析解對應(yīng)的曲線1和有限元—無限元數(shù)值模擬對應(yīng)的曲線2幾乎完全重合，基于有限元—無限元算法所產(chǎn)生的絕對誤差在0左右，說明了基于有限元—無限元算法的準(zhǔn)確性。對比有限元—無限元(曲線2)和求解域邊長為200 m的傳統(tǒng)有限元(曲線3)的數(shù)值模擬結(jié)果，在相同的區(qū)域內(nèi)，基于有限元—無限元的算法(曲線2)，在邊界處所產(chǎn)生的絕對誤差基本為0，基本避免了截斷誤差，大大改善了邊界處的計算結(jié)果，再次說明基于有限元—無限元算法相對于傳統(tǒng)有限元算法，在邊界處能夠得到更高的精度。

圖6 絕對誤差分布曲線Fig.6 The absolute error distribution curve

不同有限元算法模型的計算效率和誤差對比如表1所示?；谟邢拊獰o限元的算法需要對外部設(shè)置的無限元區(qū)域進(jìn)行網(wǎng)格剖分，所以相對初始模型會生成更多的網(wǎng)格節(jié)點。由表1可知，當(dāng)求解域邊界長度設(shè)置為400 m時，傳統(tǒng)有限元算法的平均相對誤差為8.4%，均方根誤差為86.1 nT，最大絕對誤差為120.9 nT，相對于初始模型的計算精度有了很大的提升。但是其設(shè)置邊界長度過大，求解區(qū)域的大小增大了近12倍，計算所需要的時間為240 s，占用了超過20 GB的內(nèi)存，若觀測平面范圍增大，在保證精度的情況下，求解域會進(jìn)一步增大，計算量也會進(jìn)一步增加，嚴(yán)重影響了有限元數(shù)值模擬的正演效率，大大增加了運算成本。而基于有限元—無限元算法的平均相對誤差僅有0.78%，是使用擴(kuò)邊算法精度的10倍，占用內(nèi)存為4.85 GB，僅為擴(kuò)邊算法的 1/5，而運算時間則提高了2.4倍。所以，綜合考慮剩磁和退磁的球體模型實驗結(jié)果說明：基于有限元—無限元的算法，在相同的計算區(qū)域時，邊界處能夠達(dá)到更高的精度；同時也說明在相同的精度下，基于有限元—無限元算法與傳統(tǒng)有限元算法相比能夠在相對更小的區(qū)域內(nèi)達(dá)到精度要求，使用的節(jié)點數(shù)更少，占用內(nèi)存更少，大大降低了運算成本，提高了正演模擬的精度。

表1 不同算法模型計算效率及誤差對比Table 1 Comparison about efficiency and relative error of different model

2.2 組合體模型

圖7為多個地質(zhì)體賦存的較復(fù)雜組合體模型，地磁場大小T0=50 000 nT，地磁傾角I0=45°，偏角A0=0°；球體中心坐標(biāo)為(60，0，-50)，半徑為30 m，剩余磁化強(qiáng)度大小為3.98 A/m，剩磁傾角Ir1=50°，偏角Ar1=60°，球體模型磁化率k1=0.4 SI；圓柱體沿y方向無限延伸，埋深為50 m，截面半徑為15 m，剩余磁化強(qiáng)度大小為7.96 A/m，剩磁傾角Ir2=70°，偏角Ar2=50°，磁化率k2=0.5 SI；棱柱體中心坐標(biāo)為(-60,0,-60)，長100 m，寬40 m，高40 m，方位為西偏北10°，磁化率k3=0.1 SI，剩余磁化強(qiáng)度為0。球體的退磁系數(shù)N1=1/3，而對于任意橫截面形狀的無限長柱體，當(dāng)與長度垂直方向磁化時，退磁系數(shù)N2=1/2[32]，對球體和無限延伸圓柱體進(jìn)行退磁改正，棱柱體磁化率較小且沒有剩磁，則不進(jìn)行退磁改正。

圖7所示，模型求解域邊界長度為200 m，觀測平面設(shè)置在z=0平面，為200 m×200 m的方形區(qū)域，x和y的坐標(biāo)范圍均為-100～100 m。測線位于觀測平面上，測點分布在y=-100～100 m，x=0的直線上。在相同的網(wǎng)格剖分情況下，分別使用基于有限元—無限元的算法和傳統(tǒng)有限單元法進(jìn)行計算。不同有限元算法數(shù)值模擬結(jié)果及平面誤差分布如圖8所示，計算效率和測線中的誤差對比如表2所示。

圖8中，兩種有限單元法的數(shù)值模擬結(jié)果都能較好地展示z=0平面的ΔT磁異常形態(tài)特征，但是數(shù)值模擬的結(jié)果有較大的差異。傳統(tǒng)三維有限元數(shù)值模擬共受到6個方向截斷邊界的影響，當(dāng)求解域邊界長度為600 m時，即求解域邊界長度為觀測區(qū)域邊長的3倍，相對于初始模型求解范圍增大了26倍，使平均相對誤差降低到7.62%，但是其內(nèi)存占用增加了近7倍，計算耗時增加了近9倍，圖8d中，在z=0平面上仍存在近60 nT的絕對誤差，說明傳統(tǒng)算法提升計算精度能力差。而基于有限元—無限元算法，在邊界處的誤差為10 nT左右，平均相對誤差為1.68%，占用內(nèi)存僅為3.67 G。表明基于有限元—無限元算法能在占用更少運算資源的情況下，達(dá)到遠(yuǎn)高于傳統(tǒng)有限單元法的計算精度，更適用于復(fù)雜模型體的三維磁場正演計算，進(jìn)一步說明了基于有限元—無限元算法的優(yōu)越性。

圖7 組合形體模型及網(wǎng)格剖分示意(紅色直線代表測線，黃色為觀測平面范圍)Fig.7 Diagram of combined model and mesh generation(the red line is observation line, the yellow plane is the range of observation plane)

表2 組合形體計算效率及誤差對比Table 2 Comparison about efficiency and relative error of different model

圖8 有限元數(shù)值模擬結(jié)果及平面誤差分布Fig.8 The map of the Finite element numerical simulation results and plane error distribution

在組合體模型的基礎(chǔ)上，將觀測面設(shè)置為起伏地表，如圖9所示，利用有限元—無限元方法進(jìn)行正演數(shù)值模擬，網(wǎng)格剖分參數(shù)和表2中一致。

圖10為模擬結(jié)果，圖10b中有限元—無限元耦合算法的結(jié)果與圖10a中的解析解的異常形態(tài)及幅值基本一致。數(shù)值模擬產(chǎn)生的最大絕對誤差為49.87 nT，平均相對誤差為3.04%，相對于觀測面為z=0平面時誤差有所增加，但仍處于較低的水平。說明基于有限元—無限元耦合算法適用于起伏地形情況，表明了有限元—無限元耦合算法的靈活性。

3 結(jié)論及討論

1) 在進(jìn)行磁場三維正演模擬時，有限單元法會受到6個方向的截斷邊界效應(yīng)的影響，傳統(tǒng)算法在擴(kuò)大邊界時會成倍的增加求解區(qū)域，占用大量運算資源，而計算精度的提升有限。在進(jìn)行復(fù)雜條件下、復(fù)雜形體的正演模擬時，應(yīng)該充分使用基于有限元—無限元的算法，以充分利用計算資源，提高正演效率。

圖9 起伏地表模型(黃色為觀測面范圍)Fig.9 Diagram of rugged surface model (The yellow plane is the range of observation plane)

圖10 起伏地形下有限元數(shù)值模擬結(jié)果及平面誤差分布Fig.10 The map of the Finite element numerical simulation results and plane error distribution in rugged surface

2) 通過孤立球體模型和組合模型的數(shù)值模擬，在相同測區(qū)的計算范圍內(nèi)，有限元—無限元算法與傳統(tǒng)的有限元算法相比，能夠達(dá)到更高的計算精度。尤其是距離邊界較近時，傳統(tǒng)有限元方法誤差較大，而在本文設(shè)置的模型中，基于有限元—無限元的算法能夠達(dá)到2%以內(nèi)的相對誤差。

3) 傳統(tǒng)的有限元方法為達(dá)到較高的精度，需要進(jìn)行擴(kuò)邊處理，設(shè)置較大的求解區(qū)域范圍，且隨著求解域范圍的擴(kuò)大，精度提升有限。而基于有限元—無限元算法，只需與計算區(qū)域相當(dāng)?shù)木W(wǎng)格剖分范圍就可滿足精度要求。相對于傳統(tǒng)的有限元算法，有限元—無限元算法能在占用更少運算資源的情況下，達(dá)到遠(yuǎn)高于傳統(tǒng)有限單元法的計算精度，正演效率更高。

4) 在起伏地形下，基于有限元—無限元耦合算法仍能高精度的完成數(shù)值模擬，進(jìn)一步說明了有限元—無限元算法的靈活性，同時也表明該算法在復(fù)雜地質(zhì)條件下有較好的應(yīng)用前景。

致謝：感謝中國地質(zhì)大學(xué)(武漢)馬火林副教授、朱丹博士后對本文提出的指導(dǎo)意見。