磁力張量數據的目標體邊界識別方法探討

秦秀穎， 魯光銀， 朱自強， 曹書錦

(中南大學　地球科學與信息物理學院，長沙　410083)

?

磁力張量數據的目標體邊界識別方法探討

秦秀穎， 魯光銀， 朱自強， 曹書錦

(中南大學地球科學與信息物理學院，長沙410083)

摘要：這里系統地介紹了磁張量數據邊界識別數值計算的幾種方法，包括利用磁張量不變量、水平梯度模法(the total horizontal derivatives，THDR)、解析信號法(The analytic signal amplitude，ASM)、Tilt梯度法、Tilt梯度的水平導數法、Theta圖法、Bspan的水平梯度模法等，給出了各種方法的張量組合方式、各自的適用范圍及優缺點，利用結果圖中的極值點、零點或其他特征點來識別邊界。通過對矩形棱柱體的模型驗算，證明深度、磁傾角對邊界識別影響較大，磁偏角幾乎無影響。利用磁梯度不變量、解析信號法、Tilt梯度法得到目標體邊界精確、完整，與實際目標體邊界吻合度高，表明利用這些方法能夠很好地完成磁張量數據的邊界識別。

關鍵詞：磁張量； 張量不變量; 磁源邊界; 邊界提取

0引言

在研究地質目標體的橫向不均勻性特別是地質目標體的邊緣位置時，磁位場有其獨特的優勢。這里的地質目標體邊緣是指斷裂構造線、不同地質體的邊界線，實際上是具有一定磁性差異的地質體的邊界線。磁力梯度測量作為一種新型的地球物理勘探方法，具有快速、經濟、范圍大等優點，是劃分大地構造單元，進行構造分區，確定斷裂構造帶的位置，區分不同巖性與地層的分布，進行物性填圖快速而有效的手段[1]。磁力梯度張量測量系統的測量對象是磁場矢量分量的梯度，磁測結果能夠反映目標體的矢量磁矩信息，可用于描述場源體的磁化方向和幾何形態，提高對磁源體的分辨率。隨著工藝技術的發展，磁張量測量已經成為現實[2]?，F今有關張量數據處理方法的研究受到人們廣泛地關注，但大多是針對重力張量數據，而對于磁張量數據的研究方法較少[3]。在地質體的邊緣附近，磁梯度是磁場中場變化最劇烈的地方，磁異常變化率較大，它反應了地質體的邊界，故所有的邊緣識別方法均利用這一特點進行設計[4]。利用磁位場識別地質體邊緣的方法分為數理統計，數值計算和其他三大類，這里主要介紹數值計算類[5]，它包括利用磁張量不變量、水平梯度模法、解析信號法、Tilt梯度法、Tilt梯度的水平導數法、Theta圖法、BZZ的水平梯度模法等，并設置多個模型來驗證各種方法的應用效果與優勢。

1邊界提取方法

磁法勘探中測量的物理量有地球磁場總場的模量、總磁場模量的空間變化率、地球磁場的三個分量、磁場三分量的空間變化率，即磁梯度張量。

磁梯度張量是磁場矢量的三個分量在相互正交的三個方向的空間變化率，磁場為矢量場，在形式上，梯度張量可以表示為兩個三元素矢量的矩陣乘法運算：

(1)

其中:U為磁標勢；Bx、By和Bz是磁場矢量B在空間三個方向上的分量。

由磁場強度旋度為零可知，磁梯度張量是對稱的，即：

Bxy=Byx,Byz=Bzy,Bxz=Bzx

(2)

而且在無源區，磁勢滿足拉普拉斯方程：

▽·U=0

(3)

則磁梯度張量對角元素之和為零，即：

Bxx+Byy+Bzz=0

(4)

根據式(2)和(4)可知，磁梯度張量的九個元素中，只有五個元素是獨立的，則可將其表示為：

(5)

綜上可得磁場三個矢量分量、總場強度以及磁梯度張量各元素在直角坐標系中的關系如圖1所示。

