多姿態(tài)具磁性未爆彈磁異?？焖僬?/h1>
2018-05-24 09:00:26吳健生蔣婉聰
吳健生， 蔣婉聰
(1.同濟(jì)大學(xué) 海洋與地球科學(xué)學(xué)院，上海 200092；2.同濟(jì)大學(xué) 海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092)
對(duì)于未爆彈(UXO)的磁異常正演,現(xiàn)有研究中多使用磁偶極子場(chǎng)連續(xù)積分模擬有限長(zhǎng)圓柱體，或使用有限元方法中的立方體單元組合模擬.Zalevsky等[1]將被地球磁場(chǎng)磁化的UXO視作磁偶極子，計(jì)算了3種不同高度的南北走向總磁異常梯度數(shù)據(jù).Churchill等[2]使用COMSOL有限元分析軟件計(jì)算了不同形態(tài)下UXO的磁異常，并將求解結(jié)果與解析法橢球體模型計(jì)算結(jié)果進(jìn)行了對(duì)比.Sinex[3]和 Sanchez等[4]采用立方磁體單元構(gòu)造了UXO和非未爆彈(non-UXO)鐵構(gòu)件磁異常數(shù)值分析模型.葛健等[5]使用三度體的圓柱體模型模擬航彈，并對(duì)航彈的模量磁異常ΔT和垂直磁梯度異常進(jìn)行了數(shù)值模擬和比較.郭志馗等[6]基于偶極線磁場(chǎng)，在三維積分的基礎(chǔ)上得到有限長(zhǎng)圓柱體改進(jìn)后的磁場(chǎng).朱慧慧等[7]使用近似立方體單元重構(gòu)地下鐵質(zhì)管線的幾何結(jié)構(gòu)，結(jié)合磁偶極子構(gòu)造法建立管道磁異常正演模型.
采用偶極線磁場(chǎng)模擬UXO磁異常，計(jì)算快速簡(jiǎn)便，但無(wú)法真實(shí)模擬未爆彈的多樣化以及不規(guī)則形狀，并且適合半徑遠(yuǎn)小于埋深情況下的近似;同時(shí)考慮到未爆彈與管狀物類似，其磁性集中在外殼，內(nèi)充物質(zhì)無(wú)磁性或弱磁性，對(duì)于近地表的精細(xì)磁測(cè)，該假設(shè)存在一定的模型誤差和理論誤差.使用有限元方法中的立方體單元在求解磁化率和剩余磁化強(qiáng)度各向異性的磁性體磁場(chǎng)時(shí)有優(yōu)勢(shì).對(duì)于均勻磁化的磁性體，可使用三角面元法只針對(duì)模型外形進(jìn)行三角形面網(wǎng)格剖分計(jì)算，從而減少立方體單元剖分法巨大的計(jì)算量.
計(jì)算效率在正演中非常重要，迫切需要提高.雖然三角面元法相比其他正演方法有優(yōu)勢(shì)，但是使用三角面元法對(duì)目標(biāo)體整個(gè)表面進(jìn)行剖分計(jì)算效率低，而且三角面元法的計(jì)算量會(huì)隨著網(wǎng)格剖分精度的增加而上升.三角面元法正演的速度很大程度上依賴于網(wǎng)格的剖分效率，對(duì)于一個(gè)新的磁性體，不需要重新剖分網(wǎng)格，而是利用一個(gè)已剖分成三角網(wǎng)格的體積元，在空間上插值從而構(gòu)建磁性體.因此，本文提出采用結(jié)合三角面元法和三次樣條插值方法數(shù)值模擬UXO磁異常，可以兼顧模型的精細(xì)剖分和快速計(jì)算.
1 基于三角面元法的磁異常正演方法
1.1 三角面元法正演原理
與極化電介質(zhì)交界面上存在面束縛電荷相似，在磁化介質(zhì)交界面上也存在面磁荷，面磁荷密度
σm=μ0(M1-M2)·n
(1)
式中：σm為磁化介質(zhì)面磁荷密度；μ0為真空磁導(dǎo)率；M1為磁介質(zhì)1的磁化強(qiáng)度;M2為磁介質(zhì)2的磁化強(qiáng)度；n為界面法線方向的單位矢量，方向從磁介質(zhì)1指向磁介質(zhì)2.式(1)表示磁化介質(zhì)交界面上的面磁荷密度等于2種磁介質(zhì)磁化強(qiáng)度的法向分量之差.已知磁化介質(zhì)體磁荷密度ρm，可以用直接積分法求得磁化介質(zhì)的磁標(biāo)勢(shì),如下所示:
(2)
根據(jù)位場(chǎng)關(guān)系式，可得磁場(chǎng)T的面磁荷積分形式為
(3)
式中:t是空間單位矢量.
將Um分別對(duì)x、y、z求導(dǎo)即可求得磁場(chǎng)水平x方向分量Hax，水平y(tǒng)方向分量Hay和垂直方向分量Za的積分式，如下所示：
(4)
式中:Mn為界面法線方向的磁化強(qiáng)度.最后,根據(jù)地磁場(chǎng)總強(qiáng)度的模量差ΔT與Hax、Hay、Za之間的基本關(guān)系式計(jì)算ΔT，如下所示：
ΔT=HaxcosIcosD+HaycosIsinD+ZasinI
(5)
式中：I是地磁傾角；D是地磁偏角.
從面磁荷積分表達(dá)式出發(fā)，對(duì)任意形狀三度體表面進(jìn)行三角面元剖分.對(duì)剖分后每個(gè)三角面元，根據(jù)三角形的3個(gè)頂點(diǎn)求得每個(gè)三角形磁荷面上的法向量n，得到面磁荷密度σ=M·n，然后通過(guò)積分運(yùn)算，可求得每個(gè)三角形磁荷面的磁位、磁場(chǎng)表達(dá)式(具體推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[8-9])，累加起來(lái)就可以得到磁性體的磁場(chǎng)表達(dá)式.由于磁性體表面是連續(xù)變化的，觀測(cè)點(diǎn)越密集，劃分重建后的模型越精確.
1.2 以球體模型為基礎(chǔ)的正演實(shí)驗(yàn)
球體模型是最基本的、具代表性的有限大小的地質(zhì)體模型，均勻磁化球體的外部磁場(chǎng)與位于球心的偶極子磁場(chǎng)等效，因而球體磁場(chǎng)表達(dá)式有解析公式.對(duì)球體模型進(jìn)行三角剖分，計(jì)算其產(chǎn)生的磁異常，對(duì)比解析解以驗(yàn)證三角面元法的可行性，并對(duì)不同剖分形式下的三角面元計(jì)算結(jié)果進(jìn)行分析對(duì)比.
1.2.1球體模型的剖分方法
圖1 球體三角面元剖分示意圖
(6)
在對(duì)每層薄片的表面進(jìn)行擬合時(shí)，不同三角形剖分方式對(duì)球體表面的擬合程度不一樣，進(jìn)而計(jì)算出的球體磁異常的精度也不同.根據(jù)在每個(gè)薄片與球體相交的圓周上設(shè)立的P個(gè)節(jié)點(diǎn)組成的多邊形是圓外切多邊形還是內(nèi)接多邊形，可將三角面元分為外切面元和內(nèi)接面元.圖2a為一個(gè)球體被切割成9個(gè)薄片，在每個(gè)薄片與球體相交的圓周上設(shè)立的12個(gè)節(jié)點(diǎn)形成的外切三角面元在xOy面上的投影，圖2b則為球體內(nèi)接面元在xOy面上的投影.
a 外切三角面元
b 內(nèi)接三角面元
Fig.2Projectionschematicdiagramofspheredivisionintriangularsurfaceelement
設(shè)球體模型的鉛垂線與z軸重合，z軸向下為正方向.球體的上頂點(diǎn)埋深為h,球體半徑為R.球體模型的上頂點(diǎn)坐標(biāo)為(0,0,h)，下頂點(diǎn)坐標(biāo)為(0,0,h+2R).對(duì)于外切三角面元的剖分方式，中間每一層各個(gè)節(jié)點(diǎn)的坐標(biāo)可以表示為
zij=h+(i-1)Δd
(7)
對(duì)于內(nèi)接三角面元的剖分方式，中間每一層各個(gè)節(jié)點(diǎn)的坐標(biāo)可以表示為
zij=h+(i-1)Δd
(8)
其中,i=2，…,N+1,j=1,…,P+1，Δθ=2π/P，Δd=2R/N.
1.2.2球體模型的正演結(jié)果
通過(guò)球體模型對(duì)三角面元法進(jìn)行驗(yàn)證.僅考慮磁性體的感應(yīng)磁化強(qiáng)度，并且在均勻磁化的情況下，設(shè)地磁場(chǎng)強(qiáng)度為49 600 nT,磁傾角為44.5°，磁偏角為-2.5°，磁化率為2.設(shè)球體半徑為1 m，中心點(diǎn)埋深為2 m，中心位于坐標(biāo)原點(diǎn)處.圖3所示為N=24、P=12情況下，球體模型產(chǎn)生的ΔT磁異常剖面曲線.從圖3可見(jiàn)，使用外切三角面元計(jì)算得到的球體磁異常曲線比使用內(nèi)接三角面元計(jì)算得到的球體磁異常曲線更接近理論磁異常曲線，與理論曲線重合.外切和內(nèi)接三角面元正演得到球體磁異常的均方根誤差(RMSE)分別為0.22和0.83，這表明不同節(jié)點(diǎn)位置構(gòu)建的外切或內(nèi)接多面體模型，會(huì)對(duì)目標(biāo)磁性體磁異常正演產(chǎn)生一定影響.同樣網(wǎng)格剖分下，外切和內(nèi)接三角面元對(duì)物體表面的逼近速度不同，對(duì)球體模型而言，外切三角面元的逼近速度更快.同時(shí)，進(jìn)一步細(xì)化剖分網(wǎng)格，可以進(jìn)一步提高正演精度.
圖3 不同方法球體磁異常剖面曲線對(duì)比
2 不同姿態(tài)UXO模型的正演方法
UXO大體上可以分為4類：集束彈藥、普通陸戰(zhàn)武器彈藥、航空炸彈、地雷[10].UXO外殼多為鐵磁性物質(zhì)，磁性參數(shù)與周圍環(huán)境的差別較大，在一定條件下能夠被磁法勘探儀器所識(shí)別.然而,UXO外形差別較大.航空炸彈一般體形較大，長(zhǎng)度較長(zhǎng)，外形與圓柱體比較接近，因此可使用均勻有限長(zhǎng)的三度體圓柱體模擬航空炸彈；其他幾類彈藥體形較小，體長(zhǎng)較短，長(zhǎng)度和埋深的比值很小，理論和實(shí)踐已經(jīng)證明這幾類彈藥產(chǎn)生的磁場(chǎng)近似于球體所產(chǎn)生的磁場(chǎng).常見(jiàn)的航空炸彈的外形往往不是標(biāo)準(zhǔn)的有限長(zhǎng)圓柱體，而是變徑有限長(zhǎng)圓柱體，即沿長(zhǎng)軸方向半徑不同，用單一有限長(zhǎng)圓柱體的模型無(wú)法很好地近似，但可以用不同半徑的圓柱體元進(jìn)行組合.
2.1 有限長(zhǎng)圓柱體模型的剖分方法
與球體的剖分方式類似，沿圓柱體長(zhǎng)軸方向切割成N個(gè)等厚薄片，在每個(gè)薄片與柱體相交的圓周上等弧長(zhǎng)地設(shè)立P個(gè)節(jié)點(diǎn)，在每層薄片上用P個(gè)三角形來(lái)擬合薄片表面，如圖4所示.對(duì)于圓柱體的頂/底面，頂/底面圓心與第一層各節(jié)點(diǎn)的連線共形成P個(gè)等腰三角形，而對(duì)于沿圓柱體長(zhǎng)軸切割出的N個(gè)薄片，則能組成2P(N-1)個(gè)三角面元.同樣，在對(duì)每層薄片的表面進(jìn)行擬合時(shí)，也可根據(jù)在每個(gè)薄片與柱體相交的圓周上設(shè)立P個(gè)節(jié)點(diǎn)組成的多邊形是圓外切多邊形還是內(nèi)接多邊形，將三角面元分為外切面元和內(nèi)接面元.
除此之外，如果圓柱體各個(gè)薄片上的對(duì)應(yīng)節(jié)點(diǎn)Dij和Di+1,j(i=1,2,…,N,j=1,2,…,P)位于同一條軸線上，那么按照斜線DijDi+1,j+1分割出的三角形則為直角三角形(見(jiàn)圖4a).如果圓柱體各個(gè)薄片按照奇偶層交錯(cuò)選取節(jié)點(diǎn)，則可形成等腰三角形(見(jiàn)圖4b).
a 直角三角面元
b 等腰三角面元
2.2 全空間任意姿態(tài)有限長(zhǎng)圓柱體磁異常計(jì)算方法
利用三角面元法計(jì)算有限長(zhǎng)圓柱體磁場(chǎng)的優(yōu)勢(shì)在于，不通過(guò)旋轉(zhuǎn)坐標(biāo)系，即不改變觀測(cè)點(diǎn)坐標(biāo)，即可快速計(jì)算出有限長(zhǎng)圓柱體在空間任意姿態(tài)(任意水平和垂直位置、任意走向、任意傾角)下的磁場(chǎng)，減少了計(jì)算量，并能保證快速穩(wěn)定計(jì)算.
假設(shè)起始位置處走向0°為x軸正方向(x軸正方向?yàn)榇疟保蛭鳛樨?fù)，向東為正)，傾角從水平投影面起算為0°(向下為正)，中心點(diǎn)空間位置為(0,0,0)，在目標(biāo)位置處設(shè)走向角度為β，傾角為α，中心點(diǎn)空間位置為(x1,y1,z1).對(duì)圓柱體表面剖分出的任意一個(gè)節(jié)點(diǎn)Dij的三維坐標(biāo)(xij,yij,zij)進(jìn)行投影，先投影到傾角方向，再投影到走向方向，最后平移到目標(biāo)位置，投影運(yùn)算實(shí)際上是投影矩陣與坐標(biāo)相乘.走向方向和傾角方向的投影矩陣分別為Tβ和Tα，如下所示:
(9)
最終可以得到目標(biāo)位置處節(jié)點(diǎn)Dij的坐標(biāo)為
(10)
計(jì)算出每個(gè)節(jié)點(diǎn)的坐標(biāo)后，按照一定順序輸出三角形的3個(gè)角點(diǎn)坐標(biāo)，各個(gè)角點(diǎn)的輸出順序應(yīng)統(tǒng)一為順時(shí)針?lè)较颍?jì)算三角面元的平面法向量時(shí)遵循右手螺旋法則.對(duì)每個(gè)三角面元產(chǎn)生的磁異常進(jìn)行計(jì)算，將所有三角面元產(chǎn)生的磁異常在每一個(gè)觀測(cè)點(diǎn)處累加，即可得到該觀測(cè)點(diǎn)處目標(biāo)體產(chǎn)生的磁異常.
2.3 有限長(zhǎng)圓柱體模型的正演結(jié)果
設(shè)磁場(chǎng)參數(shù)同1.2.2節(jié)，僅考慮磁性體的感應(yīng)磁化強(qiáng)度且均勻磁化，模型體各項(xiàng)物理參數(shù)相同，即磁化率均為2,半徑均為0.1 m，長(zhǎng)度均為2 m.圖5是N=20、P=36下，采用等腰內(nèi)接三角面元的計(jì)算結(jié)果.圖5a和5b是中心點(diǎn)平面坐標(biāo)位于(1,1)、中心點(diǎn)埋深1 m的有限長(zhǎng)圓柱體在不同走向β和傾角α下的形態(tài),圖5c和5d是2種形態(tài)下對(duì)應(yīng)的磁異常.
a 圓柱體1(β=45°,α=0°)
b 圓柱體2(β=45°,α=30°)
c 圓柱體1磁場(chǎng)
d 圓柱體2磁場(chǎng)
圖5不同走向和傾角下有限長(zhǎng)圓柱體的形態(tài)及其產(chǎn)生的磁場(chǎng)
Fig.5Shapeoffinite-lengthcylinderindifferentstrikesanddipanglesandthecorrespondingmagneticfield
從圖5可以看出，磁異常平面圖中磁異常的形態(tài)與位置以及大小反映了模型體的空間位置，比較容易標(biāo)定和識(shí)別，進(jìn)而可以大致確定出UXO的走向和傾向.這一點(diǎn)對(duì)于在實(shí)際探測(cè)中根據(jù)探測(cè)到的磁異常形態(tài)對(duì)UXO進(jìn)行空間定位有重要的理論意義.
對(duì)于中心點(diǎn)坐標(biāo)(0,0,2) m的垂直有限長(zhǎng)圓柱體，使用3種剖分方法得到的三角面元分別計(jì)算出磁異常剖面曲線，對(duì)比目前研究中使用偶極線磁場(chǎng)模擬得到的圓柱體磁異常曲線，如圖6所示.其中，對(duì)于三角面元法，N=20，P=36.
圖6 不同方法垂直有限長(zhǎng)圓柱體磁異常剖面曲線對(duì)比
Fig.6Comparisonofsectioncurvesofverticalfinite-lengthcylindricalmagneticanomalyamongdifferentmethods
理論上三角面元的面積和應(yīng)盡可能接近圓柱體的表面積，則計(jì)算出的圓柱體磁異常可以無(wú)限逼近真實(shí)磁異常.若以圓柱體表面積(此例中為13 194.69 cm2)和三角面元面積和之間的差值作為評(píng)判標(biāo)準(zhǔn)，則在相同的層數(shù)和節(jié)點(diǎn)數(shù)目下3種剖分方式面積和如表1所示.
表1圓柱體三角面元面積和與實(shí)際面積差
Tab.1Differencebetweenthecalculatedareausingtriangularsurfaceelementandtheactualareaofthecylinder
三角面元形式三角面元面積和/cm2面積差絕對(duì)值/cm2等腰外切13228.3833.69等腰內(nèi)接13175.6519.04直角內(nèi)接13175.5619.13
從圖6和表1可以看出,在內(nèi)接情況下等腰和直角三角面元對(duì)磁異常正演結(jié)果的影響較小，但相比于等腰內(nèi)接三角面元等腰外切三角面元對(duì)磁異常正演結(jié)果的影響則較大.理論上外切三角面元計(jì)算得到的圓柱體磁異常比真實(shí)值大，而內(nèi)接三角面元計(jì)算出的圓柱體磁異常比真實(shí)值小，但從圖6可見(jiàn)，采用偶極線計(jì)算得到的磁異常誤差比外切三角面元計(jì)算的磁異常更大，且該誤差無(wú)法修正改進(jìn).對(duì)圓柱體模型而言，內(nèi)接三角面元的逼近速度更快，等腰內(nèi)接三角元計(jì)算出的磁異常最接近真實(shí)值.
值得注意的是，現(xiàn)有的UXO磁異常計(jì)算[5-6]認(rèn)為UXO內(nèi)部鐵磁性物質(zhì)是均勻分布的，但實(shí)際未爆物的磁性主要來(lái)源于UXO鐵磁性外殼，外殼以內(nèi)爆炸物往往是弱磁性或無(wú)磁性.因此，要進(jìn)一步提高正演精度，則必須考慮這種內(nèi)部“中空”磁性結(jié)構(gòu),只有基于表面積分法的三角面元法能夠有效解決該問(wèn)題.將未爆物產(chǎn)生的磁異常與管狀磁性物產(chǎn)生的磁異常近似，與相同空間位置上不同半徑有限長(zhǎng)圓柱體在觀測(cè)點(diǎn)處產(chǎn)生的磁場(chǎng)相減，來(lái)近似計(jì)算均勻磁化中空UXO產(chǎn)生的磁異常.
3 基于三次樣條的UXO磁異?？焖僬莘椒?/h2>
雖然使用三角面元法計(jì)算均勻磁化磁性體磁異常精度高，但是計(jì)算不同尺寸的目標(biāo)體時(shí)需要對(duì)目標(biāo)體表面重復(fù)進(jìn)行網(wǎng)格剖分.隨著三角網(wǎng)格的細(xì)化，計(jì)算量會(huì)提高，導(dǎo)致計(jì)算效率下降.因此，引入三次樣條插值方法，對(duì)不同長(zhǎng)度的模型，不需要重新剖分網(wǎng)格，而是利用一個(gè)已剖分成三角網(wǎng)格的體積元在空間上插值構(gòu)建磁性體，進(jìn)而實(shí)現(xiàn)UXO的磁異?？焖僬?UXO模型近似為有限長(zhǎng)圓柱體或變徑有限長(zhǎng)柱體.
正演思想可以分為以下幾個(gè)步驟：
步驟1將UXO模型磁性體沿長(zhǎng)軸方向均勻切割為幾個(gè)塊體，每個(gè)塊體形態(tài)和尺寸相同或相似(不同塊體可以通過(guò)通過(guò)同一塊體半徑變化獲得).對(duì)其中一塊進(jìn)行上述三角面元剖分，但對(duì)切割面不進(jìn)行剖分，因?yàn)榍懈蠲娌粚儆谀繕?biāo)體整體的表面.對(duì)邊界塊體上的端面要單獨(dú)進(jìn)行剖分，并不計(jì)入以下步驟的計(jì)算.
步驟2對(duì)每一個(gè)觀測(cè)點(diǎn)計(jì)算出小塊體在UXO腔體內(nèi)任意幾處以及邊界位置上的磁異常.
步驟3選取步驟2中幾處位置上的磁異常值作為插值節(jié)點(diǎn)，以塊體的切割間距進(jìn)行插值.
步驟4在每一個(gè)觀測(cè)點(diǎn)處，累加插值得到各個(gè)磁異常值，并補(bǔ)充邊界塊體上的端面產(chǎn)生的磁異常.
3.1 三次樣條插值方法基本原理
Lagrange插值、Newton插值及Hermite插值等高次插值中會(huì)出現(xiàn)龍格(Runge)現(xiàn)象.分段插值雖然可以避免龍格現(xiàn)象,但是分段線性插值在節(jié)點(diǎn)處不光滑。分段Hermite插值僅僅要求插值函數(shù)在節(jié)點(diǎn)上等于被插函數(shù)的一階導(dǎo)數(shù)，并且導(dǎo)數(shù)值往往并不已知.三次樣條插值先由函數(shù)值確定導(dǎo)數(shù)值，再由分段Hermite插值解決問(wèn)題，因而避免了龍格現(xiàn)象，并且插值曲線光滑.
三次樣條函數(shù)的定義如下所示[11]：
設(shè)對(duì)任意區(qū)間[a,b]有剖分Δ∶a=l0(1)S(l)∈C2[a,b]，即具有連續(xù)的一階、二階導(dǎo)數(shù)，Sk,v(lk)=Sk-1,v(lk),v=0,1,2，k=0,1,…,n-1.
(2)S(l)在每一個(gè)小區(qū)間[lk,lk+1]上是次數(shù)≤3的多項(xiàng)式，則稱S(l)為關(guān)于剖分Δ的一個(gè)三次樣條函數(shù).
(3) 設(shè)給定函數(shù)表(lq,f(lq)),q=0,1,…,n，若三次樣條函數(shù)S(l)還滿足以下插值條件：
S(lq)=f(lq),q=0,1,…,n
則稱S(l)為f(l)關(guān)于剖分Δ的三次樣條插值函數(shù).求取三次樣條插值函數(shù)表達(dá)式，首先以分段三次Hermite插值為基礎(chǔ)，滿足上述條件的情況下，從下述3種邊界條件中的一種推導(dǎo)而來(lái):
(1) 第一種邊界條件.S′(l0)、S′(ln)已知，即S′(l0)=f′(l0)，S′(ln)=f′(ln).
(2) 第二種邊界條件.S″(l0)、S″(ln)已知，即S″(l0)=f″(l0)，S″(ln)=f″(ln)，若S″(l0)=S″(ln)=0，則稱為自然邊界條件.
(4) 非扭結(jié)邊界.使兩端點(diǎn)處多項(xiàng)式三次項(xiàng)系數(shù)與這兩端點(diǎn)的鄰近點(diǎn)的三次項(xiàng)系數(shù)相等，此時(shí)不用給出邊界值.
3.2 不同模型的正演結(jié)果
對(duì)于半徑R=0.1 m、長(zhǎng)度L=2 m、埋深h=1 m、中心點(diǎn)平面坐標(biāo)為(0,0) m的水平有限長(zhǎng)圓柱體，將其沿走向均勻切割成N=40塊，每一塊都是長(zhǎng)L/N的小圓柱體，將每一塊的中心點(diǎn)橫坐標(biāo)位置作為插值點(diǎn)，在空間上對(duì)應(yīng)有N個(gè)插值點(diǎn)c1,…,cN.對(duì)每一塊體進(jìn)行截面節(jié)點(diǎn)P=12的等腰內(nèi)接三角面元剖分，在圓柱體上的對(duì)稱位置任意選擇n=5個(gè)型值點(diǎn)E1,…,E5，對(duì)每個(gè)測(cè)點(diǎn)用三角面元法計(jì)算出單個(gè)塊體在5個(gè)型值點(diǎn)處的磁場(chǎng)值f(Ek)(k=1,…,n)作為插值點(diǎn).
在每一個(gè)觀測(cè)點(diǎn)上，使用自然邊界條件下的三次樣條函數(shù)，插值出塊體在N個(gè)插值點(diǎn)處的磁場(chǎng)值S(ci)(i=1,…,N)，對(duì)得到的結(jié)果進(jìn)行疊加，即得到有限長(zhǎng)圓柱體的磁場(chǎng).圖7a展示了使用插值和不使用插值計(jì)算出的水平圓柱體磁異常剖面曲線對(duì)比，圖7b展示了兩者的絕對(duì)誤差和相對(duì)誤差分布.絕對(duì)誤差分布與磁異常曲線相似，最大值位于端點(diǎn)附近為11 nT,相對(duì)誤差平均值約2.9%.圖7c展示了地面中心觀測(cè)點(diǎn)處，根據(jù)5個(gè)型值點(diǎn)-0.9、-0.5、0、0.5、0.9 m插值得到的磁場(chǎng)值，可見(jiàn)插值曲線非常光滑.將40個(gè)插值點(diǎn)的磁場(chǎng)值累加，即得到中心觀測(cè)點(diǎn)處水平圓柱體整體的磁場(chǎng)值.
a 三次樣條插值剖面曲線
b 誤差曲線
c 中心觀測(cè)點(diǎn)處插值曲線
圖8 三角面元構(gòu)造的變徑紡錘形UXO模型
Fig.8FusiformUXOmodelwithvariablediametersmadebytriangularsurfaceelement
對(duì)于變徑有限長(zhǎng)圓柱體，逐段計(jì)算每個(gè)圓柱體元產(chǎn)生的磁異常后疊加，計(jì)算量較大，可以用上述原理進(jìn)行快速正演.插值計(jì)算過(guò)程同上述，對(duì)于在變徑有限長(zhǎng)柱體上的對(duì)稱位置任意選擇n=5個(gè)型值點(diǎn)(-0.9、-0.5、0、0.5、0.9 m)，對(duì)每個(gè)測(cè)點(diǎn)計(jì)算出單個(gè)塊體在n個(gè)型值點(diǎn)處的磁場(chǎng)值作為插值點(diǎn).在每一個(gè)觀測(cè)點(diǎn)上，使用自然邊界條件下的三次樣條函數(shù)，插值出塊體在N個(gè)插值點(diǎn)處的磁場(chǎng)值，對(duì)得到的結(jié)果進(jìn)行疊加，即得到變徑有限長(zhǎng)圓柱體的磁場(chǎng)值.圖9a展示了使用插值和不使用插值計(jì)算出的變徑有限長(zhǎng)圓柱體模型磁異常曲線對(duì)比.圖9b展示了兩者的絕對(duì)誤差和相對(duì)誤差分布，絕對(duì)誤差最大值位于端點(diǎn)附近，為12 nT,平均值為0.7 nT,相對(duì)誤差平均值為0.35%,最大值為8%.
4 結(jié)語(yǔ)
本文針對(duì)工程探測(cè)中遇到的地下多姿態(tài)和形態(tài)各異的均勻磁化UXO的磁正演問(wèn)題，基于面磁荷積分原理，采用三角面元法的和三次樣條插值方法，實(shí)現(xiàn)了高精度快速正演模擬，并進(jìn)行了相關(guān)問(wèn)題的分析.
本文提出的正演計(jì)算方法，有較好的計(jì)算速度和正演精度，可以模擬不規(guī)則外形以及任意空間位置下空心UXO的磁異常.該方法可進(jìn)一步與反演算法結(jié)合，以提高UXO定位的速度和準(zhǔn)確性，也可應(yīng)用于地下復(fù)雜金屬管道的磁正演.
a 三次樣條插值剖面曲線
b 誤差曲線
Fig.9Resultsofvariablediametercylindricalmagneticanomalycalculatedwithsplineinterpolation
參考文獻(xiàn):
[1] ZALEVSKY Z,BREGMAN Y,SALOMONSKI N,etal.Resolution enhanced magnetic sensing system for wide coverage real time UXO detection[J].Journal of Applied Geophysics,2012,84(4):70.
[2] CHURCHILL K M,LINK C,YOUMANS C C.A comparison of the finite-element method and analytical method for modeling unexploded ordnance using magnetometer[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(7):2720.
[3] SINEX D B.Advancing the state of the art of UXO discrimination for total-field magnetic data[D].Golden:Colorado School of Mines,2006.
[4] SANCHEZ V,LI Y G,NABIGHIAN M N,etal.Numerical modeling of higher order magnetic moments in UXO discrimination[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(9):2568.
[5] 葛健，陸承達(dá)，董浩斌，等. 基于Overhauser傳感器的近地表UXO磁梯度法探測(cè)技術(shù)[J]. 儀器儀表學(xué)報(bào)，2015，36(5)：961.
GE Jian,LU Chengda,DONG Haobin,etal.The detection technology of near-surface UXO based on magnetic gradient method and Overhauser sensor[J].Chinese Journal of Scientific Instrument,2015，36(5)：961.
[6] 郭志馗，陳超，陶春輝，等. 有限長(zhǎng)圓柱體磁異常場(chǎng)全空間正演方法[J]. 地球物理學(xué)報(bào), 2017,60(4):1557.
GUO Zhikui,CHEN Chao,TAO Chunhui,etal.Forward modeling of magnetic anomaly of finite-length cylinder in 3D space[J].Chinese Journal of Geophysics,2017,60(4):1557.
[7] 朱慧慧，劉得軍，馮碩,等. 基于磁偶極子構(gòu)造法幾何建模的鐵管磁異常正演[J]. 計(jì)算機(jī)測(cè)量與控制,2017,25(3):201.
ZHU Huihui,LIU Dejun,FENG Shuo,etal.Magnetic dipole reconstruction geometry modeling for underground ferromagnetic pipe magnetic abnormal detection[J].Computer Measurement and Control,2017,25(3):201.
[8] 吳文鸝, 管志寧, 高艷芳,等. 重磁異常數(shù)據(jù)三維人機(jī)聯(lián)作模擬[J]. 物探化探計(jì)算技術(shù), 2005, 27(3):227.
WU Wenli,GUAN Zhining,GAO Yanfang,etal.Interactive man/computer modeling 3D body of gravity and magnetic anomaly data[J].Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(3):227.
[9] 管志寧.地磁場(chǎng)與磁力勘探[M].北京：地質(zhì)出版社，2005.
GUAN Zhining.Geomagnetic field and magnetic exploration[M].Beijing:Geological Publishing House,2005.
[10] 李小康，楊磊，劉磊.未爆彈藥問(wèn)題及地球物理解決方案綜述[J].中國(guó)礦業(yè)，2010，19:187.
LI Xiaokang,YANG Lei,LIU Lei.Overview on unexploded ordnance problem and a solution: the geophysical scheme[J].China Mining Magzine,2010，19:187.
[11] 關(guān)冶，陸金甫.數(shù)值方法[M].北京：清華大學(xué)出版社,2011.
GUAN Ye,LU Jinfu.Numerical method[M].Beijing:Tsinghua University House,2011.
猜你喜歡
隨機(jī)粗糙面散射中遮蔽效應(yīng)算法的改進(jìn)遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版)(2022年4期)2022-12-27 12:27:28 關(guān)于橢圓外切平行四邊形的一個(gè)幾何不變量中等數(shù)學(xué)(2022年1期)2022-06-05 07:50:08 基于重心剖分的間斷有限體積元方法南京大學(xué)學(xué)報(bào)(數(shù)學(xué)半年刊)(2020年1期)2020-03-19 02:24:48 二元樣條函數(shù)空間的維數(shù)研究進(jìn)展赤峰學(xué)院學(xué)報(bào)·自然科學(xué)版(2018年9期)2018-10-18 07:15:52 探究拋物線內(nèi)接、外切三角形的性質(zhì)中等數(shù)學(xué)(2018年4期)2018-08-01 06:36:34 橢圓內(nèi)接外切六邊形的幾何特性研討中學(xué)數(shù)學(xué)研究(廣東)(2018年23期)2018-03-05 07:54:32 圓外切三角形與圓的關(guān)系初中生世界·九年級(jí)(2017年9期)2017-10-13 10:27:04 基于改進(jìn)Gordon方程的RCS快速算法裝備學(xué)院學(xué)報(bào)(2016年5期)2016-11-14 06:26:42 一種實(shí)時(shí)的三角剖分算法池州學(xué)院學(xué)報(bào)(2015年3期)2016-01-05 01:13:07 復(fù)雜地電模型的非結(jié)構(gòu)多重網(wǎng)格剖分算法地震地質(zhì)(2015年3期)2015-12-25 03:29:42
吳健生， 蔣婉聰
(1.同濟(jì)大學(xué) 海洋與地球科學(xué)學(xué)院，上海 200092；2.同濟(jì)大學(xué) 海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092)
對(duì)于未爆彈(UXO)的磁異常正演,現(xiàn)有研究中多使用磁偶極子場(chǎng)連續(xù)積分模擬有限長(zhǎng)圓柱體，或使用有限元方法中的立方體單元組合模擬.Zalevsky等[1]將被地球磁場(chǎng)磁化的UXO視作磁偶極子，計(jì)算了3種不同高度的南北走向總磁異常梯度數(shù)據(jù).Churchill等[2]使用COMSOL有限元分析軟件計(jì)算了不同形態(tài)下UXO的磁異常，并將求解結(jié)果與解析法橢球體模型計(jì)算結(jié)果進(jìn)行了對(duì)比.Sinex[3]和 Sanchez等[4]采用立方磁體單元構(gòu)造了UXO和非未爆彈(non-UXO)鐵構(gòu)件磁異常數(shù)值分析模型.葛健等[5]使用三度體的圓柱體模型模擬航彈，并對(duì)航彈的模量磁異常ΔT和垂直磁梯度異常進(jìn)行了數(shù)值模擬和比較.郭志馗等[6]基于偶極線磁場(chǎng)，在三維積分的基礎(chǔ)上得到有限長(zhǎng)圓柱體改進(jìn)后的磁場(chǎng).朱慧慧等[7]使用近似立方體單元重構(gòu)地下鐵質(zhì)管線的幾何結(jié)構(gòu)，結(jié)合磁偶極子構(gòu)造法建立管道磁異常正演模型.
采用偶極線磁場(chǎng)模擬UXO磁異常，計(jì)算快速簡(jiǎn)便，但無(wú)法真實(shí)模擬未爆彈的多樣化以及不規(guī)則形狀，并且適合半徑遠(yuǎn)小于埋深情況下的近似;同時(shí)考慮到未爆彈與管狀物類似，其磁性集中在外殼，內(nèi)充物質(zhì)無(wú)磁性或弱磁性，對(duì)于近地表的精細(xì)磁測(cè)，該假設(shè)存在一定的模型誤差和理論誤差.使用有限元方法中的立方體單元在求解磁化率和剩余磁化強(qiáng)度各向異性的磁性體磁場(chǎng)時(shí)有優(yōu)勢(shì).對(duì)于均勻磁化的磁性體，可使用三角面元法只針對(duì)模型外形進(jìn)行三角形面網(wǎng)格剖分計(jì)算，從而減少立方體單元剖分法巨大的計(jì)算量.
計(jì)算效率在正演中非常重要，迫切需要提高.雖然三角面元法相比其他正演方法有優(yōu)勢(shì)，但是使用三角面元法對(duì)目標(biāo)體整個(gè)表面進(jìn)行剖分計(jì)算效率低，而且三角面元法的計(jì)算量會(huì)隨著網(wǎng)格剖分精度的增加而上升.三角面元法正演的速度很大程度上依賴于網(wǎng)格的剖分效率，對(duì)于一個(gè)新的磁性體，不需要重新剖分網(wǎng)格，而是利用一個(gè)已剖分成三角網(wǎng)格的體積元，在空間上插值從而構(gòu)建磁性體.因此，本文提出采用結(jié)合三角面元法和三次樣條插值方法數(shù)值模擬UXO磁異常，可以兼顧模型的精細(xì)剖分和快速計(jì)算.
1 基于三角面元法的磁異常正演方法
1.1 三角面元法正演原理
與極化電介質(zhì)交界面上存在面束縛電荷相似，在磁化介質(zhì)交界面上也存在面磁荷，面磁荷密度
σm=μ0(M1-M2)·n
(1)
式中：σm為磁化介質(zhì)面磁荷密度；μ0為真空磁導(dǎo)率；M1為磁介質(zhì)1的磁化強(qiáng)度;M2為磁介質(zhì)2的磁化強(qiáng)度；n為界面法線方向的單位矢量，方向從磁介質(zhì)1指向磁介質(zhì)2.式(1)表示磁化介質(zhì)交界面上的面磁荷密度等于2種磁介質(zhì)磁化強(qiáng)度的法向分量之差.已知磁化介質(zhì)體磁荷密度ρm，可以用直接積分法求得磁化介質(zhì)的磁標(biāo)勢(shì),如下所示:
(2)
根據(jù)位場(chǎng)關(guān)系式，可得磁場(chǎng)T的面磁荷積分形式為
(3)
式中:t是空間單位矢量.
將Um分別對(duì)x、y、z求導(dǎo)即可求得磁場(chǎng)水平x方向分量Hax，水平y(tǒng)方向分量Hay和垂直方向分量Za的積分式，如下所示：
(4)
式中:Mn為界面法線方向的磁化強(qiáng)度.最后,根據(jù)地磁場(chǎng)總強(qiáng)度的模量差ΔT與Hax、Hay、Za之間的基本關(guān)系式計(jì)算ΔT，如下所示：
ΔT=HaxcosIcosD+HaycosIsinD+ZasinI
(5)
式中：I是地磁傾角；D是地磁偏角.
從面磁荷積分表達(dá)式出發(fā)，對(duì)任意形狀三度體表面進(jìn)行三角面元剖分.對(duì)剖分后每個(gè)三角面元，根據(jù)三角形的3個(gè)頂點(diǎn)求得每個(gè)三角形磁荷面上的法向量n，得到面磁荷密度σ=M·n，然后通過(guò)積分運(yùn)算，可求得每個(gè)三角形磁荷面的磁位、磁場(chǎng)表達(dá)式(具體推導(dǎo)過(guò)程見(jiàn)文獻(xiàn)[8-9])，累加起來(lái)就可以得到磁性體的磁場(chǎng)表達(dá)式.由于磁性體表面是連續(xù)變化的，觀測(cè)點(diǎn)越密集，劃分重建后的模型越精確.
1.2 以球體模型為基礎(chǔ)的正演實(shí)驗(yàn)
球體模型是最基本的、具代表性的有限大小的地質(zhì)體模型，均勻磁化球體的外部磁場(chǎng)與位于球心的偶極子磁場(chǎng)等效，因而球體磁場(chǎng)表達(dá)式有解析公式.對(duì)球體模型進(jìn)行三角剖分，計(jì)算其產(chǎn)生的磁異常，對(duì)比解析解以驗(yàn)證三角面元法的可行性，并對(duì)不同剖分形式下的三角面元計(jì)算結(jié)果進(jìn)行分析對(duì)比.
1.2.1球體模型的剖分方法
圖1 球體三角面元剖分示意圖
(6)
在對(duì)每層薄片的表面進(jìn)行擬合時(shí)，不同三角形剖分方式對(duì)球體表面的擬合程度不一樣，進(jìn)而計(jì)算出的球體磁異常的精度也不同.根據(jù)在每個(gè)薄片與球體相交的圓周上設(shè)立的P個(gè)節(jié)點(diǎn)組成的多邊形是圓外切多邊形還是內(nèi)接多邊形，可將三角面元分為外切面元和內(nèi)接面元.圖2a為一個(gè)球體被切割成9個(gè)薄片，在每個(gè)薄片與球體相交的圓周上設(shè)立的12個(gè)節(jié)點(diǎn)形成的外切三角面元在xOy面上的投影，圖2b則為球體內(nèi)接面元在xOy面上的投影.
a 外切三角面元
b 內(nèi)接三角面元
Fig.2Projectionschematicdiagramofspheredivisionintriangularsurfaceelement
設(shè)球體模型的鉛垂線與z軸重合，z軸向下為正方向.球體的上頂點(diǎn)埋深為h,球體半徑為R.球體模型的上頂點(diǎn)坐標(biāo)為(0,0,h)，下頂點(diǎn)坐標(biāo)為(0,0,h+2R).對(duì)于外切三角面元的剖分方式，中間每一層各個(gè)節(jié)點(diǎn)的坐標(biāo)可以表示為
zij=h+(i-1)Δd
(7)
對(duì)于內(nèi)接三角面元的剖分方式，中間每一層各個(gè)節(jié)點(diǎn)的坐標(biāo)可以表示為
zij=h+(i-1)Δd
(8)
其中,i=2，…,N+1,j=1,…,P+1，Δθ=2π/P，Δd=2R/N.
1.2.2球體模型的正演結(jié)果
通過(guò)球體模型對(duì)三角面元法進(jìn)行驗(yàn)證.僅考慮磁性體的感應(yīng)磁化強(qiáng)度，并且在均勻磁化的情況下，設(shè)地磁場(chǎng)強(qiáng)度為49 600 nT,磁傾角為44.5°，磁偏角為-2.5°，磁化率為2.設(shè)球體半徑為1 m，中心點(diǎn)埋深為2 m，中心位于坐標(biāo)原點(diǎn)處.圖3所示為N=24、P=12情況下，球體模型產(chǎn)生的ΔT磁異常剖面曲線.從圖3可見(jiàn)，使用外切三角面元計(jì)算得到的球體磁異常曲線比使用內(nèi)接三角面元計(jì)算得到的球體磁異常曲線更接近理論磁異常曲線，與理論曲線重合.外切和內(nèi)接三角面元正演得到球體磁異常的均方根誤差(RMSE)分別為0.22和0.83，這表明不同節(jié)點(diǎn)位置構(gòu)建的外切或內(nèi)接多面體模型，會(huì)對(duì)目標(biāo)磁性體磁異常正演產(chǎn)生一定影響.同樣網(wǎng)格剖分下，外切和內(nèi)接三角面元對(duì)物體表面的逼近速度不同，對(duì)球體模型而言，外切三角面元的逼近速度更快.同時(shí)，進(jìn)一步細(xì)化剖分網(wǎng)格，可以進(jìn)一步提高正演精度.
圖3 不同方法球體磁異常剖面曲線對(duì)比
2 不同姿態(tài)UXO模型的正演方法
UXO大體上可以分為4類：集束彈藥、普通陸戰(zhàn)武器彈藥、航空炸彈、地雷[10].UXO外殼多為鐵磁性物質(zhì)，磁性參數(shù)與周圍環(huán)境的差別較大，在一定條件下能夠被磁法勘探儀器所識(shí)別.然而,UXO外形差別較大.航空炸彈一般體形較大，長(zhǎng)度較長(zhǎng)，外形與圓柱體比較接近，因此可使用均勻有限長(zhǎng)的三度體圓柱體模擬航空炸彈；其他幾類彈藥體形較小，體長(zhǎng)較短，長(zhǎng)度和埋深的比值很小，理論和實(shí)踐已經(jīng)證明這幾類彈藥產(chǎn)生的磁場(chǎng)近似于球體所產(chǎn)生的磁場(chǎng).常見(jiàn)的航空炸彈的外形往往不是標(biāo)準(zhǔn)的有限長(zhǎng)圓柱體，而是變徑有限長(zhǎng)圓柱體，即沿長(zhǎng)軸方向半徑不同，用單一有限長(zhǎng)圓柱體的模型無(wú)法很好地近似，但可以用不同半徑的圓柱體元進(jìn)行組合.
2.1 有限長(zhǎng)圓柱體模型的剖分方法
與球體的剖分方式類似，沿圓柱體長(zhǎng)軸方向切割成N個(gè)等厚薄片，在每個(gè)薄片與柱體相交的圓周上等弧長(zhǎng)地設(shè)立P個(gè)節(jié)點(diǎn)，在每層薄片上用P個(gè)三角形來(lái)擬合薄片表面，如圖4所示.對(duì)于圓柱體的頂/底面，頂/底面圓心與第一層各節(jié)點(diǎn)的連線共形成P個(gè)等腰三角形，而對(duì)于沿圓柱體長(zhǎng)軸切割出的N個(gè)薄片，則能組成2P(N-1)個(gè)三角面元.同樣，在對(duì)每層薄片的表面進(jìn)行擬合時(shí)，也可根據(jù)在每個(gè)薄片與柱體相交的圓周上設(shè)立P個(gè)節(jié)點(diǎn)組成的多邊形是圓外切多邊形還是內(nèi)接多邊形，將三角面元分為外切面元和內(nèi)接面元.
除此之外，如果圓柱體各個(gè)薄片上的對(duì)應(yīng)節(jié)點(diǎn)Dij和Di+1,j(i=1,2,…,N,j=1,2,…,P)位于同一條軸線上，那么按照斜線DijDi+1,j+1分割出的三角形則為直角三角形(見(jiàn)圖4a).如果圓柱體各個(gè)薄片按照奇偶層交錯(cuò)選取節(jié)點(diǎn)，則可形成等腰三角形(見(jiàn)圖4b).
a 直角三角面元
b 等腰三角面元
2.2 全空間任意姿態(tài)有限長(zhǎng)圓柱體磁異常計(jì)算方法
利用三角面元法計(jì)算有限長(zhǎng)圓柱體磁場(chǎng)的優(yōu)勢(shì)在于，不通過(guò)旋轉(zhuǎn)坐標(biāo)系，即不改變觀測(cè)點(diǎn)坐標(biāo)，即可快速計(jì)算出有限長(zhǎng)圓柱體在空間任意姿態(tài)(任意水平和垂直位置、任意走向、任意傾角)下的磁場(chǎng)，減少了計(jì)算量，并能保證快速穩(wěn)定計(jì)算.
假設(shè)起始位置處走向0°為x軸正方向(x軸正方向?yàn)榇疟保蛭鳛樨?fù)，向東為正)，傾角從水平投影面起算為0°(向下為正)，中心點(diǎn)空間位置為(0,0,0)，在目標(biāo)位置處設(shè)走向角度為β，傾角為α，中心點(diǎn)空間位置為(x1,y1,z1).對(duì)圓柱體表面剖分出的任意一個(gè)節(jié)點(diǎn)Dij的三維坐標(biāo)(xij,yij,zij)進(jìn)行投影，先投影到傾角方向，再投影到走向方向，最后平移到目標(biāo)位置，投影運(yùn)算實(shí)際上是投影矩陣與坐標(biāo)相乘.走向方向和傾角方向的投影矩陣分別為Tβ和Tα，如下所示:
(9)
最終可以得到目標(biāo)位置處節(jié)點(diǎn)Dij的坐標(biāo)為
(10)
計(jì)算出每個(gè)節(jié)點(diǎn)的坐標(biāo)后，按照一定順序輸出三角形的3個(gè)角點(diǎn)坐標(biāo)，各個(gè)角點(diǎn)的輸出順序應(yīng)統(tǒng)一為順時(shí)針?lè)较颍?jì)算三角面元的平面法向量時(shí)遵循右手螺旋法則.對(duì)每個(gè)三角面元產(chǎn)生的磁異常進(jìn)行計(jì)算，將所有三角面元產(chǎn)生的磁異常在每一個(gè)觀測(cè)點(diǎn)處累加，即可得到該觀測(cè)點(diǎn)處目標(biāo)體產(chǎn)生的磁異常.
2.3 有限長(zhǎng)圓柱體模型的正演結(jié)果
設(shè)磁場(chǎng)參數(shù)同1.2.2節(jié)，僅考慮磁性體的感應(yīng)磁化強(qiáng)度且均勻磁化，模型體各項(xiàng)物理參數(shù)相同，即磁化率均為2,半徑均為0.1 m，長(zhǎng)度均為2 m.圖5是N=20、P=36下，采用等腰內(nèi)接三角面元的計(jì)算結(jié)果.圖5a和5b是中心點(diǎn)平面坐標(biāo)位于(1,1)、中心點(diǎn)埋深1 m的有限長(zhǎng)圓柱體在不同走向β和傾角α下的形態(tài),圖5c和5d是2種形態(tài)下對(duì)應(yīng)的磁異常.
a 圓柱體1(β=45°,α=0°)
b 圓柱體2(β=45°,α=30°)
c 圓柱體1磁場(chǎng)
d 圓柱體2磁場(chǎng)
圖5不同走向和傾角下有限長(zhǎng)圓柱體的形態(tài)及其產(chǎn)生的磁場(chǎng)
Fig.5Shapeoffinite-lengthcylinderindifferentstrikesanddipanglesandthecorrespondingmagneticfield
從圖5可以看出，磁異常平面圖中磁異常的形態(tài)與位置以及大小反映了模型體的空間位置，比較容易標(biāo)定和識(shí)別，進(jìn)而可以大致確定出UXO的走向和傾向.這一點(diǎn)對(duì)于在實(shí)際探測(cè)中根據(jù)探測(cè)到的磁異常形態(tài)對(duì)UXO進(jìn)行空間定位有重要的理論意義.
對(duì)于中心點(diǎn)坐標(biāo)(0,0,2) m的垂直有限長(zhǎng)圓柱體，使用3種剖分方法得到的三角面元分別計(jì)算出磁異常剖面曲線，對(duì)比目前研究中使用偶極線磁場(chǎng)模擬得到的圓柱體磁異常曲線，如圖6所示.其中，對(duì)于三角面元法，N=20，P=36.
圖6 不同方法垂直有限長(zhǎng)圓柱體磁異常剖面曲線對(duì)比
Fig.6Comparisonofsectioncurvesofverticalfinite-lengthcylindricalmagneticanomalyamongdifferentmethods
理論上三角面元的面積和應(yīng)盡可能接近圓柱體的表面積，則計(jì)算出的圓柱體磁異常可以無(wú)限逼近真實(shí)磁異常.若以圓柱體表面積(此例中為13 194.69 cm2)和三角面元面積和之間的差值作為評(píng)判標(biāo)準(zhǔn)，則在相同的層數(shù)和節(jié)點(diǎn)數(shù)目下3種剖分方式面積和如表1所示.
表1圓柱體三角面元面積和與實(shí)際面積差
Tab.1Differencebetweenthecalculatedareausingtriangularsurfaceelementandtheactualareaofthecylinder
三角面元形式三角面元面積和/cm2面積差絕對(duì)值/cm2等腰外切13228.3833.69等腰內(nèi)接13175.6519.04直角內(nèi)接13175.5619.13
從圖6和表1可以看出,在內(nèi)接情況下等腰和直角三角面元對(duì)磁異常正演結(jié)果的影響較小，但相比于等腰內(nèi)接三角面元等腰外切三角面元對(duì)磁異常正演結(jié)果的影響則較大.理論上外切三角面元計(jì)算得到的圓柱體磁異常比真實(shí)值大，而內(nèi)接三角面元計(jì)算出的圓柱體磁異常比真實(shí)值小，但從圖6可見(jiàn)，采用偶極線計(jì)算得到的磁異常誤差比外切三角面元計(jì)算的磁異常更大，且該誤差無(wú)法修正改進(jìn).對(duì)圓柱體模型而言，內(nèi)接三角面元的逼近速度更快，等腰內(nèi)接三角元計(jì)算出的磁異常最接近真實(shí)值.
值得注意的是，現(xiàn)有的UXO磁異常計(jì)算[5-6]認(rèn)為UXO內(nèi)部鐵磁性物質(zhì)是均勻分布的，但實(shí)際未爆物的磁性主要來(lái)源于UXO鐵磁性外殼，外殼以內(nèi)爆炸物往往是弱磁性或無(wú)磁性.因此，要進(jìn)一步提高正演精度，則必須考慮這種內(nèi)部“中空”磁性結(jié)構(gòu),只有基于表面積分法的三角面元法能夠有效解決該問(wèn)題.將未爆物產(chǎn)生的磁異常與管狀磁性物產(chǎn)生的磁異常近似，與相同空間位置上不同半徑有限長(zhǎng)圓柱體在觀測(cè)點(diǎn)處產(chǎn)生的磁場(chǎng)相減，來(lái)近似計(jì)算均勻磁化中空UXO產(chǎn)生的磁異常.
3 基于三次樣條的UXO磁異?？焖僬莘椒?/h2>
雖然使用三角面元法計(jì)算均勻磁化磁性體磁異常精度高，但是計(jì)算不同尺寸的目標(biāo)體時(shí)需要對(duì)目標(biāo)體表面重復(fù)進(jìn)行網(wǎng)格剖分.隨著三角網(wǎng)格的細(xì)化，計(jì)算量會(huì)提高，導(dǎo)致計(jì)算效率下降.因此，引入三次樣條插值方法，對(duì)不同長(zhǎng)度的模型，不需要重新剖分網(wǎng)格，而是利用一個(gè)已剖分成三角網(wǎng)格的體積元在空間上插值構(gòu)建磁性體，進(jìn)而實(shí)現(xiàn)UXO的磁異?？焖僬?UXO模型近似為有限長(zhǎng)圓柱體或變徑有限長(zhǎng)柱體.
正演思想可以分為以下幾個(gè)步驟：
步驟1將UXO模型磁性體沿長(zhǎng)軸方向均勻切割為幾個(gè)塊體，每個(gè)塊體形態(tài)和尺寸相同或相似(不同塊體可以通過(guò)通過(guò)同一塊體半徑變化獲得).對(duì)其中一塊進(jìn)行上述三角面元剖分，但對(duì)切割面不進(jìn)行剖分，因?yàn)榍懈蠲娌粚儆谀繕?biāo)體整體的表面.對(duì)邊界塊體上的端面要單獨(dú)進(jìn)行剖分，并不計(jì)入以下步驟的計(jì)算.
步驟2對(duì)每一個(gè)觀測(cè)點(diǎn)計(jì)算出小塊體在UXO腔體內(nèi)任意幾處以及邊界位置上的磁異常.
步驟3選取步驟2中幾處位置上的磁異常值作為插值節(jié)點(diǎn)，以塊體的切割間距進(jìn)行插值.
步驟4在每一個(gè)觀測(cè)點(diǎn)處，累加插值得到各個(gè)磁異常值，并補(bǔ)充邊界塊體上的端面產(chǎn)生的磁異常.
3.1 三次樣條插值方法基本原理
Lagrange插值、Newton插值及Hermite插值等高次插值中會(huì)出現(xiàn)龍格(Runge)現(xiàn)象.分段插值雖然可以避免龍格現(xiàn)象,但是分段線性插值在節(jié)點(diǎn)處不光滑。分段Hermite插值僅僅要求插值函數(shù)在節(jié)點(diǎn)上等于被插函數(shù)的一階導(dǎo)數(shù)，并且導(dǎo)數(shù)值往往并不已知.三次樣條插值先由函數(shù)值確定導(dǎo)數(shù)值，再由分段Hermite插值解決問(wèn)題，因而避免了龍格現(xiàn)象，并且插值曲線光滑.
三次樣條函數(shù)的定義如下所示[11]：
設(shè)對(duì)任意區(qū)間[a,b]有剖分Δ∶a=l0 (1)S(l)∈C2[a,b]，即具有連續(xù)的一階、二階導(dǎo)數(shù)，Sk,v(lk)=Sk-1,v(lk),v=0,1,2，k=0,1,…,n-1. (2)S(l)在每一個(gè)小區(qū)間[lk,lk+1]上是次數(shù)≤3的多項(xiàng)式，則稱S(l)為關(guān)于剖分Δ的一個(gè)三次樣條函數(shù). (3) 設(shè)給定函數(shù)表(lq,f(lq)),q=0,1,…,n，若三次樣條函數(shù)S(l)還滿足以下插值條件： S(lq)=f(lq),q=0,1,…,n 則稱S(l)為f(l)關(guān)于剖分Δ的三次樣條插值函數(shù).求取三次樣條插值函數(shù)表達(dá)式，首先以分段三次Hermite插值為基礎(chǔ)，滿足上述條件的情況下，從下述3種邊界條件中的一種推導(dǎo)而來(lái): (1) 第一種邊界條件.S′(l0)、S′(ln)已知，即S′(l0)=f′(l0)，S′(ln)=f′(ln). (2) 第二種邊界條件.S″(l0)、S″(ln)已知，即S″(l0)=f″(l0)，S″(ln)=f″(ln)，若S″(l0)=S″(ln)=0，則稱為自然邊界條件. (4) 非扭結(jié)邊界.使兩端點(diǎn)處多項(xiàng)式三次項(xiàng)系數(shù)與這兩端點(diǎn)的鄰近點(diǎn)的三次項(xiàng)系數(shù)相等，此時(shí)不用給出邊界值. 對(duì)于半徑R=0.1 m、長(zhǎng)度L=2 m、埋深h=1 m、中心點(diǎn)平面坐標(biāo)為(0,0) m的水平有限長(zhǎng)圓柱體，將其沿走向均勻切割成N=40塊，每一塊都是長(zhǎng)L/N的小圓柱體，將每一塊的中心點(diǎn)橫坐標(biāo)位置作為插值點(diǎn)，在空間上對(duì)應(yīng)有N個(gè)插值點(diǎn)c1,…,cN.對(duì)每一塊體進(jìn)行截面節(jié)點(diǎn)P=12的等腰內(nèi)接三角面元剖分，在圓柱體上的對(duì)稱位置任意選擇n=5個(gè)型值點(diǎn)E1,…,E5，對(duì)每個(gè)測(cè)點(diǎn)用三角面元法計(jì)算出單個(gè)塊體在5個(gè)型值點(diǎn)處的磁場(chǎng)值f(Ek)(k=1,…,n)作為插值點(diǎn). 在每一個(gè)觀測(cè)點(diǎn)上，使用自然邊界條件下的三次樣條函數(shù)，插值出塊體在N個(gè)插值點(diǎn)處的磁場(chǎng)值S(ci)(i=1,…,N)，對(duì)得到的結(jié)果進(jìn)行疊加，即得到有限長(zhǎng)圓柱體的磁場(chǎng).圖7a展示了使用插值和不使用插值計(jì)算出的水平圓柱體磁異常剖面曲線對(duì)比，圖7b展示了兩者的絕對(duì)誤差和相對(duì)誤差分布.絕對(duì)誤差分布與磁異常曲線相似，最大值位于端點(diǎn)附近為11 nT,相對(duì)誤差平均值約2.9%.圖7c展示了地面中心觀測(cè)點(diǎn)處，根據(jù)5個(gè)型值點(diǎn)-0.9、-0.5、0、0.5、0.9 m插值得到的磁場(chǎng)值，可見(jiàn)插值曲線非常光滑.將40個(gè)插值點(diǎn)的磁場(chǎng)值累加，即得到中心觀測(cè)點(diǎn)處水平圓柱體整體的磁場(chǎng)值. a 三次樣條插值剖面曲線 b 誤差曲線 c 中心觀測(cè)點(diǎn)處插值曲線 圖8 三角面元構(gòu)造的變徑紡錘形UXO模型 Fig.8FusiformUXOmodelwithvariablediametersmadebytriangularsurfaceelement 對(duì)于變徑有限長(zhǎng)圓柱體，逐段計(jì)算每個(gè)圓柱體元產(chǎn)生的磁異常后疊加，計(jì)算量較大，可以用上述原理進(jìn)行快速正演.插值計(jì)算過(guò)程同上述，對(duì)于在變徑有限長(zhǎng)柱體上的對(duì)稱位置任意選擇n=5個(gè)型值點(diǎn)(-0.9、-0.5、0、0.5、0.9 m)，對(duì)每個(gè)測(cè)點(diǎn)計(jì)算出單個(gè)塊體在n個(gè)型值點(diǎn)處的磁場(chǎng)值作為插值點(diǎn).在每一個(gè)觀測(cè)點(diǎn)上，使用自然邊界條件下的三次樣條函數(shù)，插值出塊體在N個(gè)插值點(diǎn)處的磁場(chǎng)值，對(duì)得到的結(jié)果進(jìn)行疊加，即得到變徑有限長(zhǎng)圓柱體的磁場(chǎng)值.圖9a展示了使用插值和不使用插值計(jì)算出的變徑有限長(zhǎng)圓柱體模型磁異常曲線對(duì)比.圖9b展示了兩者的絕對(duì)誤差和相對(duì)誤差分布，絕對(duì)誤差最大值位于端點(diǎn)附近，為12 nT,平均值為0.7 nT,相對(duì)誤差平均值為0.35%,最大值為8%. 本文針對(duì)工程探測(cè)中遇到的地下多姿態(tài)和形態(tài)各異的均勻磁化UXO的磁正演問(wèn)題，基于面磁荷積分原理，采用三角面元法的和三次樣條插值方法，實(shí)現(xiàn)了高精度快速正演模擬，并進(jìn)行了相關(guān)問(wèn)題的分析. 本文提出的正演計(jì)算方法，有較好的計(jì)算速度和正演精度，可以模擬不規(guī)則外形以及任意空間位置下空心UXO的磁異常.該方法可進(jìn)一步與反演算法結(jié)合，以提高UXO定位的速度和準(zhǔn)確性，也可應(yīng)用于地下復(fù)雜金屬管道的磁正演. a 三次樣條插值剖面曲線 b 誤差曲線 Fig.9Resultsofvariablediametercylindricalmagneticanomalycalculatedwithsplineinterpolation 參考文獻(xiàn): [1] ZALEVSKY Z,BREGMAN Y,SALOMONSKI N,etal.Resolution enhanced magnetic sensing system for wide coverage real time UXO detection[J].Journal of Applied Geophysics,2012,84(4):70. [2] CHURCHILL K M,LINK C,YOUMANS C C.A comparison of the finite-element method and analytical method for modeling unexploded ordnance using magnetometer[J].IEEE Transactions on Geoscience and Remote Sensing,2012,50(7):2720. [3] SINEX D B.Advancing the state of the art of UXO discrimination for total-field magnetic data[D].Golden:Colorado School of Mines,2006. [4] SANCHEZ V,LI Y G,NABIGHIAN M N,etal.Numerical modeling of higher order magnetic moments in UXO discrimination[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(9):2568. [5] 葛健，陸承達(dá)，董浩斌，等. 基于Overhauser傳感器的近地表UXO磁梯度法探測(cè)技術(shù)[J]. 儀器儀表學(xué)報(bào)，2015，36(5)：961. GE Jian,LU Chengda,DONG Haobin,etal.The detection technology of near-surface UXO based on magnetic gradient method and Overhauser sensor[J].Chinese Journal of Scientific Instrument,2015，36(5)：961. [6] 郭志馗，陳超，陶春輝，等. 有限長(zhǎng)圓柱體磁異常場(chǎng)全空間正演方法[J]. 地球物理學(xué)報(bào), 2017,60(4):1557. GUO Zhikui,CHEN Chao,TAO Chunhui,etal.Forward modeling of magnetic anomaly of finite-length cylinder in 3D space[J].Chinese Journal of Geophysics,2017,60(4):1557. [7] 朱慧慧，劉得軍，馮碩,等. 基于磁偶極子構(gòu)造法幾何建模的鐵管磁異常正演[J]. 計(jì)算機(jī)測(cè)量與控制,2017,25(3):201. ZHU Huihui,LIU Dejun,FENG Shuo,etal.Magnetic dipole reconstruction geometry modeling for underground ferromagnetic pipe magnetic abnormal detection[J].Computer Measurement and Control,2017,25(3):201. [8] 吳文鸝, 管志寧, 高艷芳,等. 重磁異常數(shù)據(jù)三維人機(jī)聯(lián)作模擬[J]. 物探化探計(jì)算技術(shù), 2005, 27(3):227. WU Wenli,GUAN Zhining,GAO Yanfang,etal.Interactive man/computer modeling 3D body of gravity and magnetic anomaly data[J].Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(3):227. [9] 管志寧.地磁場(chǎng)與磁力勘探[M].北京：地質(zhì)出版社，2005. GUAN Zhining.Geomagnetic field and magnetic exploration[M].Beijing:Geological Publishing House,2005. [10] 李小康，楊磊，劉磊.未爆彈藥問(wèn)題及地球物理解決方案綜述[J].中國(guó)礦業(yè)，2010，19:187. LI Xiaokang,YANG Lei,LIU Lei.Overview on unexploded ordnance problem and a solution: the geophysical scheme[J].China Mining Magzine,2010，19:187. [11] 關(guān)冶，陸金甫.數(shù)值方法[M].北京：清華大學(xué)出版社,2011. GUAN Ye,LU Jinfu.Numerical method[M].Beijing:Tsinghua University House,2011.3.2 不同模型的正演結(jié)果
4 結(jié)語(yǔ)
猜你喜歡
遼寧師范大學(xué)學(xué)報(bào)(自然科學(xué)版)(2022年4期)2022-12-27 12:27:28
中等數(shù)學(xué)(2022年1期)2022-06-05 07:50:08
南京大學(xué)學(xué)報(bào)(數(shù)學(xué)半年刊)(2020年1期)2020-03-19 02:24:48
赤峰學(xué)院學(xué)報(bào)·自然科學(xué)版(2018年9期)2018-10-18 07:15:52
中等數(shù)學(xué)(2018年4期)2018-08-01 06:36:34
中學(xué)數(shù)學(xué)研究(廣東)(2018年23期)2018-03-05 07:54:32
初中生世界·九年級(jí)(2017年9期)2017-10-13 10:27:04
裝備學(xué)院學(xué)報(bào)(2016年5期)2016-11-14 06:26:42
池州學(xué)院學(xué)報(bào)(2015年3期)2016-01-05 01:13:07
地震地質(zhì)(2015年3期)2015-12-25 03:29:42
