非線性共軛梯度三維反演在山西河津鐵礦音頻大地電磁數(shù)據(jù)勘察中的應(yīng)用

鄒宗霖，翁愛華，周子琨，連鑫葆，郭俊豪

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130026

0 引言

河津市位于山西省的西南部，在大地構(gòu)造上處于秦嶺構(gòu)造帶北緣，西側(cè)是祁呂賀蘭山字型構(gòu)造前弧，東側(cè)是新華夏系第三隆起帶[1]。

研究區(qū)在河津市義唐村的東北部(圖1)，位于新絳至河津含鎳富鐵礦成礦遠(yuǎn)景區(qū)內(nèi)。地層上部為新蓋層(寒武系、奧陶系、石炭系和二疊系地層)，在古老地層中侵入有基性--酸性巖體，這些侵入活動(dòng)伴隨大量礦液，形成熱液型富鐵礦[2]。山西省地球物理化學(xué)勘查院于2008年在河津市外圍開展了詳細(xì)的勘察研究工作，發(fā)現(xiàn)了董村、劉西咀和陽山莊磁異常，經(jīng)過鉆探驗(yàn)證均為富鐵礦。河津市東北部義唐村附近后期發(fā)現(xiàn)了新的航磁異常，為了確定異常的電性特征及分布規(guī)律，山西省地球物理化學(xué)勘查院在異常區(qū)布置了兩條正交的音頻大地電磁剖面。

(a)研究區(qū)在華北克拉通上的位置；(b)研究區(qū)在新絳至河津熱液型鐵礦成礦遠(yuǎn)景區(qū)中的位置；(c)音頻大地電磁測(cè)點(diǎn)分布。圖1 工區(qū)位置圖Fig.1 Location of study area

音頻大地電磁法(AMT)是利用較高頻率的電磁信號(hào)進(jìn)行地下電性結(jié)構(gòu)探測(cè)的一種大地電磁方法[3]，主要用于較淺的工程勘察和礦產(chǎn)勘察之中[4--6]。過去由于反演理論和技術(shù)的約束，對(duì)于剖面性的大地電磁數(shù)據(jù)，主要進(jìn)行二維(2D)反演。但進(jìn)行二維反演的前提是異常體呈線狀沿走向發(fā)育，并且整條剖面具有相同的構(gòu)造走向[7--8]。然而，該地區(qū)位于斷陷盆地內(nèi)，發(fā)育較多的褶皺和斷裂構(gòu)造，異常體呈現(xiàn)非常明顯的三維特征，因此，二維反演很難滿足這種條件。對(duì)于三維體引起的異常，胡祖志等[9]的研究表明，進(jìn)行剖面觀測(cè)并借助二維大地電磁反演，獲得的剖面下方電阻異常特征卻不能有效反映出地下異常體的真實(shí)分布。

目前大地電磁三維(3D)反演方法和理論已基本成熟[10]，并廣泛應(yīng)用到三維數(shù)據(jù)的反演研究中。林昌洪等[11]討論了三維反演用于單條剖面或稀疏測(cè)線中理論模型數(shù)據(jù)反演中，發(fā)現(xiàn)三維反演能夠較可靠恢復(fù)出剖面域內(nèi)和剖面外附近的三維電性異常體。

因此，為獲取研究區(qū)深部真實(shí)電阻率結(jié)構(gòu)，圈定更為可靠的鐵礦體分布有利區(qū)，筆者選用非線性共軛梯度(NLCG)法，針對(duì)兩條正交的大地電磁剖面主阻抗數(shù)據(jù)進(jìn)行3D反演。通過三維反演，電性結(jié)果不僅確認(rèn)了剖面正交位置下方根據(jù)航磁異常推測(cè)的礦體，還發(fā)現(xiàn)了新的低阻異常，并推測(cè)可能為礦體的反映。

1 音頻大地電磁測(cè)深數(shù)據(jù)

音頻大地電磁測(cè)深法(AMT)觀測(cè)天然電磁場(chǎng)的時(shí)間序列信號(hào)，然后將時(shí)間域數(shù)據(jù)轉(zhuǎn)化為頻率域數(shù)據(jù)，進(jìn)而計(jì)算出每個(gè)頻點(diǎn)的電阻率值和阻抗相位[12]。根據(jù)Maxwell方程可得，在均勻大地和水平層狀大地情況下，波阻抗是電場(chǎng)E和磁場(chǎng)H的水平分量的比值[13--14]。

(1)

(2)

(3)