地質體邊界的圈定一直是磁場數據解釋中的重要內容，可以用很多方法對地質體邊界進行識別[6]，數值計算類邊界識別方法是研究最多、應用最廣的邊界識別方法。邊界識別結果是一種定性或半定量解釋結果，與定量解釋結果有一定區別，識別結果可作為邊界位置定量反演的初值。直接根據原始磁異常不能準確地劃分出地質體的邊界，但對某個張量數據單獨進行處理解釋，忽略有用信息，容易造成信號的不合理解釋[7]。張量分量的不同組合所產生的不變量形式簡單、容易處理解釋，并且與磁異常體的大小、形狀相關，集中把多個分量的信息增強，能夠更好地處理和解釋數據[8]。目前常用的磁張量數據邊界識別數值計算的方法，包括利用水平梯度模法、解析信號法、Tilt梯度法、Tilt梯度的水平導數法、Theta圖法、BZZ的水平梯度模法等，下面研究它們的張量組合方式、各自的適用范圍及優缺點。

由實對稱矩陣的性質可知，磁梯度張量矩陣具有三個旋轉不變量，這三個不變量不受坐標系統旋轉變換影響，分別表示為式(6)[8]：

BxyBzz+Bxz(BxyByz-BxzByy))1/3

(6)

這三個張量組合同時包含了五個獨立分量，上述磁力張量不變量I1本身具有很好的性質，如不變量等值線圖不受地磁場方向影響，能夠很好地繪出場源邊界，而I2能很好地突出淺部異常體的特征。

圖1　磁場矢量分量、總場強度以及磁梯度張量各元素關系圖Fig．1　The relationship diagram of magnetic field vector components, the total field strength and the magnetic gradient tensor elements(a)磁場分量、總場強度；(b)磁梯度張量各元素間關系

水平梯度模也是一種常用的邊界識別方法，余欽范[4]指出了可以利用水平梯度極大值法識別磁性體邊界，利用重力位與磁位的泊松關系，把磁異常換算為偽重力異常(磁源重力異常)，然后計算偽重力異常水平梯度的模(強度)，經過處理后的水平梯度模便可用來確定物性邊界的位置。在水平梯度模極大值位置處繪出標志以表示邊界位置，其表達式為式(7)[9]。

(7)

最近幾年，國內、外文獻有不少關于利用磁異?？偺荻饶＿M行磁測資料解釋的報導，它能迅速推斷磁源的水平邊界位置。Roest等[11]人提出了利用三維解析信號法(總梯度模)進行磁異常解釋?？偺荻饶７?，又稱為解析信號法，其表達式為式(8)[10]：

(8)

解析信號振幅的極大值出現在磁性差異處，利用這些極大值可以確定目標體邊界的位置。總梯度模法在磁異常方面研究較多，它受磁化方向影響小，能夠準確地識別有一定走向的垂直邊界、無限延深厚板等磁性體的位置。吳文賢[10]用模型試驗表明了該方法在低緯度地區磁異常分析與解釋方面的有效性,并將基于總梯度模的磁異常分析方法應用于老撾爬立山鐵礦區的磁測資料處理，由此排除了由化極不穩定導致的虛假磁異常，由該方法確定的磁異常區域與區內磁鐵礦體分布位置的吻合性較好,圈定了礦體異常的分布范圍，進一步指導礦區鐵礦勘探。

Tilt梯度是用垂直梯度比上水平梯度的絕對值的arctan角度(式(9))[11]：

(9)

Hood[15]探討了應用垂直梯度法和水平梯度法來識別磁性體的垂直邊界；Grauch[16]提出用最大水平梯度法探測傾斜臺階的傾斜邊界；余欽范[4]利用水平梯度極大值法識別密度和磁性體邊界，假定磁源邊界為單個垂直邊界，磁異常水平梯度模在邊界上恰好取得極大值，由此便可確定磁源邊界的位置。

Theta圖法利用解析信號振幅對總水平導數進行歸一化，其θ角的計算公式如式(10)所示[12]：

(10)

式中cosθ的取值在0～1之間，在場源邊界處的取值趨近于1。Wijins[14]完成了二度體磁力異常和三度體磁力異常的模型試算以及實際資料處理，指出此法不受磁異常分量和磁化方向的影響，比解析信號振幅的分辨能力更強。

