尚曉榮 岳明鑫* 楊曉冬 吳小平 李 勇

(①中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院，安徽合肥 230026； ②中國科學(xué)技術(shù)大學(xué)安徽蒙城地球物理國家野外科學(xué)觀測研究站，安徽蒙城 233500； ③中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所自然資源部地球物理電磁法探測技術(shù)重點實驗室，河北廊坊 065000)

0 引言

可控源電磁法(Controlled-source electromagnetic method)采用人工場源激發(fā)交流電磁波信號(10-3～105Hz)，通過測量目標區(qū)交變電磁場達到研究地下介質(zhì)電性結(jié)構(gòu)的目的。由于成本低、勘探深度大、抗干擾能力強等特點，可控源電磁法在油氣和礦產(chǎn)資源勘探[1-2]、地?zé)豳Y源開發(fā)[3-4]及環(huán)境工程[5-6]等領(lǐng)域發(fā)揮了重要作用?？煽卦措姶欧ò〞r間域電磁法和頻率域電磁法。時間域電磁法利用不接地回線源或接地線源向地下發(fā)射一次脈沖磁場，在一次場間歇期間，利用線圈或接地電極觀測二次渦流場[7-8]。頻率域電磁法常見的是可控源音頻大地電磁法，主要使用接地水平電偶源作為人工信號發(fā)射源[9]，本文研究的即是頻率域可控源電磁法。近年來，為適應(yīng)“深地”國家重大探測戰(zhàn)略中電磁精細結(jié)構(gòu)探測的需求，可控源電磁儀器的研發(fā)及數(shù)據(jù)采集技術(shù)獲得了迅猛發(fā)展，數(shù)據(jù)的反演解釋正逐步邁入三維實用化階段[10-13]。與實際需求相比，三維可控源電磁數(shù)據(jù)的反演技術(shù)尚未成熟，計算效率和準確性均有待提高，特別是考慮起伏地形的三維反演尚不多見[14-15]。正演是反演的基礎(chǔ)，反演解釋的準確性和計算效率很大程度上依賴于正演算法，因此研究可控源電磁數(shù)據(jù)的三維高精度正演方法十分必要。

三維可控源電磁法常用的數(shù)值模擬方法主要包括積分方程法[16-18]、有限差分法[19-20]、有限體積法[21-22]以及有限單元法[23-26]等。任政勇等[27]基于有限元方法提出一種新的非結(jié)構(gòu)化網(wǎng)格局部節(jié)點加密方法，實現(xiàn)了復(fù)雜模型完全非結(jié)構(gòu)化四面體全自動剖分； Puzyrev等[28]基于節(jié)點有限元方法，采用復(fù)雜迭代求解器和基于稀疏近似逆的有效預(yù)調(diào)節(jié)器對大型稀疏線性方程組進行求解，實現(xiàn)了各向異性介質(zhì)中三維可控源電磁場的正演模擬； 蔡紅柱等[29]提出一種基于非結(jié)構(gòu)化網(wǎng)格的海洋電磁有限單元正演算法，該算法適用于井中電磁、航空電磁、環(huán)境地球物理等非均勻且各向異性介質(zhì)的電磁感應(yīng)基礎(chǔ)研究； 湯文武等[30]比較了基于節(jié)點有限元與矢量有限元的三維可控源電磁正演，結(jié)果表明相同條件下矢量有限元法的計算精度更高，但速度較慢； He等[31]基于棱邊有限元方法實現(xiàn)了各向異性介質(zhì)中三維張量可控源電磁法的正演。

常規(guī)可控源電磁數(shù)據(jù)解釋多是基于平坦地形，這會導(dǎo)致反演異常體的位置與形態(tài)常發(fā)生畸變，甚至產(chǎn)生虛假異常。隨著近年來地質(zhì)勘查目標和野外實地勘測場景越來越復(fù)雜，學(xué)者們開始研究起伏地形下的三維電磁正、反演。Nam等[32]采用棱邊有限元方法研究了大地電磁法在三維地形下的響應(yīng)特征； Lin等[33]分析了可控源音頻大地電磁測深數(shù)據(jù)的三維地形效應(yīng)； 殷長春等[34]采用自適應(yīng)非結(jié)構(gòu)有限元法進行了帶地形的時間域海洋電磁法正演。積分方程法、有限差分法和有限體積法大多是基于結(jié)構(gòu)化網(wǎng)格實現(xiàn)的，其中有限單元法對網(wǎng)格剖分比較靈活，非結(jié)構(gòu)化網(wǎng)格更適于模擬真實地形起伏[35]?；诜墙Y(jié)構(gòu)化網(wǎng)格有限元法的二維和三維海洋電磁問題已經(jīng)取得了一些研究成果[36-38]。相比之下，起伏地形下的陸地三維可控源電磁數(shù)值模擬研究相較滯后，目前關(guān)于起伏地形對各個場分量的影響特征的研究尚不多見。