式中：f為頻率；ρ為電阻率；E為電場(chǎng)強(qiáng)度；H為磁場(chǎng)強(qiáng)度；φE為電場(chǎng)相位；φH為磁場(chǎng)相位。

本次音頻大地電磁工作共布設(shè)2條測(cè)線，呈十字相交，其中2#線33個(gè)測(cè)點(diǎn)，點(diǎn)距100 m，測(cè)線長(zhǎng)3.2 km；3#線36個(gè)測(cè)點(diǎn)，點(diǎn)距80 m，測(cè)線長(zhǎng)3.0 km(204點(diǎn)至236點(diǎn)間隔320 m)，測(cè)點(diǎn)共計(jì)69個(gè)，具體分布見圖1。采集儀器采用加拿大鳳凰公司的V5--2000型大地電磁測(cè)深儀，有效響應(yīng)頻率范圍為10 400～0.1 Hz。數(shù)據(jù)采集采用張量測(cè)量方式，平均采集時(shí)間約30 min。

數(shù)據(jù)處理主要是將采集到的時(shí)間域信號(hào)經(jīng)傅里葉變換為頻率域信號(hào)，并通過頻譜計(jì)算視電阻率及相位[15--17]。對(duì)于在村莊、公路等干擾較大的地方可采用磁場(chǎng)遠(yuǎn)參考技術(shù)[18--19]以降低相關(guān)噪聲的干擾，并基于Nyquist圖對(duì)阻抗進(jìn)行編輯[20]。圖2為本次反演數(shù)據(jù)316號(hào)點(diǎn)原始阻抗Zxx、Zxy、Zyx和Zyy的實(shí)部和虛部響應(yīng)曲線。采集到的數(shù)據(jù)主阻抗質(zhì)量較好，輔阻抗質(zhì)量較差，全部分量反演會(huì)造成擬合效果不好故只選擇主阻抗進(jìn)行反演。經(jīng)過前人研究，在正交測(cè)線的三維反演中，主阻抗對(duì)于三維反演結(jié)果的貢獻(xiàn)是非常大的，可以較好的確定異常體的位置，加入輔阻抗后結(jié)果可能會(huì)有改善但變化不大[21]。所以本次三維反演只針對(duì)主阻抗進(jìn)行反演是可行的。

由于三維反演需要在統(tǒng)一的觀測(cè)坐標(biāo)中進(jìn)行，而剖面測(cè)量時(shí)，是在獨(dú)立的測(cè)線坐標(biāo)系中進(jìn)行大地電磁矢量測(cè)量站點(diǎn)布置。因此需要把圖1中兩條測(cè)線中沿各自剖面方向作為北方向測(cè)量的阻抗旋轉(zhuǎn)到統(tǒng)一的真實(shí)北方向上[22--23]，之后再進(jìn)行三維反演。

2 三維反演

三維反演采用的NLCG法特點(diǎn)是在求取目標(biāo)函數(shù)的極小值時(shí)沿著共軛梯度方向進(jìn)行一維搜索，迭代中只需計(jì)算一階導(dǎo)數(shù)[24--25]。該算法并不是在全局求解最小值，而是分解成計(jì)算量較小的極小化處理[26]，計(jì)算速度快且較為穩(wěn)定。

采集得到的數(shù)據(jù)本質(zhì)上是由兩條二維剖面數(shù)據(jù)組成，目前針對(duì)二維大地電磁數(shù)據(jù)主要還是使用一維、二維反演計(jì)算手段，常規(guī)大地電磁二維反演利用不同極化模式的數(shù)據(jù)可以獲得多種不同的反演結(jié)果，各個(gè)結(jié)果往往相差較大。前人對(duì)于極化模式的選擇已經(jīng)做了許多工作，認(rèn)為TE極化模式對(duì)于淺部的三維異常體易產(chǎn)生低阻假異常，而TM極化模式可得到較好的效果[27]。但TM極化模式與TE極化模式相比更易受到旁側(cè)異常體的影響，體現(xiàn)的是測(cè)線下方及附近區(qū)域的宏觀電性特征[28]。TETM聯(lián)合模式反演結(jié)果的穩(wěn)定性又弱于TM極化模式，效果也只是一般[29]。對(duì)于這種二維剖面數(shù)據(jù)如果選用三維反演效果則會(huì)好很多，可以獲得測(cè)線附近較可靠的三維地電結(jié)構(gòu)，由于數(shù)據(jù)的限制仍不能夠獲得較遠(yuǎn)處的地質(zhì)體信息，這在許多實(shí)測(cè)資料中得到了驗(yàn)證[30]。

圖2 316號(hào)點(diǎn)全阻抗數(shù)據(jù)(具體位置見圖1)Fig.2 Full impedance data at point 316