HDT法(horizontaldirectivetendencyordifferentialcurvaturemagnitude)，也可稱作水平方向性或是微分曲率幅值[13]，其表達式為式(11)。

(11)

2模型試算

為了檢驗上述幾種邊界識別方法的應用效果，現對矩形棱柱體模型引起的磁張量異常進行處理。

模型一：場源空間大小為1 000m×1 000m×500m，場源空間存在一個400m×400m×200m的矩形棱柱體，棱柱體頂部埋深為50m，底部埋深為450m，長方體的水平中心位置為(500m，500m)，磁化強度為1nT/m，x、y方向點距為50 m，測線長度為1 000 m。

當磁傾角為90°，磁偏角為0°時，對于上面提到的6種磁梯度測量邊界提取方式分別進行了計算，所得等值線圖如圖2～圖7所示。

圖2為磁張量不變量I2的結果等值線圖，圖2中白色虛線表示實際目標體的邊界?？梢钥闯觯纱艔埩坎蛔兞縄2等值線圖，可以很好地勾繪出實際目標體邊界。

圖2　磁張量不變量I2等值線圖及邊界提取結果Fig．2　Magnetic tensor invariant I2contour 　　　map and boundary recognition results

圖3　水平梯度模等值線圖及邊界提取結果Fig．3　Horizontal gradient mode contour map 　　　and boundary recognition results

圖3為水平梯度法的等值線圖及邊界提取結果。其中白色虛線代表實際目標體邊界。從圖3中可以看出，在水平梯度模的極大值處描繪的標志，與實際目標體邊界吻合較好。

圖4　解析信號法的等值線圖及邊界提取結果Fig．4　Analytical method of signal contour 　　　map and boundary recognition results

圖4是利用解析信號法得到的邊界提取結果等值線示意圖。圖4中白色虛線表示實際目標體的邊界。對于有一定水平尺度的磁性體，磁異?？偺荻饶O大值與磁性體淺部邊界有較好的對應關系，利用ASM可以較準確地圈定磁性地質體的邊界位置。其表達式相對簡單，具有較強的分辨迭加異常的特性。從圖4中可以看出，利用總梯度模能夠很好地確定磁源邊界。

圖5　Tilt梯度等值線圖及邊界提取結果Fig．5　Tilt gradient contour map and 　　　boundary recognition results

圖5為Tilt梯度等值線圖，圖5中紅色虛線表示實際目標體的邊界。Miller[17]等指出，Tilt梯度在場源上為正值，場源邊界處取為零值，在場源外部為負值。Tilt梯度的零等值線位置對應物體源的邊界，從圖5中可以看出，Tilt梯度的零等值線與目標體的理論邊界剛好重合。

圖6　Theta圖的等值線圖及邊界提取結果Fig．6　Theta map contour map and 　　　boundary recognition results

圖6顯示了利用Theta圖計算的等值線圖及邊界提取結果。其中白色虛線代表實際目標體邊界。在cosθ取值趨近于“1”時繪制目標體邊界。從圖6中可以看出，所得到結果與實際邊界吻合度較高。

圖7　HDT法的等值線圖及邊界提取結果Fig．7　HDT method contour map and 　　　boundary recognition results

圖7顯示了HDT法的等值線圖及邊界提取結果。其中白色虛線代表實際目標體邊界。從圖7中可以看出，由磁異常的微分曲率幅值的最大值可以確定邊界的四個角點，從而確定目標體邊界。

圖8　模型二邊界提取結果等值線圖Fig．8　Edge recognition results based on the model 2(a)磁張量不變量I2; (b)水平梯度模;(c)解析信號法;(d)Tilt梯度;(e)Theta圖;(f)HDT法

為驗證深度對邊界識別效果的影響，現設置模型二：場源空間大小 1 000m×1 000m×500m，場源空間存在一個400m×400m×200m的矩形棱柱體，改變棱柱體頂部埋深為100m，底部埋深為500m，長方體的水平中心位置為(500m，500m)，磁化強度為1nT/m，x、y方向點距為50 m，測線長度為1 000 m。計算結果如圖8所示。