因此，本文基于矢量有限元算法開展了起伏地形下的三維可控源電磁法多分量數(shù)值模擬研究。采用非結(jié)構(gòu)化四面體網(wǎng)格對計算區(qū)域進行剖分，以準確模擬起伏地形，并基于棱邊基函數(shù)實現(xiàn)電磁場的插值和求解。首先，構(gòu)建層狀介質(zhì)模型，將其三維矢量有限元計算結(jié)果與一維解析解[39]進行對比，驗證算法的正確性；其次，構(gòu)建簡單異常體模型，驗證算法對探測低阻異常體的有效性； 然后，分別討論了測點處為山脊地形和山谷地形、發(fā)射源處于山脊地形、發(fā)射源與測點間為山脊地形等情況下頻率域可控源電磁法的響應(yīng)特征； 最后，研究了復(fù)雜三維地質(zhì)模型的電磁響應(yīng)特征，驗證該算法的實用性。

1 數(shù)值模擬理論

從麥克斯韋方程組出發(fā)，假設(shè)時諧因子為e-iω t的平面電磁波入射到均勻介質(zhì)中，則計算區(qū)域電磁場滿足

?×E=iωμH

(1)

?×H=σE-iωεE+J

(2)

式中：E為電場強度；H為磁場強度；ω為角頻率；μ為介質(zhì)磁導(dǎo)率；σ為電導(dǎo)率；ε為相對介電常數(shù)；J為外加電流密度。

由式(1)和式(2)可得整個計算區(qū)域內(nèi)可控源電磁法滿足的電場波動方程

(3)

利用Garlerkin方法，有

(4)

式中：k2=εμω2+iσμω；V表示積分空間。

對式(4)采用基于棱邊的非結(jié)構(gòu)四面體矢量有限元法進行求解，電場矢量場可表示為六個矢量棱邊插值基函數(shù)的線性組合。非結(jié)構(gòu)四面體單元如圖1所示。

圖1 非結(jié)構(gòu)四面體單元e示意圖e1～e6為棱邊編號

圖1中每一個四面體單元的電場強度可表示為[39]

(5)

由式(4)和式(5)并應(yīng)用Garlerkin方法，可得電場強度關(guān)于單元e的有限元方程

AeEe-iωWeEe-ω2CeEe=iωSe

(6)

式中：Ae、We、Ce是單元e的剛度矩陣；Se為單元e的質(zhì)量矩陣。矩陣Ae、We、Ce的第(n，j)個元素和矩陣Se的第n個元素可分別表示為

(7)

(8)

(9)

(10)

式中Je為四面體單元e內(nèi)電偶極子的電流密度。為計算總場方法下的單元源項矩陣Se，將水平線源置于四面體網(wǎng)格的棱邊上，假設(shè)線源的中點坐標為(x0,y0,z0)，則單元e內(nèi)的電流密度Je可表示為

Je=Iδ(x-x0)δ(y-y0)δ(z-z0)

(11)

式中：I表示穿過四面體單元的源電流強度； δ為Dirac-δ函數(shù)。

對起伏地形模型的正演方程組(式(3))求解時，右端項J不僅包含地下電阻率異常體引起的散射電流密度，還包括起伏地形產(chǎn)生的散射電流密度。使用非結(jié)構(gòu)有限元法剖分模型，同時對地形起伏區(qū)域進行網(wǎng)格加密，以更好地對地形進行擬合。

將單元矩陣進行組裝，可得到大型稀疏復(fù)線性方程組

KE=S

(12)