反演選用1～10 000 Hz的頻率范圍，在所有頻點(diǎn)中等對(duì)數(shù)間隔抽出24個(gè)頻點(diǎn)，測(cè)點(diǎn)共計(jì)69個(gè)，合計(jì)頻率數(shù)為1 656個(gè)。反演初始模型為200 Ω·m的均勻半空間；水平x和y方向的網(wǎng)格大小均為40 m，兩個(gè)方向均有100個(gè)網(wǎng)格；x、y、z3個(gè)方向的擴(kuò)邊網(wǎng)格都為5個(gè)，x、y方向網(wǎng)格尺寸按照2.5的冪次擴(kuò)大；z方向的首層厚度是15 m，層厚以1.1系數(shù)遞增，單方向以5的系數(shù)擴(kuò)邊5個(gè)網(wǎng)格，一共30層，最深達(dá)到17 km，誤差門限設(shè)置為5%的噪聲值，反演數(shù)據(jù)類型選為主阻抗張量。

經(jīng)過120次的迭代計(jì)算，擬合差由203.48降到了6.81。圖3a、b是反演迭代中擬合差和λ的變化，擬合差開始下降較快，然后趨于平穩(wěn)。圖3b為λ變化曲線，當(dāng)擬合差變化量小于閾值時(shí)，降低λ，降低模型擬合在目標(biāo)函數(shù)中的權(quán)重。圖3c為本次反演的數(shù)據(jù)擬合交匯圖，理論上，如果完全擬合數(shù)據(jù)，反演數(shù)據(jù)應(yīng)該等于或者非常接近觀測(cè)數(shù)據(jù)，在交匯圖上應(yīng)為一條45°的直線。擬合越好，這種規(guī)律越好。圖3c中數(shù)據(jù)基本沿45°直線分布。因此，本次三維反演的數(shù)據(jù)擬合效果較好。

(a) 擬合差變化曲線; (b) 正則化參數(shù)λ變化曲線; (c) 擬合數(shù)據(jù)交匯圖。圖3 三維反演擬合結(jié)果評(píng)價(jià)Fig.3 Evaluation of fitting results by 3D inversion

3 結(jié)果與討論

3.1 反演結(jié)果

為了更清晰地表達(dá)深部的信息，將三維反演結(jié)果按照剖面位置分別做出其垂直斷面圖(圖4)。圖4a中測(cè)線3下方，在縱向上呈現(xiàn)明顯的三層電阻率結(jié)構(gòu)，第一層為高阻，第二層為低阻，第三層為中高阻。其中淺部100 m內(nèi)顯示的高阻特征，推測(cè)為淺成侵入巖或變質(zhì)片麻巖的體現(xiàn)。低阻異常在點(diǎn)號(hào)236、284和364--412的下方，在500 m及更深處呈層狀分布。第三層的中高阻可能是研究區(qū)的花崗閃長(zhǎng)巖基底。

圖4b給出了測(cè)線2下方電阻率斷面，可以明顯看出存在著與圖4a相似的電阻率分布。 在其中點(diǎn)號(hào)60--110、140--290和320--380下方200～600 m的深部中含有電阻率在10 Ω·m以下的低阻異常，推測(cè)這些低阻異常是鐵礦及礦化圍巖的反映，并且在500 m處這些低阻異常相連接，表現(xiàn)為水平狀的低阻層，厚度在100 m±。

(a)三維反演結(jié)果沿測(cè)線3電阻率斷面;(b)三維反演結(jié)果沿測(cè)線2電阻率斷面;(c)測(cè)線2二維反演電阻率斷面。涂黑的倒三角代表兩條測(cè)線交點(diǎn)的位置。圖4 電性結(jié)構(gòu)模型切片圖Fig.4 Section diagrams of electrical structure model

對(duì)測(cè)線2做二維非線性共軛梯度法反演進(jìn)行對(duì)比。由于影響二維反演結(jié)果的參數(shù)設(shè)置非常多，所得到的二維反演結(jié)果也大不相同[31--32]。 在對(duì)照工區(qū)的地質(zhì)地球物理認(rèn)識(shí)和三維反演結(jié)果的基礎(chǔ)上，最終選用TETM聯(lián)合模式的二維反演結(jié)果，這也是二維反演中綜合TE和TM極化模式優(yōu)點(diǎn)的較好選擇(圖4c)。