從圖 8中可以看出識別效果受異常體埋深影響，圖9為選用解析信號法提取邊界結果隨深度變化的誤差分布曲線。中心位置深度變化在(50 m，240 m),從圖中可看出計算誤差在(0，0.2)區間范圍內，并且隨深度的增加先降后升，即識別能力隨著地下目標體埋深的增大先升后降，中心位置深度在(130 m，190 m)區間內可以得到較好的邊界識別結果。

圖9　解析信號法提取邊界結果隨深度變化誤差分布曲線Fig．9　The error distribution curves of the results 　　　of analytic signal method vary with depth

為驗證不同磁傾角、磁偏角磁異常體邊界識別效果，現設置模型三：模型三(A)當磁偏角為45°，磁傾角為90°；模型三(B)當磁偏角為0°，磁傾角為75°。對上述邊界識別方法進行計算得到等值線圖如圖10所示。

由圖10、圖11中可以看出，改變磁偏角大小，這幾種方法的計算結果等值線圖幾乎沒有改變，與實際邊界吻合度仍舊較高，說明這幾種邊界識別方法不受磁偏角的影響；當磁傾角改變時，幾種邊界提取方法的計算結果與實際目標體的邊界有一定的偏差，圖12為解析信號法提取邊界結果隨磁傾角變化誤差分布曲線，磁傾角變化范圍(0°，90°)，邊界識別結果受磁傾角變化的影響較為復雜，誤差曲線前半部分呈下降趨勢，傾角越大，誤差越小，在磁傾角為60°時可取得較為理想的邊界識別效果；后半部分隨著磁傾角的增大誤差先增大后減小。

圖10　模型三(A)邊界提取結果等值線圖Fig．10　Edge recognition results based on the model 3(A)(a)磁張量不變量I2；(b)水平梯度模；(c)解析信號法；(d)Tilt梯度；(e)Theta圖；(f)HDT法

圖11　模型三(B)邊界識別結果等值線圖Fig．11　Edge recognition results based on the model 3(B)(a)磁張量不變量I2；(b)水平梯度模；(c)解析信號法；(d)Tilt梯度；(e)Theta圖；(f)HDT法

圖12　解析信號法提取邊界結果隨磁傾角　　　　　　　變化誤差分布曲線Fig．12　The error distribution curves of the results 　　　of analytic signal method vary with 　　　magnetic inclination

為了驗證上述幾種邊界識別方法在多個地質體模型中的應用效果，現設置模型四(A)：場源空間大小為 1 000 m×1 000 m×500 m，如圖所示，場源空間存在兩個矩形棱柱體，磁化強度為1 nT/m，左側棱柱體頂部埋深為50 m，頂部中心位置為(300 m，500 m)，頂部尺寸為400 m×200 m,底部埋深為450 m。右側棱柱體頂部埋深為50 m，頂部中心位置為(700 m,500 m)，頂部尺寸為400 m×200 m,底部埋深為450 m。x、y方向點距為50 m，測線長度為1 000 m。磁偏角為0°，磁傾角為90°。通過計算所得模型邊界提取結果等值線圖如圖13所示。

模型四(B)：場源空間大小為1 000 m×1 000 m×500 m，場源空間存在三個矩形棱柱體，磁化強度為1 nT/m，左側棱柱體頂部埋深為50 m，頂部中心位置為(200 m，700 m)，頂部尺寸為300 m×150 m,底部埋深為350 m。中間棱柱體頂部埋深為50 m，頂部中心位置為(500 m，500 m)，頂部尺寸為300 m×150 m,底部埋深為350 m。右側棱柱體頂部埋深為50 m，頂部中心位置為(800 m,300 m)，頂部尺寸為300 m×150 m,底部埋深為350 m。x、y方向點距為50 m，測線長度為1 000 m。磁偏角為0°，磁傾角為90°。通過計算所得模型邊界提取結果等值線圖如圖14所示。

從圖13、圖14中可以看出，當出現多個地質體模塊時，上述幾種數值邊界識別方法均可標記出目標體邊界，與磁異常體排列順序無關。其中直接利用磁張量不變量、解析信號、Tilt梯度方法所得結果，與實際邊界吻合度更高，能夠更好地反映目標體的位置和大小，相比于其他幾種方法更為精準。