式中：K=Ae-iωWe-ω2Ce，為正演系數(shù)矩陣；S=iωSe，反映源的作用。這里的E即為待求電場強度矢量。通常，K是一大型對稱稀疏矩陣，采用CSR(Compressed Sparse Row，稀疏矩陣按行壓縮)格式進行存儲以減小存儲空間。

最后，使用直接求解器Pardiso求解控制方程，可得到單元網(wǎng)格棱邊上的電場強度和磁感應(yīng)強度分量。

2 數(shù)值模擬算例

2.1 層狀模型

為驗證算法的正確性，構(gòu)建一個三層地層的層狀模型。對模型的電、磁響應(yīng)特征進行分析，并根據(jù)模擬結(jié)果與一維解析解的相對誤差(計算公式為(數(shù)值解-解析解)/解析解)分析本文算法的精度。

層狀模型的空間尺寸為100km×100km×100km，如圖2a所示。發(fā)射源位于y=-10km處，測點位于原點，發(fā)射頻率為0.1～5000Hz，共10個頻點：0.1、0.5、1、5、10、50、100、500、1000、5000Hz。用非結(jié)構(gòu)四面體網(wǎng)格剖分生成的最終網(wǎng)格分布見圖2b。采用Key[40]提出的一維可控源電磁正演算法與本文基于非結(jié)構(gòu)矢量有限元算法進行模擬，并對比響應(yīng)曲線，結(jié)果見圖3。由圖可見，本文算法模擬結(jié)果與一維解析解吻合較好，電場的分量Ex實部和虛部的數(shù)值解與解析解的相對誤差的絕對值分別小于1%和2%，精度較高，驗證了本文算法的正確性。

圖2 層狀模型幾何結(jié)構(gòu)(a)和網(wǎng)格剖分示意圖(b)

圖3 層狀模型Ex分量實部和虛部數(shù)值解與解析解對比(左)及其相對誤差(右)

在保證精度的前提下，分別測試了層狀地層模型不同剖分網(wǎng)格單元數(shù)情況下數(shù)值模擬所需要的時間，并統(tǒng)計了數(shù)值解與解析解的誤差，結(jié)果見表1。

從表1可以看出，網(wǎng)格數(shù)增加約1倍時，數(shù)值模擬耗時大約增加2倍； 網(wǎng)格數(shù)大于33萬后，誤差下降較慢。因此，在綜合考慮耗時與精度的情況下，本文選擇33萬網(wǎng)格數(shù)進行后文模型剖分。

表1 層狀模型數(shù)值模擬時間和精度統(tǒng)計

2.2 簡單組合異常體模型

為了分析本文算法對異常體電阻率的敏感性，構(gòu)建兩組簡單組合異常體模型，背景為均勻半空間。模型包含4個相同的低阻(10Ω·m)/高阻(1000Ω·m)異常體，異常體尺寸為500m×500m×200m(圖4)。模型空間為100km×100km×100km； 發(fā)射源點位于y=-3km處； 測量點位于x軸的-1～1km、y軸-1～1km所圍區(qū)域的地面，x、y方向的點距均為200m； 發(fā)射頻率同前文層狀模型。對接收點附近進行了網(wǎng)格加密，在遠離觀測點處，網(wǎng)格剖分比較稀疏，模型最終網(wǎng)格剖分如圖5所示。

圖4 簡單組合異常體模型

圖5 組合異常體模型網(wǎng)格剖分示意圖左：x=0剖面； 右：z=100m剖面。黃色方塊是異常體

分別對低阻異常體和高阻異常體兩個模型進行模擬，結(jié)果見圖6。從圖中可以看出，低阻異常體產(chǎn)生的電場強度異常較高阻異常體明顯，且不同頻率時的相對異常(絕對值)也不同。

圖6 簡單組合異常體模型的Ex振幅曲線(左)及其相對異常曲線(右)

為了進一步分析背景模型、包含高阻異常體和低阻異常體模型的電場特征，繪制了0.5Hz的Ex振幅圖(圖7)。可以看出，異常體的賦存區(qū)域在有、無異常體的情況下差異較大； 對比低阻和高阻異常體的情況可知，頻率域可控源電磁法對于低阻異常體更敏感。

圖7 0.5Hz時組合異常體模型地表Ex分量振幅圖(a)不存在異常體(背景模型)； (b)存在低阻異常體； (c)存在高阻異常體； (d)圖b與圖a的相對異常； (e)圖c與圖a的相對異常