比較三維結(jié)果與二維結(jié)果可以看出，在三維反演結(jié)果中，中心區(qū)域低阻異常的位置及電阻率能夠較好的與二維反演結(jié)果相對(duì)應(yīng)；但在320--380號(hào)點(diǎn)下方的低阻體在二維反演剖面中變得更深，達(dá)到了1 km，這可能是二維反演中將測(cè)線附近的低阻異常響應(yīng)劃分為測(cè)線下方更大規(guī)模低阻異常的反映；對(duì)于高阻體來說，在二維反演剖面中最下層的高阻體比三維反演結(jié)果中的高阻體電阻率高，超過2 000 Ω·m，研究區(qū)的基底是電阻率1 000 Ω·m±的花崗閃長(zhǎng)巖，三維結(jié)果中深部是較可信的，而二維反演在這里的幾個(gè)點(diǎn)低頻段擬合并不是很好，造成了偏高的電阻率；另一個(gè)明顯區(qū)別是在二維結(jié)果大號(hào)點(diǎn)下方淺部是高阻體而深部是低阻，這與三維反演結(jié)果完全相反。

三維反演與二維反演結(jié)果之間存在差異是正常的，因?yàn)閷?shí)際的礦體更適合三維模型描述，而二維模型忽略了實(shí)際數(shù)據(jù)的三維性，將剖面旁側(cè)的異常效應(yīng)強(qiáng)行用剖面下方的異常進(jìn)行平衡，造成反演結(jié)果與實(shí)際模型的偏差，從而在500 m深處出現(xiàn)低阻異常，而這個(gè)異?？赡苁遣豢尚诺模欢鴾\部的低阻異常附近可能有電阻率較高的高阻體，造成二維反演結(jié)果中異常位置和電阻率的不準(zhǔn)確。

3.2 討論與分析

兩條剖面相交處的低阻異常與義唐村航磁異常位置相近，航磁異常強(qiáng)度在500 nT±，該異常位于新絳至河津熱液型鐵礦成礦遠(yuǎn)景區(qū)內(nèi)，周邊已完成的鉆孔中在739～1 859 m范圍內(nèi)發(fā)現(xiàn)了多層磁鐵礦體[33]。推斷反演結(jié)果的電阻率模型中低阻異常為磁鐵礦及磁鐵礦化花崗閃長(zhǎng)巖引起的。

通過對(duì)鉆孔中各類巖性物性的測(cè)定統(tǒng)計(jì)，鐵礦體是低電阻率的地質(zhì)體，電阻率平均值在10 Ω·m±，而圍巖及沉積蓋層電阻率均在1 000 Ω·m以上。所以為了更直觀地展示鐵礦的三維形態(tài)，拋去磁鐵礦化圍巖的影響，用10 Ω·m的電阻率等值面反映測(cè)鐵礦的立體圖像。由于淺部異常分辨力不足，所以將地表附近約50 m的信息舍棄，得到的異常體形態(tài)如圖5所示。這也是富存鐵礦的有利區(qū)域。

電性結(jié)果顯示，除了剖面中心500 m處的層狀異常，在3號(hào)剖面的東端以及2號(hào)剖面的北端，可能在剖面的外側(cè)存在更大規(guī)模的異常，并且延伸到1 km或者更深處，這是二維反演中并沒有發(fā)現(xiàn)的。在平面上該異常位于義唐村的東北部，在航磁探測(cè)中這里屬于航磁異常的邊緣，沒有引起大的磁異常，這樣的深部電性異常似乎更值得重視，極可能是深部鐵礦的體現(xiàn)。

圖5 三維反演得到的磁鐵礦體空間分布圖(異常體由10 Ω·m電阻率等值面構(gòu)成)Fig.5 Spatial distribution map of iron deposit obtained by 3D inversion

4 結(jié)論

(1) 相對(duì)于二維反演，三維反演可以將局部異?？臻g分布特征刻畫得更客觀，并且能克服剖面旁側(cè)局部電阻率異常對(duì)觀測(cè)數(shù)據(jù)的影響。

(2) 測(cè)區(qū)地下電阻率特征呈明顯的3層結(jié)構(gòu)，第一層高阻反映淺成侵入巖或變質(zhì)片麻巖，下方的低阻可能指示鐵致異常，而高阻的第三層則反映了研究區(qū)花崗閃長(zhǎng)巖基底的特征。

(3) 圈定了研究區(qū)剖面域鐵礦及礦化圍巖的有利分布區(qū)，主要分布于測(cè)區(qū)地下500 m深度附近，厚度200 m±。

(4) 發(fā)現(xiàn)了新的電阻率異常，異常位于義唐村的東北部，埋深超過1 km，可能是深部鐵礦的體現(xiàn)。