為了檢驗上述幾種邊界識別方法的抗干擾能力,在模型四(B)的基礎上加入12%的隨機噪聲，各邊界識別方法計算結果如圖15所示。

圖15為對上述各個張量組合分別加入12%的隨機噪聲的邊界提取結果，從圖15中可以看出在隨機噪聲使得各等值線圖不平滑，這些不平滑影響了異常體邊界識別結果。但仍舊可以得到與實際目標磁性體吻合的邊界，可看出幾種邊界識別方法有較好的抗干擾能力。

圖13　模型四(A)邊界識別結果等值線圖Fig．13　Edge recognition results based on the model 4(A)(a)磁張量不變量I2；(b)水平梯度模；(c)解析信號法；(d)Tilt梯度；(e)Theta圖；(f)HDT法

圖14　模型四(B)邊界識別結果等值線圖Fig．14　Edge recognition results based on the model 4(B)(a)磁張量不變量I2；(b)水平梯度模；(c)解析信號法； (d)Tilt梯度；(e)Theta圖；(f)HDT法

圖15　模型四(C)邊界識別結果等值線圖Fig．15　Edge recognition results based on the model 4(C)(a)磁張量不變量I2；(b)水平梯度模；(c)解析信號法； (d)Tilt梯度；(e)Theta圖；(f)HDT法

3結論

作者系統地闡述了磁梯度張量數據邊界提取的各種方法，并均進行計算得出結果，針對解析信號法計算出隨深度及磁傾角變化的誤差分布曲線。識別方法隨著地質體埋深的增大而識別能力先上升后下降；在埋深相同的情況下，識別能力不受磁偏角的影響；受磁傾角影響情況較為復雜。通過對矩形棱柱體的結果檢驗，說明利用磁梯度不變量、解析信號法、Tilt梯度相對于其他邊界提取算法所得的邊界，能更好地反映磁源的邊界位置，特別是對于埋深較淺的垂直棱柱體，其邊界與實際模型吻合得非常好。在多個矩形棱柱體模型驗算中，幾種方法所得邊界簡單，清晰，易于識別，精確度依舊很高，抗干擾能力較強，證明方法可行。磁梯度張量數據，對于直接勾劃異常體的邊界、圈定地下異常體具有更高的精度，將其應用于實際磁梯度張量測量時，有一定的指導意義。

參考文獻：

[1]馬國慶, 李麗麗, 杜曉娟. 磁張量數據的邊界識別和解釋方法[J]. 石油地球物理勘探, 2012, 47(5): 815-821.

MA G Q, LI L L,DU X J.Boundary detection and interpretation by magnetic gradient tensor data[J]. OGP, 2012, 47(5): 815-821.(In Chinese)

[2]吳招才, 劉天佑. 磁力梯度張量測量及應用[J]. 地質科技情報, 2008, 27(3): 107-110.

WU Z C, LIU T Y. Magnetic gradient tensor: its properties and uses in geophysics[J]. Geological Science and Technology Information, 2008, 27(3): 107-110.(In Chinese)

[3]BEIKI M, PEDERSEN L B, NAZI H. Interpretation of aeromagnetic data using eigenvector analysis of pseudo gravity gradient tensor[J]. Geophysics, 2011, 76(3): L1-L10.

[4]余欽范, 樓海. 水平梯度法提取重磁源邊界位置[J]. 物探化探計算技術, 1994, 16(4): 363-367.

YU Q F, LOU H. Locating the boundaries of magnetic or gravity sources using horizon gradient anomalies[J]. Computing Techniques for Geophysical and Geochemical Explosion,1994, 16(4): 363-367.(In Chinese)

[5]王萬銀, 邱之云, 楊永. 位場邊緣識別方法研究進展[J].地球物理學進展, 2010, 25(1): 196-210.

WANG W Y, QIU Z Y, YANG Y. Some advances in the edge recognition of the potential field[J].Progress in Geophysics,2010, 25(1): 196-210.(In Chinese)

[6]BEIKI M, PEDERSEN L B. Eigenvector analysis of gravity gradient tensor to locate geologic bodies[J]. Geophysics, 2010, 75(6): 137-149.