2.3 三維起伏地形模型

2.3.1 測點處山脊地形模型

建立一個含有山脊地形的三維簡單異常體模型(圖8)。山脊頂部長、寬均為800m，底部長、寬均為1600m，高度為400m，其中心點位于(0，0，200m)。假設(shè)大地背景電阻率為100Ω·m，異常體為不規(guī)則形狀(圖8中的藍色區(qū)域)，厚度為200m，長、寬均為500m，電阻率為1Ω·m，中心點位于(0，0，750m)。采用非結(jié)構(gòu)化四面體單元對模型進行網(wǎng)格剖分，對發(fā)射源和接收點處進行網(wǎng)格加密，模型網(wǎng)格剖分結(jié)果見圖8。采用線源，源中點水平位置為(3km，-3km)，發(fā)射頻率同前文層狀模型采用的頻點； 接收點位于地面-1～1km(x)、-1～1km(y)區(qū)域內(nèi)，x、y方向測點間距均為200m。

圖9為圖8所示三個模型在測點(0，0，0)處的電場強度分量(Ex)和磁感應(yīng)強度分量(By)振幅計算結(jié)果。由圖可見，模型M02與模型M03的正演Ex和By振幅曲線較一致，而M01模型與二者相差較大，說明山脊地形會削弱異常體的存在，同時也說明正演時如果忽略山脊地形，會對正演結(jié)果產(chǎn)生較大影響。對比圖9b與圖9d可以發(fā)現(xiàn)，山脊地形對Ex振幅的影響大于對By振幅的影響。

圖8 山脊地形低阻異常體模型網(wǎng)格剖分示意圖(a)水平地形的低阻異常體模型(M01)； (b)含山脊地形的低阻異常體模型(M02)； (c)含山脊地形的無異常體模型(M03)

圖9 山脊地形模型Ex和By振幅正演曲線及相對異常(a)Ex振幅； (b)Ex振幅相對異常； (c)By振幅;(d)By振幅相對異常

圖10為圖8中三個模型在頻率為1Hz時沿x方向(y=0)的正演分量Ex和By的振幅及其相對異常?？梢钥闯?，山脊根部位置的正演Ex和By振幅均高于水平地形下的場分量幅值，山脊頂部位置的正演Ex和By振幅均低于水平地形下的場分量幅值，這是由于在地形高處(即山脊頂部)，供電電流是發(fā)散的，會呈現(xiàn)低電流密度，而在地形相對低處(即山脊根部)，供電電流是收斂的，呈現(xiàn)高電流密度。此現(xiàn)象說明山脊地形會削弱異常體所產(chǎn)生的異常。

圖10 1Hz時山脊地形模型Ex和By振幅曲線(a)Ex振幅； (b)Ex振幅相對異常； (c)By振幅;(d)By振幅相對異常

2.3.2 測點處山谷地形模型

建立一個含有山谷地形的簡單低阻異常體模型(圖11)。山谷頂部長、寬均為1600m，底部長、寬均為800m，高度為400m，其中心點位于(0，0，200m)。異常體參數(shù)與前文的山脊模型低阻體相同。采用與山脊地形模型相同的加密剖分網(wǎng)格，在相同位置布設(shè)發(fā)射源、接收點，采用同樣的發(fā)射頻率。模型網(wǎng)格剖分結(jié)果見圖11。

圖11 山谷地形模型網(wǎng)格剖分示意圖(a)低阻異常體模型(M04)； (b)無異常體模型(M05)

圖12為水平地形的低阻異常體模型(M01)和圖11中的山谷地形模型M04和M05的Ex和By正演振幅曲線及相對異常曲線。由圖可見，M01、M04和M05三個模型的電、磁場強度振幅曲線差異較大； 與山脊地形的By正演振幅曲線(圖9c和圖9d)相比，山谷地形低阻體模型的差異更大，可見頻率域可控源電磁法的電場強度、磁感應(yīng)強度分量受山谷地形的影響更大。

圖12 山谷地形模型Ex和By振幅正演曲線及相對異常(a)Ex振幅； (b)Ex分量振幅相對異常； (c)By振幅; (d)By振幅相對異常