[7]朱自強, 曾思紅, 魯光銀. 重力張量數據的目標邊緣檢測方法探討[J].石油地球物理勘探, 2011, 46(3): 482-488.

ZHU Z Q, ZENG S H, LU G Y. Discussion on the edge detection technique for gravity gradient tensor data[J]. OGP,2011, 46(3): 482-488.(In Chinese)

[8]AFIF H. SAAD. Understanding gravity gradient statutorily[J]. The Leading Edge,2006, 25(8): 942-947.

[9]PEDERSEN L B, RASMUSSEN T M The gradient tensor of Potential field anomalies: some imp- lications on data collection and data processing of maps [J]. Geophysics, 1990, 55(12): 1558-1566.

[10]吳文賢，范文玉，王永華. 基于總梯度模分析的磁源體邊界信息提取及找礦意義[J].地質與勘探,2013,49(6):1164-1169.

WU W X,FAN W Y,WANG Y H.Based on the analysis of total gradient mode of magnetic source body boundary information extraction and prospecting significance[J].Geology and Prospecting,2013,49(6):1164-1169.(In Chinese)

[11]WR ROEST J, VERHOEF, M PILKINGTON, et al. Magnetic interpretation using the 3-D analytic signal [J]. Geophysics, 1992, 57(l): 116-125.

[12]GRAUCH, V.J.S.,CORDELL, LINDRITH.Mapping basement magnetization zones from aeromagnetic data in the San Juan basin,New Mexico[A].//Hinzc W J.The utility of regional and magnetic anomaly[C].Society of Exploration Geophysicists,1985:181- 197.

[13]BRUNO VERDUZEO J, DEREK FAIRHEAD, CHRIS M. GREEN, et al. The meter reader-new in- sights into magnetic derivatives for structural mapping [J]. The Leading Edge, 2004,23(2):116-119.

[14]CHRIS WIJNS, CARLOS PEREZ, PETER KOWALCZYK. Theta map: Edge detection in magnetic data[J]. Geophysies,2005, 70(4): L39-L43.

[15]HOOD P. Aeromagnetic surveying[M].US:Geophysics. Springer US, 1989.

[16]GRAUCH V J S. Limitations of determining density or magnetic boundaries from the horizontal gradient of gravity or pseudogravity data[J]. Geophysics, 1987, 52(1):118-121.

[17]HUGH G MILLER, VIJAY SINGH. Potential field tilt-a new concept for location of potential fied sources[J].Journal of Applied Geophysics, 1994, 32(2-3):213-217.

Discussion of magnetic tensor data target body boundary identification methods

QIN Xiu-ying, LU Guang-yin, ZHU Zi-qiang, CAO Shu-jin

(School of geosciences and info-physics, Central South University, Changsha410083,China)

Abstract:Several methods of the magnetic tensor data boundary identification of numerical calculation are systematically introduced in this paper, including using the magnetic tensor invariant, horizontal gradient modulus method (the total horizontal derivatives, THDR), analytic signal method, the analytic signal amplitude(ASM), Tilt gradient method, the horizontal derivative of Tilt gradient method, Theta map method and the level of the BZZgradient modulus method. Tensor combination of various methods, their respective applicable scope and the advantages and disadvantages are given. Using the results the extreme value point, zero point, or other feature points in the graph to identify boundary. Through the calculation of rectangular prism model, proving that the depth, dip angle has greater influence on the boundary identification, magnetic declination almost had no effect. The geological edges by using the method of magnetic gradient in-variants, analytic signal and Tilt gradient method are accurate and complete. The result is in accordance with model edge. Showing that is able to complete the boundary identification of the magnetic tensor data.

Key words:magnetic tensor; tensor invariant; magnetic source boundary; boundary extraction

中圖分類號：P 631.2

文獻標志碼：A

DOI：10.3969/j.issn.1001-1749.2016.01.03

文章編號：1001-1749(2016)01-0019-11

作者簡介：秦秀穎(1989-)，女,碩士,主要從事重磁數據處理與解釋方面的研究,E-mail:iwanted@126.com。

基金項目:國家自然科學基金(41174061,41374120)

收稿日期：2015-01-04改回日期：2015-04-01