圖13為山谷模型在頻率為1Hz時沿x方向(y=0)的正演Ex和By振幅及振幅相對異常?？梢钥闯觯谏焦葍蓚?cè)Ex和By振幅均低于水平地形下的場分量幅值，在山谷底部則均相反。這一特征與山脊地形的計算結(jié)果相似，這是由于在地形高處，地表下供電電流發(fā)散，呈現(xiàn)低電流密度； 在地形相對低處(即山谷底部)，供電電流收斂，呈現(xiàn)高電流密度。此現(xiàn)象說明山谷地形會增強異常體產(chǎn)生的異常。

圖13 1Hz時山谷地形模型Ex和By振幅沿x方向變化曲線(a)Ex振幅； (b)Ex振幅相對異常； (c)By振幅; (d)By振幅相對異常

2.3.3 發(fā)射源處起伏地形模型

構(gòu)建一個發(fā)射源處為起伏地形的模型(圖14)。源的中點位于起伏地形區(qū)域的中心位置(0，-8000m，-395m)，山脊頂部長、寬均為1000m，底部長、寬均為1600m，高度為400m。其他參數(shù)設(shè)置均與測點處山脊地形模型參數(shù)相同。發(fā)射源處地形平坦模型(M06)和地形起伏模型(M07)的網(wǎng)格剖分結(jié)果見圖14。

圖14 發(fā)射源處地形平坦模型(M06，上)和地形起伏模型(M07，下)網(wǎng)格剖分示意圖

圖15為模型M06和M07在測點(0，0，0)處的電、磁場分量Ex和By振幅隨頻率變化曲線及相對異常。從圖15a可以看出，頻率高于100Hz時，模型M07的Ex振幅高于模型M06； 頻率低于100Hz時，情況相反。這是因為頻率越低，探測深度越大，發(fā)射源處的地形為山脊地形，相對增加了介質(zhì)厚度； 而且，對于高阻異常體，在相同的頻率下，電場衰減較快，振幅衰減也更快。高頻信號的探測深度較小，發(fā)射源處向下傳導(dǎo)的電流在相同深度接收點處的電場信號強度會更大。另外，從圖15b可以看出，兩個模型的Ex相對振幅差較大，最大值達48%。由圖15b和圖15d可見兩個模型的By分量相對異常與Ex分量有相似特征。

圖15 發(fā)射源處地形起伏模型Ex和By振幅正演曲線及相對異常(a)Ex振幅； (b)Ex分量振幅相對異常； (c)By振幅; (d)By振幅相對異常

2.3.4 發(fā)射源與測點間地形起伏模型

建立一個源點與測點間地形起伏模型(圖16)，分析這類地形對可控源電、磁場響應(yīng)的影響特征。起伏地形山脊頂部長、寬均為800m，底部長、寬均為1600m，高度為400m，中心點為(0，-2000m，0)。發(fā)射源中點位于地表(8km，-8km)處。其他參數(shù)設(shè)置同2.3.1節(jié)模型。圖16為發(fā)射源和測點間地形平坦模型(M06)和地形起伏模型(M08)網(wǎng)格剖分圖。

圖16 發(fā)射源與測點間地形平坦模型(M06，上)和地形起伏模(M08，下)網(wǎng)格剖分示意圖

圖17為模型M06和M08在測點(0，0，0)處的電、磁場分量Ex和By振幅曲線及振幅相對異常曲線。從圖17a和圖17b可見，兩個模型的Ex振幅較接近，差異較小。另外，發(fā)射源與起伏地形的距離、地形起伏程度、地形與異常體的距離等因素也會影響電、磁場。圖17c和圖17d是兩個模型By分量振幅曲線和相對異常曲線，可見By分量振幅相對異常曲線與Ex分量有相似規(guī)律。

圖17 發(fā)射源與測點間地形起伏模型Ex和By振幅正演曲線及相對異常(a)Ex振幅； (b)Ex分量振幅相對異常； (c)By振幅; (d)By振幅相對異常

從上述分析可知，由于地形的存在，尤其是發(fā)射源之下起伏地形和測點處的起伏地形，采用頻率域可控源電磁法對地下異常體進行探測時，電、磁場分量均發(fā)生了畸變，出現(xiàn)了虛假異常。因此在實際數(shù)據(jù)的處理解釋中，需要考慮地表地形對可控源電磁法電、磁場的影響。

3 復(fù)雜金屬礦應(yīng)用實例

為了驗證本文算法在實際應(yīng)用中的可靠性，以湘西某銅鉛鋅多金屬礦床的成礦模式為基礎(chǔ)建立圖18所示地質(zhì)模型，所包含的主要巖(礦)石的電阻率信息如表2所示。

圖18 湘西銅鉛鋅多金屬礦床成礦模式圖

表2 湘西銅鉛鋅多金屬礦床主要巖(礦)石電阻率統(tǒng)計

根據(jù)表中巖(礦)石電阻率范圍，取其平均值，給出圖19所示的金屬礦床模型二維電阻率填圖結(jié)果及其非結(jié)構(gòu)四面體模型網(wǎng)格剖分結(jié)果。對原模型四周區(qū)域進行延拓，其中空—地邊界沿原模型邊界兩端分別向兩側(cè)延伸，而延拓區(qū)域的地層電阻率設(shè)為與原模型相鄰區(qū)域的電阻率一致。延拓模型空間范圍為40km×40km×40km，包含四個長度為1km的線源，發(fā)射電流為1A，源的中點坐標分別為(0，-3000m，105m)、(-3000m，0，115m)、(0，3000m，162m)和(3000m，0，115m)(圖20)。模型中異常體分布區(qū)域為-750m≤x≤750m，-750m≤y≤750m，0≤z≤400m，低阻目標體為圖19中填圖顏色為藍紫色的小規(guī)模片狀塊體及青綠色的塊狀異常體。測點沿地形起伏均勻分布，地表范圍為-750m≤x≤750m，-750≤y≤750m，測點間距為100m。發(fā)射頻率為0.1～5000Hz，頻點分布同2.1節(jié)層狀模型。

圖19 湘西銅鉛鋅多金屬礦床模型二維電阻率填圖(上)及其非結(jié)構(gòu)四面體單元剖分局部示意圖(下)

圖20 湘西銅鉛鋅多金屬礦床延拓模型發(fā)射源布置示意圖

分別正演計算有、無異常體情況下，金屬礦床模型中心點在x=0、y=-200m處Ex振幅隨頻率的變化曲線，結(jié)果見圖21a，二者的相對異常見圖21b。由于模型深部電阻率均一，因此低頻時Ex振幅較穩(wěn)定，且相對異常基本保持為常數(shù)； 在高頻段，由于地質(zhì)異常體分布復(fù)雜，圖21b中的相對異常曲線出現(xiàn)了明顯的波動。

圖21 金屬礦模型偏移距x=0、y=-200m處Ex振幅隨頻率變化曲線(a)及相對異常曲線(b)

選取正演計算的頻率1Hz時的Ex和By振幅數(shù)據(jù)，繪制圖22所示的等值線圖?？梢园l(fā)現(xiàn)，Ex和By振幅異常區(qū)域均與低阻地層的平面位置大致吻合(圖22中虛線橢圓區(qū)域)； 同時，Ex和By振幅相對異常圖中相對異常較大的區(qū)域(圖22中實線橢圓區(qū)域)也與實際低阻異常體的賦存區(qū)域存在很強的相關(guān)性。

圖22 金屬礦模型頻率1Hz時正演計算的地表位置Ex(上)和By(下)振幅等值線圖及相對異常圖(a)包含異常體； (b)不包含異常體； (c)圖b與圖a的相對異常

4 結(jié)論

本文基于非結(jié)構(gòu)四面體網(wǎng)格和矢量有限元法實現(xiàn)面向復(fù)雜地質(zhì)模型的可控源電磁法三維多分量數(shù)值模擬。通過對比層狀模型的數(shù)值解與一維解析解，以及組合異常體模型的電、磁場強度響應(yīng)特征分析，驗證了本文算法的正確性和有效性，說明了模型剖分方案的合理性。通過測點處山脊和山谷地形模型、發(fā)射源下山脊地形模型、發(fā)射源和測點間山脊地形模型的數(shù)值模擬，說明了地形起伏對可控源電、磁場會產(chǎn)生顯著的影響。因此，在地形復(fù)雜的區(qū)域開展可控源電磁工作必須考慮地形的影響。實際復(fù)雜帶地形金屬礦模型的計算結(jié)果也進一步驗證了該算法的實用性和可靠性。