• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于交叉梯度交替結(jié)構(gòu)約束的二維地震走時(shí)與全通道直流電阻率聯(lián)合反演

    2016-11-23 05:59:50高級(jí)張海江
    地球物理學(xué)報(bào) 2016年11期
    關(guān)鍵詞:斑巖

    高級(jí),張海江

    1 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實(shí)驗(yàn)室,合肥 2300262 安徽萬泰地球物理技術(shù)有限公司,合肥 230026

    林方麗1,2,王光杰1,楊曉勇3

    1 中國科學(xué)院地質(zhì)與地球物理研究所,北京 1000292 中國科學(xué)院大學(xué),北京 1000493 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院,合肥 230026

    ?

    基于交叉梯度交替結(jié)構(gòu)約束的二維地震走時(shí)與全通道直流電阻率聯(lián)合反演

    高級(jí)1,2,張海江1*

    1 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院 地震與地球內(nèi)部物理實(shí)驗(yàn)室,合肥 2300262 安徽萬泰地球物理技術(shù)有限公司,合肥 230026

    在利用不同的地球物理勘探方法對(duì)地下復(fù)雜介質(zhì)成像時(shí),因觀測系統(tǒng)的非完備性及數(shù)據(jù)本身對(duì)某些巖石物性的不敏感性,單獨(dú)成像的結(jié)果存在較大的不確定性和不一致性.對(duì)于地震體波走時(shí)成像與直流電阻率成像,均面臨著成像陰影區(qū)問題.對(duì)于地震走時(shí)成像,地震射線對(duì)低速區(qū)域覆蓋較差形成陰影區(qū),造成低速區(qū)域分辨率降低.對(duì)于電阻率成像,電場線在高阻區(qū)域分布較少,造成高阻區(qū)域分辨率較低.為了提高地下介質(zhì)成像的精度,Gallado和Meju(2003)提出了基于交叉梯度結(jié)構(gòu)約束的聯(lián)合地球物理成像方法.在要求不同的物性模型擬合各自對(duì)應(yīng)的數(shù)據(jù)同時(shí),模型之間的結(jié)構(gòu)要求一致,即交叉梯度趨于零.為了更有效地實(shí)現(xiàn)基于交叉梯度的結(jié)構(gòu)約束,我們提出了一種新的交替結(jié)構(gòu)約束的聯(lián)合反演流程,即交替反演不同的數(shù)據(jù)而且在反演一種數(shù)據(jù)時(shí)要求對(duì)應(yīng)的模型與另一個(gè)模型結(jié)構(gòu)一致.新的算法能夠更容易地把單獨(dú)的反演系統(tǒng)耦合在一起,而且也更容易建立結(jié)構(gòu)約束和數(shù)據(jù)擬合之間的平衡.基于新的聯(lián)合反演流程,我們測試了基于交叉梯度結(jié)構(gòu)約束的二維跨孔地震走時(shí)和直流電阻率聯(lián)合成像.合成數(shù)據(jù)測試表明,我們提出的交替結(jié)構(gòu)約束流程能夠很好地實(shí)現(xiàn)基于交叉梯度結(jié)構(gòu)約束的聯(lián)合成像.與單獨(dú)成像結(jié)果相比,地震走時(shí)和全通道電阻率聯(lián)合成像更可靠地確定了速度和電阻率異常.

    體波走時(shí)成像;電阻率成像;結(jié)構(gòu)相似性;交叉梯度;聯(lián)合反演

    1 引言

    不同的地球物理勘探方法對(duì)地下介質(zhì)的成像各有其優(yōu)缺點(diǎn),而且都存在某種程度的不確定性和多解性(胡祥云等,2006).引起地球物理勘探不確定性和多解性的原因包括:有限的觀測數(shù)據(jù)、由于環(huán)境干擾和儀器精度存在的觀測誤差、地球物理正演方程和真實(shí)物理過程存在的偏差、不同地質(zhì)模型觀測的地球物理場相同情況的存在、反演算法存在的偏差等因素.因此為了更可靠地刻畫地下介質(zhì),多種地球物理數(shù)據(jù)之間的聯(lián)合反演解釋就成為一種有效的途徑(Gallardo and Meju,2003,2004;Gallardo et al.,2012;Abubakar et al.,2012).

    地下介質(zhì)具有不同的物理性質(zhì),如密度、電阻率、極化率和磁性.一類的聯(lián)合反演是基于多屬性參數(shù)的耦合即建立地下介質(zhì)不同物理性質(zhì)之間的經(jīng)驗(yàn)關(guān)系 (如Berge et al.,2000;于鵬等,2006).例如因?yàn)樗俣群兔芏戎g存在較好的經(jīng)驗(yàn)關(guān)系,Maceira和Ammon (2009)發(fā)展了面波和重力的聯(lián)合成像方法.但是對(duì)于某些地球物理屬性,例如地震速度和電阻率,很難建立起二者之間普遍適用的經(jīng)驗(yàn)關(guān)系.在這種情況下,發(fā)展了基于結(jié)構(gòu)一致性的另一類聯(lián)合反演方法(如,Lines and Lines,1988;Haber and Oldenburg,1997;Musil et al.,2003;Gallardo-Delgado et al.,2003;Hu et al.,2009).上述兩類聯(lián)合反演方法均存在一定的問題(王俊,2014):第一類方法需要預(yù)先知道不同物性參數(shù)基于巖石物理學(xué)的聯(lián)系關(guān)系才能實(shí)現(xiàn)聯(lián)合反演,而多數(shù)情況下兩種屬性之間的聯(lián)系難以確定;對(duì)于第二類方法,如何解決模型結(jié)構(gòu)的融合一直是制約該方法發(fā)展的重大障礙.Gallardo和Meju(2003,2004)提出了基于模型交叉梯度的結(jié)構(gòu)約束實(shí)現(xiàn)直流電法與地震走時(shí)的聯(lián)合反演,取得了較好的效果.該方法利用梯度場來刻畫模型的結(jié)構(gòu),如果兩個(gè)模型的梯度一致即結(jié)構(gòu)也一致.在反演的過程中,一方面要求各個(gè)模型擬合數(shù)據(jù),同時(shí)通過改變模型讓二者之間的交叉梯度趨于零實(shí)現(xiàn)結(jié)構(gòu)的一致.該聯(lián)合反演方法無需明確模型之間的物性關(guān)系,目前得到了地球物理學(xué)界廣泛的關(guān)注.例如,Abubakar等(2012)進(jìn)一步結(jié)合巖石物理性質(zhì)及結(jié)構(gòu)相似性實(shí)現(xiàn)了大地電磁與面波的聯(lián)合反演,Bennington等(2015)發(fā)展了基于歸一化交叉梯度約束的三維天然地震和二維大地電磁的聯(lián)合反演算法,李兆祥等(2015)通過交叉梯度約束實(shí)現(xiàn)了電阻率與極化率的聯(lián)合反演,并取得了較好的效果.

    對(duì)于不同地球物理數(shù)據(jù)的反演,模型的離散化和數(shù)據(jù)對(duì)模型的敏感度是不一樣的,因此把對(duì)應(yīng)于不同數(shù)據(jù)的反演系統(tǒng)與交叉梯度結(jié)構(gòu)約束完全耦合成一個(gè)大的聯(lián)合反演系統(tǒng)是一個(gè)復(fù)雜的具有挑戰(zhàn)性的工程(Gallardo and Meju,2007).為了降低實(shí)現(xiàn)基于交叉梯度結(jié)構(gòu)約束的難度,Gallardo和Meju(2007)提出了分別進(jìn)行數(shù)據(jù)擬合和交叉梯度結(jié)構(gòu)約束的聯(lián)合反演流程,通過迭代交替進(jìn)行數(shù)據(jù)反演和交叉梯度反演實(shí)現(xiàn)聯(lián)合反演.雖然該流程降低了實(shí)現(xiàn)聯(lián)合反演的難度,但是由于數(shù)據(jù)擬合與結(jié)構(gòu)約束完全分開,因此迭代交替反演時(shí)很難平衡數(shù)據(jù)擬合和結(jié)構(gòu)約束的程度.例如在一次迭代中,利用不同的數(shù)據(jù)首先對(duì)各自對(duì)應(yīng)的模型進(jìn)行了更新,但再用交叉梯度結(jié)構(gòu)約束對(duì)兩個(gè)模型進(jìn)行改變使二者結(jié)構(gòu)一致時(shí),可能會(huì)導(dǎo)致對(duì)數(shù)據(jù)的擬合在很大程度上變差;另外利用數(shù)據(jù)反演對(duì)模型進(jìn)行更新之后,二者的結(jié)構(gòu)一致性可能會(huì)得到破壞.所以利用Gallardo和Meju(2007)提出的聯(lián)合反演流程,因?yàn)閿?shù)據(jù)擬合和結(jié)構(gòu)約束是分開進(jìn)行的,因此如何平衡二者對(duì)模型的更新是一個(gè)挑戰(zhàn).

    針對(duì)于此問題,本文提出了一種新的基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演流程,避免了Gallardo和Meju(2007)提出的數(shù)據(jù)擬合和結(jié)構(gòu)約束完全分開的缺點(diǎn),同時(shí)還避免了把對(duì)應(yīng)于不同數(shù)據(jù)的反演系統(tǒng)放在一起的困難.在我們新提出的反演流程中,對(duì)一種數(shù)據(jù)進(jìn)行反演的時(shí)候,同時(shí)要求對(duì)應(yīng)的模型與另外一個(gè)模型結(jié)構(gòu)一致,兩種數(shù)據(jù)反演交替進(jìn)行直至收斂.結(jié)構(gòu)的約束還是通過交叉梯度實(shí)現(xiàn),但不同的是只改變一個(gè)模型而另一個(gè)模型不變使二者的結(jié)構(gòu)一致.本文首先研究了全通道2.5D跨孔直流電阻率層析成像(潘紀(jì)順等,2010;潘克家等,2012;雷旭友等,2009;Loke and Barker,1996;Loke and Dahlin,2002;Loke,2003),并比較了全通道接收與傳統(tǒng)四極法采集反演結(jié)果的不同.其次介紹了基于MSFM(Multi-stencils Fast Marching Methods)射線追蹤方法的跨孔地震走時(shí)反演.然后提出了新的基于交叉梯度交替結(jié)構(gòu)約束的聯(lián)合直流電阻率和地震走時(shí)反演流程.最后基于合成數(shù)據(jù),驗(yàn)證了本文提出的聯(lián)合反演流程的有效性.

    2 全通道2.5D跨孔直流電阻率層析成像

    本文在實(shí)現(xiàn)電阻率正反演時(shí),首先采用有限差分方法求解2.5D靜電場方程(范翠松等,2012;劉斌等,2012),求取觀測位置處的電位值,并基于伴隨矩陣(Adjoint method)方法計(jì)算非線性靈敏度矩陣(Pidlisecky et al.,2007;吳小平和徐果明,2000;Mcgillivray and Oldenburg,1990),最后利用非線性牛頓共軛梯度反演方法(劉云鶴和殷長春,2013)實(shí)現(xiàn)全通道電阻率層析成像(吳小平等,2015).

    二維跨孔直流電阻率層析成像裝置形式為在兩個(gè)鉆孔內(nèi)布置供電與接收電極.其供電形式可以采用單極、偶極、三極、四極等裝置形式,跑極方式與地表采集相同,也可以采用全通道接收形式.本文在實(shí)現(xiàn)直流電阻率與地震波走時(shí)聯(lián)合反演時(shí)采用全通道采集方式.

    圖1為跨孔電阻率裝置及正反演結(jié)果圖.其中圖1a為電極布置方式,分別在兩個(gè)鉆孔中布置32個(gè)供電與采集電極,孔距15 m;圖1b為電阻率模型,其中背景電阻率為100 Ωm,低高阻電阻率異常值分別為10 Ωm和1000 Ωm,異常大小為4 m×4 m;圖1c為兩個(gè)電極供電時(shí)觀測區(qū)域2D電場分布等值線圖,從中可以看出跨孔電阻率電場在空間分布較為均勻,電場線穿過孔間觀測區(qū)域,具有對(duì)孔間異常進(jìn)行成像的基礎(chǔ);圖1d為電阻率反演結(jié)果,從反演結(jié)果圖中可以看出高阻及低阻異常均具有較好的成像及較少的假異常,但受到電法勘探體積效應(yīng)的影響反演結(jié)果分辨率相對(duì)較低.

    3 地震體波初至走時(shí)成像

    地震體波走時(shí)層析成像是利用初至波到時(shí)構(gòu)建地下速度結(jié)構(gòu)的一種重要的方法(Zelt and Smith,1992),該方法可以利用主動(dòng)源、被動(dòng)源或同時(shí)使用兩種震源 (Di Stefano and Chiarabba 2002;Ritter et al.,2001).對(duì)于體波走時(shí)成像,需要追蹤震源和檢波器之間的射線和計(jì)算對(duì)應(yīng)的走時(shí).射線追蹤方法主要包括基于程函方程(Vidale,1988,1990)及基于射線方程(Zelt and Smith,1992;張風(fēng)雪等,2010)方法.本文采用不等間距節(jié)點(diǎn)法對(duì)模型區(qū)域網(wǎng)格進(jìn)行剖分(Sambridge,1990;Sambridge and Rawlinson,2004),走時(shí)計(jì)算采用有限差分法(Vidale,1988,1990)求解程函方程的MSFM(Multi-Stencils Fast Marching Methods)方法(Osher and Fedkiw,2003;Popovici and Sethian,2002).反演時(shí)為提高反演結(jié)果的穩(wěn)定性,采用二階吉洪諾夫正則化(Tikhonov,1977)約束(Calvetti et al.,2000),并采用阻尼最小二乘LSQR法迭代求解(Paige and Saunders,1982).

    圖1 二維跨孔全通道電阻率成像(a) 電極分布;(b) 真實(shí)電阻率模型;(c) 二維電位場,V為電位,I為電流;(d) 電阻率反演結(jié)果.Fig.1 Two-dimensional cross borehole full channel resistivity tomography(a) Distribution of electrodes;(b) Synthetic resistivity model;(c) 2D electric field,V is electric potential and I is electric current;(d) Inverted resistivity model.

    圖2為跨孔地震波走時(shí)觀測系統(tǒng)及正反演結(jié)果圖,對(duì)應(yīng)于圖1所示的跨孔電阻率觀測系統(tǒng).其中圖2a為觀測系統(tǒng)圖,檢波器和震源分別放置在兩個(gè)鉆孔中,每個(gè)孔中放置16炮30個(gè)檢波器,炮距2 m、檢波距1 m,在一個(gè)孔中放炮時(shí)另一個(gè)孔中30個(gè)檢波器同時(shí)接收;圖2b為速度模型,其中背景速度為2000 m·s-1、異常速度值分別為2500 m·s-1和1500 m·s-1、異常大小為4 m×4 m,其中低速對(duì)應(yīng)低阻而高速對(duì)應(yīng)高阻異常;圖2c為基于MSFM求解的波前走時(shí)場等值線圖;圖2d為走時(shí)反演結(jié)果,從反演結(jié)果圖中可以看出走時(shí)成像能確定高低速異常體所在位置,與電阻率成像結(jié)果(圖1d)比較,具有較高的縱橫向分辨率,但在異常周圍產(chǎn)生較多的假異常,產(chǎn)生假異常的主要原因是模型測試采用了符合實(shí)際工程施工的跨孔觀測系統(tǒng),即炮點(diǎn)和接收點(diǎn)位于兩個(gè)孔中,因此觀測系統(tǒng)覆蓋角度具有一定的局限性,導(dǎo)致反演的速度模型在x方向的分辨率低于z方向分辨率.

    4 基于交叉梯度約束的聯(lián)合反演方法

    基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演由Gallardo和Meju (2003)首次提出,后續(xù)被廣泛應(yīng)用于地震數(shù)據(jù)和大地電磁數(shù)據(jù)、地震與直流電法數(shù)據(jù)的聯(lián)合反演(Gallardo and Meju,2003,2004,2007;Gallardo et al.,2005,2012).國內(nèi)基于交叉梯度約束的聯(lián)合反演研究也逐步展開(彭淼,2012;周麗芬,2012;王俊,2014).對(duì)于三維地震速度模型ms和電阻率模型mr,二者之間的交叉梯度定義為速度梯度和電阻率梯度的叉積:

    =txi+tyj+tzk.

    (1)

    在二維模型情況下,模型沿y軸方向變化為0,模型之間沿x與z方向變化的交叉梯度只存在y方向分量:

    (2)

    圖3顯示了兩個(gè)不同二維模型如何計(jì)算交叉梯度.從圖中可以看出,當(dāng)兩個(gè)模型梯度方向相同或相反時(shí),即結(jié)構(gòu)一致時(shí)交叉梯度值為零,反之不為零.兩個(gè)模型交叉梯度的這個(gè)性質(zhì)為利用交叉梯度約束反演奠定了基礎(chǔ),即在反演時(shí)要求速度與電阻率模型交叉梯度值趨于零,即模型之間具有一致的結(jié)構(gòu).因此基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演要求不同的模型在擬合對(duì)應(yīng)的數(shù)據(jù)同時(shí)它們之間的結(jié)構(gòu)要一致.另外考慮到聯(lián)合反演系統(tǒng)存在一定的病態(tài)性,因此對(duì)模型加入一定的平滑約束,據(jù)此建立的聯(lián)合反演目標(biāo)函數(shù)為(Gallardo and Meju,2003;Bennington et al.,2015):

    圖2 二維跨孔地震走時(shí)層析成像(a) 兩個(gè)孔中震源和檢波器的分布;(b) 真實(shí)速度模型;(c) 單炮地震波走時(shí)場等值線分布;(d) 地震波走時(shí)反演速度模型.Fig.2 Two-dimensional cross borehole seismic travel time tomography(a) Distribution of sources and receivers in two boreholes;(b) Real velocity model;(c) Travel time field using the MSFM ray tracing method;(d) Inversion velocity model.

    圖3 兩個(gè)模型交叉梯度的測試(a) 模型1;(b) 模型1的梯度場;(c) 模型2;(d) 模型2的梯度場;(e) 模型1與模型2的交叉梯度.Fig.3 Test of cross gradient for two models(a) Model 1;(b) Gradient field of model 1;(c) Model 2;(d) Gradient field of model 2;(e) Cross gradients between model 1 and model 2.

    (3)

    式中:fs(ms)為根據(jù)地震速度模型ms計(jì)算的理論地震走時(shí)數(shù)據(jù),ds為觀測走時(shí)數(shù)據(jù),μs為地震走時(shí)擬合權(quán)重,Cs為觀測走時(shí)數(shù)據(jù)誤差的協(xié)方差矩陣,Ls為速度模型平滑因子,λs為速度模型平滑權(quán)重;fr(mr)為根據(jù)電阻模型mr計(jì)算的理論電位數(shù)據(jù),dr為觀測電位數(shù)據(jù),μr為電位數(shù)據(jù)擬合權(quán)重,Cr為觀測電位數(shù)據(jù)誤差的協(xié)方差矩陣,Lr為電阻率模型平滑因子,λr為電阻率模型平滑權(quán)重.t(ms,mr)為速度模型與電阻率模型之間的交叉梯度,β為交叉梯度約束權(quán)重.

    對(duì)于交叉梯度函數(shù)t(ms,mr),我們可以利用一階泰勒展開的方法進(jìn)行線性化,寫成矩陣形式為:

    (4)

    (5)

    圖4 Galldado和Meju(2007)提出的基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演流程圖Fig.4 Flow chart of joint inversion based on cross-gradient structure constraint proposed by Gallardo and Meju (2007)

    Bennington等(2015)進(jìn)行三維地震走時(shí)和二維大地電磁聯(lián)合反演時(shí),就是基于公式(5)的反演矩陣系統(tǒng),把地震走時(shí)反演和大地電磁反演以及交叉梯度結(jié)構(gòu)約束真正融合在一起.該種方式雖然在數(shù)學(xué)上具有嚴(yán)格性,但在實(shí)際編程實(shí)現(xiàn)時(shí)具有挑戰(zhàn)性,因?yàn)榘褟?fù)雜的單個(gè)反演系統(tǒng)融合成為一個(gè)大的反演系統(tǒng)容易出現(xiàn)問題.Gallardo和Meju(2007)提出了一種解決方案,避免了融合不同反演系統(tǒng)存在的挑戰(zhàn)性,如圖4所示的反演流程圖.例如對(duì)于地震走時(shí)和電阻率的聯(lián)合反演,在Gallardo和Meju(2007)提出的反演流程中,首先基于獨(dú)立的地震走時(shí)成像和電阻率成像系統(tǒng)分別利用地震走時(shí)數(shù)據(jù)和電位數(shù)據(jù)進(jìn)行反演得到新的地震速度模型和電阻率模型,然后利用交叉梯度結(jié)構(gòu)約束對(duì)兩個(gè)模型進(jìn)行更新得到新的模型,接著再重復(fù)上述的過程直至擬合數(shù)據(jù)和二個(gè)模型結(jié)構(gòu)一致.這種聯(lián)合反演方式把結(jié)構(gòu)約束與數(shù)據(jù)擬合孤立地進(jìn)行,雖然可以直接利用現(xiàn)有的反演系統(tǒng)進(jìn)行聯(lián)合反演,但是造成結(jié)構(gòu)約束權(quán)重難以選擇,即結(jié)構(gòu)約束太強(qiáng)會(huì)導(dǎo)致數(shù)據(jù)擬合發(fā)散,若結(jié)構(gòu)約束太弱會(huì)導(dǎo)致其失去作用,難以選取合適的交叉梯度結(jié)構(gòu)約束權(quán)重因子.

    本文為克服此缺陷,提出一種新的交叉梯度結(jié)構(gòu)約束與數(shù)據(jù)擬合進(jìn)行聯(lián)合反演的方式,避免了交叉梯度結(jié)構(gòu)約束及兩種數(shù)據(jù)擬合之間權(quán)重選擇困難的問題.式(5)所示的聯(lián)合反演矩陣系統(tǒng)可以分開為兩個(gè)交替進(jìn)行的子反演系統(tǒng):

    (6)

    (7)

    圖5 本文提出的基于交叉梯度交替結(jié)構(gòu)約束的聯(lián)合反演流程圖Fig.5 Flow chart of joint inversion with alternating cross-gradient based structure constraint proposed in this study

    圖5給出了本文提出的基于交叉梯度結(jié)構(gòu)約束的地震走時(shí)和電阻率聯(lián)合反演流程:

    (1) 輸入初始速度與電阻率模型并計(jì)算二者之間的交叉梯度.

    (2) 基于公式(6)和(7)分別進(jìn)行基于交叉梯度結(jié)構(gòu)約束的地震走時(shí)和電阻率成像.

    (3) 更新速度與電阻率模型.

    (4) 計(jì)算兩個(gè)模型的交叉梯度以及地震走時(shí)數(shù)據(jù)和電位數(shù)據(jù)的殘差.

    (5) 判斷數(shù)據(jù)擬合和交叉梯度值是否滿足迭代終止條件,若滿足直接輸出結(jié)果,若不滿足跳回步驟(2).一般迭代終止條件包括,迭代次數(shù)、每次迭代殘差下降的梯度、殘差剩余量占初始?xì)埐畎俜直取⒛P透铝康?本文迭代終止及收斂條件為電阻率或地震走時(shí)殘差下降量低于某一百分比時(shí)迭代終止(如迭代殘差下降量小于1%).

    5 合成數(shù)據(jù)測試

    我們利用圖1和圖2所建立的地震速度與電阻率模型以及對(duì)應(yīng)的觀測系統(tǒng),對(duì)本文提出的基于交叉梯度的交替結(jié)構(gòu)約束聯(lián)合反演算法進(jìn)行了測試.對(duì)于公式(6)和(7)中出現(xiàn)的平滑系數(shù)λ,阻尼因子ε和數(shù)據(jù)相對(duì)權(quán)重β可以采用L-curve方式進(jìn)行選擇(Hansen and O′Leary,1993;Aster et al.,2005).例如對(duì)于基于電阻率模型結(jié)構(gòu)約束的地震走時(shí)速度成像,通過L-curve分析,得到對(duì)應(yīng)于曲線拐點(diǎn)的最優(yōu)權(quán)重因子為4×106,即在數(shù)據(jù)擬合和結(jié)構(gòu)一致性方面達(dá)到了平衡,如圖6所示.

    圖7顯示了地震走時(shí)及電阻率聯(lián)合反演得到的速度模型和電阻率模型.與單獨(dú)反演得到的速度模型比較,聯(lián)合反演得到的速度模型很大程度上消除了低速和高速異常周圍出現(xiàn)的由于拖尾效應(yīng)引起的假異常.對(duì)于電阻率模型,聯(lián)合反演得出的異常中心位置更靠近真實(shí)模型,尤其是低阻異常,而且反演結(jié)果數(shù)值更接近真實(shí)模型數(shù)值.

    圖8為地震走時(shí)及電阻率單獨(dú)反演及聯(lián)合反演走時(shí)殘差和電位殘差隨著迭代次數(shù)的變化曲線.對(duì)于聯(lián)合反演和單獨(dú)反演來說,電位殘差隨著迭代次數(shù)變化的趨勢基本一致,但最終聯(lián)合反演得出的電阻率模型能更好地?cái)M合數(shù)據(jù),也對(duì)應(yīng)著恢復(fù)更好的電阻率模型.相比較而言,聯(lián)合反演得到的速度模型對(duì)地震走時(shí)數(shù)據(jù)的擬合比單獨(dú)反演變差,最終模型對(duì)應(yīng)的走時(shí)殘差變大.從速度模型(圖7c)可以看出,由于地震走時(shí)反演系統(tǒng)的病態(tài)性,單獨(dú)地震走時(shí)反演過度地?cái)M合了數(shù)據(jù),導(dǎo)致在真正的異常周圍出現(xiàn)很多的假象,但利用聯(lián)合反演,速度模型的結(jié)構(gòu)得到了電阻率模型的約束,因此速度模型不能過于“自由”的變化.對(duì)于本理論模型測試,可以看出電阻率結(jié)構(gòu)起了相對(duì)主導(dǎo)的作用.

    圖6 交叉梯度與擬合誤差L-curve曲線Fig.6 Cross-gradient and fitting error L-curve

    圖7 單獨(dú)反演和聯(lián)合反演得出的速度和電阻率比較(a) 真實(shí)速度模型;(b) 真實(shí)電阻率模型;(c) 單獨(dú)走時(shí)反演得到的速度模型;(d) 單獨(dú)電阻率反演得到的電阻率模型;(e) 速度和電阻率聯(lián)合反演得到的速度模型;(f) 速度和電阻率聯(lián)合反演得到的電阻率模型.Fig.7 Comparison of seismic velocity and electrical resistivity from separate and joint inversions(a) True seismic velocity model;(b) True electrical resistivity model;(c) Seismic velocity model from separate travel time inversion;(d) Resistivity model from separate resistivity tomography;(e) Seismic velocity model from joint inversion;(f) Resistivity model from joint inversion.

    圖9比較了單獨(dú)反演與聯(lián)合反演速度和電阻率模型的交叉梯度.從圖中可以看出,單獨(dú)反演模型的交叉梯度值分布范圍為-5.0~+5.0,聯(lián)合反演后模型交叉梯度值分布范圍為-0.5~+0.5,即聯(lián)合反演后模型之間交叉梯度值下降到約為單獨(dú)反演的10%(圖9a,9b),同時(shí)從整個(gè)模型空間看,聯(lián)合反演之后,交叉梯度值都較小,而且數(shù)值分布較均勻(圖9c,9d),說明利用交叉梯度的結(jié)構(gòu)約束使得模型之間具有了更好的結(jié)構(gòu)一致性,符合地質(zhì)模型結(jié)構(gòu)唯一的特性,因此基于結(jié)構(gòu)一致性的聯(lián)合反演提高了反演結(jié)果的合理性.

    在某些地質(zhì)條件下,地震速度和電阻率二者之間可能存在結(jié)構(gòu)不完全一致的情況.例如當(dāng)巖層空隙局部充水時(shí),對(duì)巖層電阻率有較大影響,造成對(duì)應(yīng)區(qū)域的電阻率降低,但可能對(duì)地震波傳播速度的影響較小,造成兩種模型具有不同異常結(jié)構(gòu).為此我們也測試了算法在地震速度和電阻率具有不一致結(jié)構(gòu)情況下的聯(lián)合反演效果,同時(shí)為了檢驗(yàn)本文算法在觀測數(shù)據(jù)存在噪聲時(shí)反演的穩(wěn)定性,我們?cè)谧邥r(shí)觀測數(shù)據(jù)與電位觀測數(shù)據(jù)中分別加入了3%的隨機(jī)誤差.

    圖8 單獨(dú)反演與聯(lián)合反演地震走時(shí)和電位數(shù)據(jù)殘差隨著迭代次數(shù)變化的曲線(a) 地震走時(shí)數(shù)據(jù);(b) 電位數(shù)據(jù).Fig.8 Variation curves of residuals for seismic travel times and electrical potentials along with iteration times for separate and joint inversions(a) Seismic travel time data;(b) Electric potential data.

    圖9 單獨(dú)與聯(lián)合反演速度和電阻率模型交叉梯度的比較(a) 單獨(dú)反演模型交叉梯度在每個(gè)深度上在水平方向的變化;(b)與(a)類似但是為聯(lián)合反演情況;(c) 單獨(dú)反演模型交叉梯度在模型空間的分布;(d) 聯(lián)合反演模型交叉梯度在模型空間的分布.Fig.9 Comparison of cross-gradient variations for seismic velocity and resistivity models from separate and joint inversions(a) Cross-gradient variations along the horizontal direction at each depth for models from separate inversions;(b) Same as (a) but for joint inversion;(c) Distribution of cross-gradients in model space for models from separate inversions;(d) Same as (c) but for joint inversion.

    圖10為地震速度模型與電阻率模型具有不同結(jié)構(gòu)情況下的聯(lián)合反演結(jié)果.地震速度和電阻率模型中的高速和高阻異常具有結(jié)構(gòu)一致性,但是低速和低阻異常的結(jié)構(gòu)不一致.從圖中可以看出,與單獨(dú)反演得到的速度模型比較,聯(lián)合反演得出的速度模型很大程度上消除了高速異常周圍的假異常,同時(shí)速度反演結(jié)果在對(duì)應(yīng)低電阻率異常位置沒有產(chǎn)生假異常.其原因?yàn)楫?dāng)一個(gè)模型的梯度為零,不論另一個(gè)模型的結(jié)構(gòu)如何,交叉梯度也為零.在低電阻率異常區(qū)域,對(duì)應(yīng)的速度模型變化非常小接近常數(shù),梯度變化接近于零.因此在進(jìn)行交叉梯度約束反演時(shí),對(duì)應(yīng)的交叉梯度值亦趨近于零.所以低電阻率異常結(jié)構(gòu)不會(huì)對(duì)速度結(jié)構(gòu)起到約束的作用而導(dǎo)致假異常的出現(xiàn).同樣速度模型中的低速異常結(jié)構(gòu)也沒有在電阻率模型對(duì)應(yīng)區(qū)域帶來假電阻異常.

    在速度模型與電阻率模型具有相同結(jié)構(gòu)的區(qū)域,速度模型與電阻率模型之間得到較好的相互結(jié)構(gòu)約束.高速異常得到高電阻率異常的x方向上的結(jié)構(gòu)約束,使高速異常區(qū)域具有更好的形態(tài)恢復(fù).對(duì)于電阻率模型來說,高電阻率異常區(qū)域得到速度模型中高速異常在z方向結(jié)構(gòu)上的約束,因此聯(lián)合反演提高了電阻率模型高電阻率異常在z方向上的分辨率.

    圖10 速度模型與電阻率模型結(jié)構(gòu)不一致情況下單獨(dú)反演和聯(lián)合反演結(jié)果比較(a) 真實(shí)速度模型;(b) 真實(shí)電阻率模型;(c) 單獨(dú)走時(shí)反演得到的速度模型;(d) 單獨(dú)電阻率反演得到的電阻率模型;(e) 聯(lián)合反演得到的速度模型;(f) 聯(lián)合反演得到的電阻率模型.Fig.10 Comparison of seismic velocity and electrical resistivity models from separate and joint inversions when two models have different structures(a) True seismic velocity model;(b) True electrical resistivity model;(c) Seismic velocity model from separate travel time inversion;(d) Resistivity model from separate resistivity tomography;(e) Seismic velocity model from joint inversion;(f) Resistivity model from joint inversion.

    6 討論和結(jié)論

    多物理屬性的聯(lián)合反演,目標(biāo)是利用地下介質(zhì)的多種巖石物理性質(zhì)來減小模型的多解性(Vozoff and Jupp,1975).不同的地球物理方法對(duì)不同巖石物理性質(zhì)具有不同的分辨率,如電阻率法對(duì)低阻異常區(qū)敏感、速度成像對(duì)高速異常區(qū)敏感,同時(shí)高阻體一般為高速,低速體一般為低阻,因此電阻率成像與速度成像對(duì)介質(zhì)敏感性是互補(bǔ)的.另外不同地球物理數(shù)據(jù)采集時(shí)受不同來源噪聲的影響是不同的,如電阻率法對(duì)電性干擾較為敏感,但地震數(shù)據(jù)受環(huán)境噪聲影響較大.基于上面兩點(diǎn),電阻率成像和地震走時(shí)成像的聯(lián)合反演系統(tǒng)可以更好地刻畫地下介質(zhì).

    針對(duì)完全耦合在一起的聯(lián)合反演系統(tǒng)存在多屬性數(shù)據(jù)數(shù)量級(jí)不同以及不同反演系統(tǒng)耦合在一起存在的技術(shù)實(shí)現(xiàn)問題,我們提出基于交叉梯度交替結(jié)構(gòu)約束的聯(lián)合反演算法.新的反演方法避免了不同屬性數(shù)據(jù)直接耦合在一起的挑戰(zhàn)而且還可以兼顧數(shù)據(jù)擬合和結(jié)構(gòu)約束之間的平衡.基于一個(gè)二維的地震速度和電阻率模型,我們利用新提出的聯(lián)合反演流程進(jìn)行了測試,表明電阻率成像與速度成像是一種較為有效互補(bǔ)的成像方式.通過基于交叉梯度約束的結(jié)構(gòu)一致性聯(lián)合成像,較單獨(dú)反演能夠更好地恢復(fù)速度和電阻率模型.合成模型測試還表明,因?yàn)殡娮杪史囱菘梢缘玫捷^穩(wěn)定平滑模型,因此有助于消除地震走時(shí)成像中出現(xiàn)的次生干擾異常.同時(shí)我們還測試了地震速度和電阻率模型結(jié)構(gòu)不一致的情況,結(jié)果表明基于交叉梯度結(jié)構(gòu)約束的聯(lián)合反演只對(duì)兩個(gè)模型具有一致結(jié)構(gòu)的區(qū)域產(chǎn)生結(jié)構(gòu)約束,而沒有在結(jié)構(gòu)不一致的地方產(chǎn)生假的異常.這主要是因?yàn)樵谝粋€(gè)模型的梯度為零的情況下不論另外一個(gè)模型的異常結(jié)構(gòu)形態(tài)如何變化,兩個(gè)模型的交叉梯度也是零.

    在聯(lián)合反演迭代過程中,要適當(dāng)控制單一模型的收斂速度,即如果一個(gè)模型收斂速度較快使異常周圍模型梯度為零,兩個(gè)模型的交叉梯度在異常周圍也為零,此時(shí)交叉梯度無法起到對(duì)另一模型的結(jié)構(gòu)進(jìn)行結(jié)構(gòu)約束的作用.另外因?yàn)榈卣鹱邥r(shí)和電阻率成像具有不同的收斂性和分辨率,可以采用在聯(lián)合反演迭代第一階段選擇以地震速度模型結(jié)構(gòu)約束為主,提高電阻率反演的分辨率,隨著迭代可以逐漸增大電阻率模型結(jié)構(gòu)約束的權(quán)重,提高地震走時(shí)反演的穩(wěn)定性.

    Aster R C,Borchers B,Thurber C H.2005.Parameter Estimation and Inverse Problems.San Diego:Academic Press.

    Bennington N L,Zhang H J,Thurber C H,et al.2015.Joint inversion of seismic and magnetotelluric data in the Parkfield region of California using the normalized cross-gradient constraint.Pure and Applied Geophysics,172(5):1033-1052.

    Berge P A,Berryman J G,Bertete-Aguirre H,et al.2000.Joint inversion of geophysical data for site characterization and restoration monitoring.LLNL Rep.UCRL-ID-128343,Proj,55411.

    Calvetti D,Morigi S,Reichel L,et al.2000.Tikhonov regularization and the L-curve for large discrete ill-posed problems.Journal of Computational and Applied Mathematics,123(1-2):423-446.

    Di Stefano R,Chiarabba C.2002.Active source tomography at Mt.Vesuvius:Constraints for the magmatic system.Journal of Geophysical Research:Solid Earth,2002,107(B11),doi:10.1029/2001JB000792.

    Fan C S,Li T L,Yan J Y.2012.Research and application experiment on 2.5D SIP inversion.Chinese J.Geophys.(in Chinese),55(12):4044-4050,doi:10.6038/j.issn.0001-5733.2012.12.016.

    Feng R,Li Z M,Li Z W,et al.2004.Resistivity tomography technology.Earthquake Research in China (in Chinese),20(1):13-30.

    Gallardo L A,Meju M A.2003.Characterization of heterogeneous near-surface materials by joint 2D inversion of DC resistivity and seismic data.Geophys.Res.Lett.,30(13):1658.

    Gallardo L A,Meju M A.2004.Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints.Journal of Geophysical Research:Solid Earth,109(B3),doi:10.1029/2003JB002716.

    Gallardo L A,Pérez-Flores M A,Gómez-Trevio E.2005.Refinement of three-dimensional multilayer models of basins and crustal environments by inversion of gravity and magnetic data.Tectonophysics,397(1-2):37-54.

    Gallardo L A,Meju M A.2007.Joint two-dimensional cross-gradient imaging of magnetotelluric and seismic traveltime data for structural and lithological classification.Geophysical Journal International,169(3):1261-1272.

    Gallardo L A,Fontes S L,Meju M A,et al.2012.Robust geophysical integration through structure-coupled joint inversion and multispectral fusion of seismic reflection,magnetotelluric,magnetic,and gravity images:Example from Santos Basin,offshore Brazil.Geophysics,77(5):B237-B251.

    Gallardo-Delgado L A,Pérez-Flores M A,Gómez-Trevio E.2003.A versatile algorithm for joint 3D inversion of gravity and magnetic data.Geophysics,68(3):949-959.

    Gao G Z,Abubakar A,Habashy T M.2012.Joint petrophysical inversion of electromagnetic and full-waveform seismic data.Geophysics,77(3):WA3-WA18.

    Haber E,Oldenburg D.1997.Joint inversion:a structural approach.Inverse Problems,13(1):63-77.

    Hansen P C,O′Leary D P.1993.The use of the L-curve in the regularization of discrete ill-posed problems.SIAM Journal on Scientific Computing,14(6):1487-1503.

    Hu W Y,Abubakar A,Habashy T M.2009.Joint electromagnetic and seismic inversion using structural constraints.Geophysics,74(6):R99-R109.

    Hu X Y,Yang D K,Liu S H,et al.2006.The developing trends of environmental and engineering geophysics.Progress in Geophysics (in Chinese),21(2):598-604.

    Jupp D L B,Vozoff K.1975.Stable iterative methods for the inversion of geophysical data.Geophysical Journal International,42(3):957-976.

    Lei X Y,Li Z W,Zhe J P.2009.Applications and research of the high resolution resistivity method in exploration of caves,mined regions and Karst region.Progress in Geophysics (in Chinese),24(1):340-347.

    Li Z X,Tan H D,Fu S S,et al.2015.Two-dimensional synchronous inversion of TDIP with cross-gradient constraint.Chinese J.Geophys.(in Chinese),58(12):4718-4726,doi:10.6038/cjg20151232.

    Liu B,Li S C,Li S C,et al.2012.3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization.Chinese J.Geophys.(in Chinese),55(1):260-268,doi:10.6038/j.issn.0001-5733.2012.01.025.

    Liu Y H,Yin C C.2013.3D inversion for frequency-domain HEM data.Chinese J.Geophys.(in Chinese),56(12):4278-4287,doi:10.6038/cjg20131230.

    Loke M H,Barker R D.1996.Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method.Geophysical Prospecting,44(1):131-152.

    Loke M H,Dahlin T.2002.A comparison of the Gauss-Newton and quasi-Newton methods in resistivity imaging inversion.Journal of Applied Geophysics,49(3):149-162.

    Maceira M,Ammon C J.2009.Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure.Journal of Geophysical Research:Solid Earth (1978—2012),114(B2),doi:10.1029/2007JB005157.

    McGillivray P R,Oldenburg D W.1990.Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem:a comparative study.Geophysical Prospecting,38(5):499-524.

    Musil M,Maurer H R,Green A G.2003.Discrete tomography and joint inversion for loosely connected or unconnected physical properties:Application to crosshole seismic and georadar data sets.Geophysical Journal International,153(2):389-402.

    Osher S,Fedkiw R.2003.Level Set Methods and Dynamic Implicit Surfaces.New York:Springer.

    Paige C C,Saunders M A.1982.LSQR:An algorithm for sparse linear equations and sparse least squares.ACM Transactions on Mathematical Software (TOMS),8(1):43-71.

    Pan J S,Ge W Z,Zhe J P.2010.High-density electrical resistivity imaging technique of ground/borehole-to-surface/inter-well.Journal of North China Institute of Water Conservancy and Hydroelectric Power (in Chinese),31(2):74-78.

    Pan K J,Tang J T,Hu H L,et al.2012.Extrapolation cascadic multigrid method for 2.5D direct current resistivity modeling.Chinese J.Geophys.(in Chinese),55(8):2769-2778,doi:10.6038/j.issn.0001-5733.2012.08.028.

    Peng M.2012.Joint inversion of magnetotelluric and teleseismic data [Ph.D.thesis] (in Chinese).Beijing:China University of Geosciences (Beijing).

    Pidlisecky A,Haber E,Knight R.2007.RESINVM3D:A 3D resistivity inversion package.Geophysics,72(2):H1-H10.

    Popovici A M,Sethian J A.2002.3-D imaging using higher order fast marching traveltimes.Geophysics,67(2):604-609.

    Rawlinson N,Sambridge M.2003.Seismic traveltime tomography of the crust and lithosphere.Advances in Geophysics,46:81-199.

    Ritter J R R,Jordan M,Christensen U R,et al.2001.A mantle plume below the Eifel volcanic fields,Germany.Earth and Planetary Science Letters,186(1):7-14.

    Sambridge M S.1990.Non-linear arrival time inversion:constraining velocity anomalies by seeking smooth models in 3-D.Geophysical Journal International,102(3):653-677.

    Tikhonov A N,Arsenin V Y.1977.Solutions of Ill-Posed Problems.Washington,DC:Vh Winston.Treitel S,Lines L.1988.Geophysical examples of inversion (with a grain of salt).The Leading Edge,7(11):32-35.

    Vidale J.1988.Finite-difference calculation of travel times.Bulletin of the Seismological Society of America,78(6):2062-2076.

    Vidale J E.1990.Finite-difference calculation of traveltimes in three dimensions.Geophysics,55(5):521-526.

    Wang J.2014.Research of the theory of joint inversion of gravity and magnetic data based on cross-gradient function [Master′s thesis] (in Chinese).Beijing:China University of Geosciences (Beijing).

    Wu X P,Xu G M.2000.Study on 3-D resistivity inversion using conjugate gradient method.Chinese J.Geophys.(in Chinese),43(3):420-427.

    Wu X P,Liu Y,Wang W.2015.3D resistivity inversion incorporating topography based on unstructured meshes.Chinese J.Geophys.(in Chinese),58(8):2706-2717,doi:10.6038/cjg20150808.

    Yu P,Wang J L,Wu J S,et al.2006.Review and discussions on geophysical joint inversion.Progress in Exploration Geophysics (in Chinese),29(2):87-93.

    Zelt C A,Smith R B.1992.Seismic travel time inversion for 2-D crustal velocity structure.Geophysical Journal International,108(1):16-34.

    Zhang F X,Wu Q J,Li Y H,et a1.2010.Application of FMM ray tracing to forward and inverse problems of seismology.Progress in Geophysics (in Chinese),25(4):1197-1205,doi:10.3969/j.issn.1004-2903.2010.04.008.

    Zhou L F.2012.Two dimensional joint inversion of MT and seismic data [Master′s thesis] (in Chinese).Beijing:China University of Geosciences (Beijing).

    附中文參考文獻(xiàn)

    范翠松,李桐林,嚴(yán)加永.2012.2.5維復(fù)電阻率反演及其應(yīng)用試驗(yàn).地球物理學(xué)報(bào),55(12):4044-4050,doi:10.6038/j.issn.0001-5733.2012.12.016.

    馮銳,李智明,李志武等.2004.電阻率層析成像技術(shù).中國地震,20(1):13-30.

    胡祥云,楊迪坤,劉少華等.2006.環(huán)境與工程地球物理的發(fā)展趨勢.地球物理學(xué)進(jìn)展,21(2):598-604.

    雷旭友,李正文,折京平.2009.超高密度電阻率法在土洞、煤窯采空區(qū)和巖溶勘探中應(yīng)用研究.地球物理學(xué)進(jìn)展,24(1):340-347.李兆祥,譚捍東,付少帥等.2015.基于交叉梯度約束的時(shí)間域激發(fā)極化法二維同步反演.地球物理學(xué)報(bào),58(12):4718-4726,doi:10.6038/cjg20151232.

    劉斌,李術(shù)才,李樹忱等.2012.基于不等式約束的最小二乘法三維電阻率反演及其算法優(yōu)化.地球物理學(xué)報(bào),55(1):260-268,doi:10.6038/j.issn.0001-5733.2012.01.025.

    劉云鶴,殷長春.2013.三維頻率域航空電磁反演研究.地球物理學(xué)報(bào),56(12):4278-4287,doi:10.6038/cjg20131230.

    潘紀(jì)順,葛為中,折京平.2010.地面/井地/井間超高密度電阻率成像技術(shù).華北水利水電學(xué)院學(xué)報(bào),31(2):74-78.

    潘克家,湯井田,胡宏伶等.2012.直流電阻率法2.5維正演的外推瀑布式多重網(wǎng)格法.地球物理學(xué)報(bào),55(8):2769-2778,doi:10.6038/j.issn.0001-5733.2012.08.028.

    彭淼.2012.大地電磁與天然地震數(shù)據(jù)聯(lián)合反演研究[博士論文].北京:中國地質(zhì)大學(xué)(北京).

    王俊.2014.基于交叉梯度重磁聯(lián)合反演方法技術(shù)研究[碩士論文].北京:中國地質(zhì)大學(xué) (北京).

    吳小平,徐果明.2000.利用共軛梯度法的電阻率三維反演研究.地球物理學(xué)報(bào),43(3):420-427.

    吳小平,劉洋,王威.2015.基于非結(jié)構(gòu)網(wǎng)格的電阻率三維帶地形反演.地球物理學(xué)報(bào),58(8):2706-2717,doi:10.6038/cjg20150808.

    于鵬,王家林,吳健生等.2006.地球物理聯(lián)合反演的研究現(xiàn)狀和分析.勘探地球物理進(jìn)展,29(2):87-93.

    張風(fēng)雪,吳慶舉,李永華等.2010.FMM射線追蹤方法在地震學(xué)正演和反演中的應(yīng)用.地球物理學(xué)進(jìn)展,25(4):1197-1205,doi:10.3969/j.issn.1004-2903.2010.04.008.

    周麗芬.2012.大地電磁與地震數(shù)據(jù)二維聯(lián)合反演研究[碩士論文].北京:中國地質(zhì)大學(xué) (北京).

    (本文編輯 何燕)

    林方麗,王光杰,楊曉勇.2016.綜合電磁法在礦區(qū)深部成礦機(jī)制中的應(yīng)用研究——以皖南烏溪多金屬礦區(qū)為例.地球物理學(xué)報(bào),59(11):4323-4337,doi:10.6038/cjg20161132.

    Lin F L,Wang G J,Yang X Y.2016.Application of comprehensive electromagnetic study in deep mineralization mechanism—A case study of the Wuxi polymetallic ore deposit,south Anhui.Chinese J.Geophys.(in Chinese),59(11):4323-4337,doi:10.6038/cjg20161132.

    綜合電磁法在礦區(qū)深部成礦機(jī)制中的應(yīng)用研究——以皖南烏溪多金屬礦區(qū)為例

    林方麗1,2,王光杰1,楊曉勇3

    1 中國科學(xué)院地質(zhì)與地球物理研究所,北京 1000292 中國科學(xué)院大學(xué),北京 1000493 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院,合肥 230026

    摘要 烏溪礦區(qū)介于長江中下游多金屬成礦帶與華南成礦帶之間,是江南造山帶上的一個(gè)找礦新區(qū).本文在這一找礦新區(qū)開展了有效的電磁方法綜合勘探,試圖對(duì)該區(qū)深部成礦機(jī)制進(jìn)行研究.我們收集了該區(qū)的地質(zhì)地球化學(xué)資料,從地質(zhì)地球化學(xué)資料中分析了區(qū)域成礦背景;我們還采集了區(qū)域出露的主要巖石類型,在實(shí)驗(yàn)室開展物性測定,在此基礎(chǔ)上選擇了三種電磁方法開展研究區(qū)的野外測量.通過地面高精度磁測、激發(fā)極化法和可控源音頻大地電磁法(CSAMT)三種地球物理方法,開展了研究區(qū)的磁化率、極化率和電阻率的分布特征研究.深入分析了地質(zhì)、地球化學(xué)和地球物理三種資料與成礦的關(guān)系,相互約束,探討成礦模型、成礦機(jī)制和成礦的可能性,推測可能的礦體賦存位置和深度范圍.在地球物理研究結(jié)果基礎(chǔ)上,結(jié)合地質(zhì)和地球化學(xué)背景資料,構(gòu)建了研究區(qū)的成礦動(dòng)力學(xué)模型,推測了研究區(qū)成礦機(jī)制,揭示了礦區(qū)的成礦潛力.綜合所得結(jié)果布設(shè)了鉆孔,鉆探結(jié)果揭示了研究區(qū)深部存在強(qiáng)烈礦化蝕變和強(qiáng)蝕變斑巖,初步確認(rèn)為斑巖型礦床,與CSAMT剖面的解釋結(jié)果基本一致,也與推測的測區(qū)內(nèi)深部成礦機(jī)制相吻合.鉆孔結(jié)果和地球物理結(jié)果的一致性以及對(duì)已構(gòu)建的成礦動(dòng)力學(xué)模型的支持,充分證明了綜合電磁方法在斑巖型多金屬礦床的發(fā)現(xiàn)和預(yù)測中的重要作用,同時(shí)斑巖型礦床的確認(rèn)進(jìn)一步印證了華南成礦帶與俯沖作用形成的成礦帶的相似性,從而推動(dòng)整個(gè)華南地區(qū)的成礦地質(zhì)研究.

    關(guān)鍵詞 成礦機(jī)制;綜合電磁法;烏溪多金屬礦;斑巖

    作者簡介 林方麗,女,1985年生,中國科學(xué)院地質(zhì)與地球物理研究所在讀博士研究生,主要從事電磁法理論及應(yīng)用研究.

    E-mail:lgg1287@163.com

    *通訊作者 王光杰,男,1966年生,副研究員,主要從事電磁法勘探研究和應(yīng)用工作.E-mail:gjwang@mail.iggcas.ac.cn

    Abstract The Wuxi polymetallic ore deposit which lies between the middle and lower reaches of the Yangtze River polymetallic metallogenic belt and South China metallogenic belt,is a new prospecting area in the Jiangnan orogen.This article is trying to study the mineralization mechanism in the deep subsurface of this area with the effective electromagnetic integrated exploration.For a new prospecting area,the collection of regional geological and geophysical background is very important.Geological data in this area is collected,which permits to analyze the regional metallogenic background.The physical property of the main rock types exposed in this area is measured in the laboratory.On the basis of these results three electromagnetic field measurements are carried out.By the measurements of three kinds of geophysical methods:ground precise magnetic survey,induced polarization and controlled source audio magnetotelluric method (CSAMT),the distribution characteristics of magnetic susceptibility,polarization and resistivity are studied,to bound mineralization possibilities and possible ore body position and depth range together.Based on the results of the geophysical studies,combined with geological and geochemical data,the mineralization dynamics model of the study area is built and then the mineralization mechanism is suggested.Drillings are deployed according to all these results.Drilling results prove the presence of strong alteration,mineralization and porphyry alteration at depth in the study area,which is consistent with the judgment of the CSAMT results.The consistency of drilling and geophysical results fully demonstrate the important role of the integrated electromagnetic method in the discovery and prediction of porphyry type polymetallic deposits.Meanwhile,it confirms the similarity between South China metallogenic belt and the mineralization belt formed by the subduction,thus promoting mineralization geological studies throughout South China.

    Keywords Metallogenic Mechanism;Integrated Approach of Electric Magnetic Method;Wuxi Polymetallic Ore;Porphyry

    1 引言

    斑巖型礦床以其巨大的經(jīng)濟(jì)價(jià)值和重要的學(xué)術(shù)意義一向受到人們的重視.近年研究表明很多大型、超大型斑巖型銅、金礦都與洋脊俯沖作用所產(chǎn)生的中酸性巖漿巖(多為斑巖)密切相關(guān)(Cooker et al.,2005),因此斑巖型礦床的成礦機(jī)制與板塊俯沖的關(guān)系得到廣泛的研究,并且取得了豐富的研究成果(Sun et al.,2003,2007,2013;Liu et al.,2010;Wang et al.,2013).孫衛(wèi)東等認(rèn)為在中國東部地區(qū)尋找洋脊俯沖的跡象,將有助于國內(nèi)大型、超大型斑巖礦床以及其他相關(guān)礦床的發(fā)現(xiàn)(Ling et al.,2009;Sun et al.,2010).本次工作就將目標(biāo)區(qū)域選在介于長江中下游多金屬成礦帶與華南成礦帶之間的烏溪礦區(qū),這里是江南造山帶上尋找斑巖型礦床的理想場所.烏溪礦區(qū)位于安徽省涇縣榔橋鎮(zhèn),基于目前勘探開發(fā)的地質(zhì)資料,應(yīng)該屬于熱液脈型成礦類型.根據(jù)區(qū)域地質(zhì)調(diào)查工作,發(fā)現(xiàn)了地表存在多個(gè)金礦化點(diǎn),并圈定了成礦有利區(qū)域.前人在此基礎(chǔ)上對(duì)該區(qū)金礦成礦條件和成礦后期構(gòu)造進(jìn)行了分析(劉惠華,2003;劉惠華和朱寧,2004),總結(jié)了該區(qū)的地質(zhì)特征及控礦因素(趙永利等,2013).

    雖然斑巖型礦床研究已經(jīng)在許多領(lǐng)域取得了較大的進(jìn)展,但許多關(guān)鍵性的問題現(xiàn)在仍沒有得到妥善解決.最使人感到棘手的問題就是在同一區(qū)域內(nèi)幾乎是同年齡同源同成份的許多巖體,含礦性往往差異極大.如何區(qū)分在成礦有利區(qū)域的圍巖和礦體成為斑巖型礦床研究中亟待解決的地質(zhì)問題.而巖石的蝕變和礦化蝕變,會(huì)引起巖石的電性、磁性等物理性質(zhì)的改變,因此在金屬礦的成礦有利區(qū)域開展有效的綜合電磁法工作,能夠從巖石的電磁性質(zhì)的差異中分辨出深部可能的含礦構(gòu)造(鄧居智等,2015;張昆等,2015).

    由于不同的地球物理方法只能反映巖礦石一個(gè)物理參數(shù)的差異,因此在實(shí)際應(yīng)用中通常同時(shí)使用多種不同的方法(郝天珧和江為為,1998;孫興國等,2007;陳偉軍等,2008;崔敏利等,2010;呂慶田等,2015;邵陸森等,2015;萬漢平等,2015;徐興倩等,2015),共同約束地下礦物的性質(zhì),盡可能的減少多解性.近年來,地球物理綜合找礦方法已在金屬礦勘探中發(fā)揮了主要作用(Bastani et al.,2009;張壯等,2013;萬芬等,2014;丁高明等,2015;張光之等,2015;王顯瑩等,2015).磁法勘探根據(jù)不同巖礦石的磁性差異,一般用于控制不同巖性的邊界及深大斷裂(婁德波等,2008),圈定在金屬礦成礦過程中形成的伴生礦物的磁異常,如磁鐵礦和磁黃鐵礦.在金屬礦的形成過程中,成礦物質(zhì)的匯聚也需要相應(yīng)的成礦通道和成礦空間,因此金屬礦勘查的關(guān)鍵就是尋找相應(yīng)的成礦通道和成礦空間,所以磁法勘探也被廣泛地應(yīng)用于金屬礦的勘查(吳國學(xué),2007;盧焱等,2008).激發(fā)極化法是借助多金屬硫化物礦物普遍存在著激發(fā)極化特性(何繼善,2006),大量的實(shí)際工作已經(jīng)證明了激電法的找礦效果(武煒等,2009;楊振威等,2013;Smith,2014;柳建新等,2004),激電測量也將有利于成礦機(jī)制的判斷.可控源音頻大地電磁法(CSAMT)則是根據(jù)電磁場的趨膚效應(yīng),通過不同的頻率的電磁信號(hào)的差異,計(jì)算和反演不同深度的電阻率分布差異(王若等,2014;孫博等,2015;湯井田等,2015;周聰?shù)龋?015).已在各種地質(zhì)勘查(陳玉玲等,2015)和深部探測(于昌明,1998;李帝銓等,2008;時(shí)彬,2012;譚章坤,2013)中成功應(yīng)用.綜上所述可知,在多金屬礦床的勘查中,地面高精度磁測、激發(fā)極化法和CSAMT電阻測深三種方法的綜合使用,對(duì)于成礦動(dòng)力學(xué)模型的建立,可以得到多參數(shù)多研究角度的立體模型.

    雖然烏溪礦區(qū)研究程度較低,作為研究中國斑巖型礦床與板塊俯沖關(guān)系的最有利區(qū)域,前人也進(jìn)行了大量的地球化學(xué)的研究(李雙等,2012,2014,2015).地球化學(xué)結(jié)果表明,烏溪礦區(qū)不論從成礦年齡還是成礦元素分析,都與德興斑巖型礦床有較好的一致性,進(jìn)一步支持了斑巖型礦床的推測,但具體工作中一直沒有找到深部斑巖證據(jù).為了更深入的分析烏溪礦區(qū)的礦床類型和深部成礦,我們?cè)诘刭|(zhì)和地球化學(xué)工作的基礎(chǔ)上,采集了區(qū)內(nèi)的巖石標(biāo)本進(jìn)行了巖石物性結(jié)構(gòu)特征分析,據(jù)此選取了上述三種方法,進(jìn)行了一系列綜合電磁法勘探.本文依據(jù)當(dāng)?shù)氐膸r石物性結(jié)構(gòu)特征,根據(jù)這三種不同的電磁方法的結(jié)果,結(jié)合已有的地質(zhì)地球化學(xué)資料構(gòu)建了地質(zhì)地球物理成礦動(dòng)力學(xué)模型,并在已有研究基礎(chǔ)上初步給出了探測區(qū)的深部成礦機(jī)制.成礦機(jī)制的研究精確定位了礦區(qū)可能的成礦位置,同時(shí)進(jìn)行了多個(gè)深鉆驗(yàn)證,鉆井巖芯結(jié)果表明,礦區(qū)深部出現(xiàn)多個(gè)蝕變斑巖脈和普遍的金屬礦化,可以確定深部存在較大的斑巖礦體,具有很好的成礦前景.該結(jié)論也支持了我國東南部可能存在的與板塊俯沖有關(guān)的斑巖型礦床,在相同的成礦有利區(qū)域內(nèi)可以進(jìn)一步尋找可能的斑巖型礦床.

    2 研究區(qū)地質(zhì)地球物理背景

    2.1 區(qū)域地質(zhì)概況

    烏溪礦區(qū)位于揚(yáng)子地臺(tái)江南古陸北側(cè),江南大斷裂與東西向周王斷裂交匯部位的南側(cè),北北東向湯口斷裂束分支在區(qū)內(nèi)通過(圖1).區(qū)域上,從震旦系至三疊系的地層均有出露.研究區(qū)內(nèi)出露地層主要為志留系和泥盆系.研究區(qū)東南部出露的榔橋巖體,是皖南地區(qū)燕山期一個(gè)重要的花崗閃長巖基,成巖時(shí)代被限定為137~139 Ma之間(李雙等,2014).

    2.2 礦區(qū)地質(zhì)背景

    礦區(qū)內(nèi)主要發(fā)育近北北東向區(qū)域斷裂構(gòu)造-即湯口斷裂(圖1),為構(gòu)造熱液型成礦的主要的控礦、含礦構(gòu)造.區(qū)內(nèi)局部可見花崗斑巖脈出露,走向以北東向?yàn)橹鳎瑢?5~25 m,長100~300 m.研究區(qū)內(nèi)地表普遍發(fā)育黃鐵礦化等強(qiáng)烈蝕變,并局部具有低品位的金礦化.靠近巖脈或斷裂,巖石較破碎,硅化、絹云母化、黃鐵礦化、黃鐵絹英巖化等蝕變發(fā)育.其它方向的構(gòu)造,含礦性較差,大多為成礦后期破壞性構(gòu)造,對(duì)早期形成的礦化有一定程度的破壞作用.礦區(qū)內(nèi)出露地層主要為志留系粉砂巖、泥質(zhì)粉砂巖和泥盆系石英細(xì)砂巖.地層產(chǎn)狀較平緩,傾角10°~25°,節(jié)理較發(fā)育.

    2.3 區(qū)域地球物理背景

    航磁異常結(jié)果顯示研究區(qū)屬于皖南強(qiáng)磁場區(qū),區(qū)域以一系列較緊密排列的NE向、NNE向條帶狀磁異常帶分布為主(王建偉等,2009).研究區(qū)處于正負(fù)磁異常的交界位置,推測與區(qū)內(nèi)斷裂分布相關(guān).巖體出露區(qū)域顯示團(tuán)塊狀正磁異常,說明區(qū)域巖體磁化率相對(duì)較高.

    安徽省布格重力異常結(jié)果與大地構(gòu)造分區(qū)的地質(zhì)礦產(chǎn)要素對(duì)應(yīng)分析,可以將三個(gè)重力異常級(jí)別對(duì)應(yīng)于不同的構(gòu)造分區(qū)(蘭學(xué)毅等,2012).研究區(qū)位于揚(yáng)子陸塊—江南地塊—皖南褶皺帶中,區(qū)域重力異常低值中心指向皖南花崗巖集中區(qū).烏溪礦區(qū)位于布格重力異常高低交變區(qū)域,因此在此區(qū)域內(nèi)重力處于不均衡狀態(tài),易產(chǎn)生地殼的均衡運(yùn)動(dòng),而地殼的均衡運(yùn)動(dòng)會(huì)進(jìn)一步帶動(dòng)區(qū)內(nèi)構(gòu)造薄弱地帶的巖石破碎,從而有利于深部物質(zhì)的帶入.

    2.4 測區(qū)巖石地球物理特征

    根據(jù)區(qū)內(nèi)巖石的出露和鉆探得到的巖芯情況,我們分別選取了泥巖、砂巖、花崗巖及含硫化物的巖石進(jìn)行磁化率、極化率及電阻率的測量.從測量結(jié)果來看,蝕變泥巖由于含水較多,電阻率普遍較小.綜合來看存在礦化蝕變的各類巖石電阻率相對(duì)沒有蝕變的巖石電阻率低,極化率相對(duì)偏高,花崗巖類是礦區(qū)中唯一磁化率較高的巖石類型,而圍巖如砂巖、泥巖等,不論是極化率還是磁化率都很小,說明區(qū)域圍巖不具有磁性也不具有激電特性.

    圖1 研究區(qū)區(qū)域構(gòu)造簡圖(安徽省地質(zhì)礦產(chǎn)局,1987)1.元古-古生代地層;2.晚古生代-中生代;3.中生代-第三紀(jì);4.花崗巖;5.花崗閃長巖;6.二長花崗巖;7.斷裂;8.測區(qū).Fig.1 Regional tectonic map (Anhui Bureau of Geology and Mineral Resources,1987)1.Proterozoic-Paleozoic strata;2.Late Paleozoic-Mesozoic;3.Mesozoic-Tertiary;4.Granite;5.Granodiorite;6.Adamellite;7.Fault;8.Study area.

    圖2 烏溪測區(qū)地質(zhì)與電磁測點(diǎn)分布圖1.第四紀(jì);2.石炭系;3.泥盆系;4.志留系;5.花崗閃長巖;6.花崗斑巖;7.見礦鉆孔;8.金礦化點(diǎn);9.CSAMT測點(diǎn);10.磁法測點(diǎn);11.IP測點(diǎn).Fig.2 Geology map with electrical and magnetic prospecting site distribution in Wuxi survey area1.Quaternary;2.Carboniferous;3.Devonian;4.Silurian;5.Granodiorite;6.Porphyry;7.Borehole;8.Gold ore locality;9.CSAMT Site;10.Magnetic Site;11.IP Site.

    圖3 烏溪測區(qū)磁異常分布平面圖1.道路;2.磁場不均勻區(qū)域.Fig.3 Magnetic anomaly distribution of Wuxi survey area1.Road;2.Nonuniform magnetic field area.

    圖4 烏溪測區(qū)極化率異常平面圖1.道路;2.金礦化點(diǎn);3.見礦鉆孔;4.磁測區(qū)域Fig.4 Polarization anomaly distribution of Wuxi survey area1.Road;2.Gold mineralization locality;3.Borehole;4.Magnetic survey area.

    圖5 雙頻激電剖面極化率和電阻率與CSAMT反演電阻率斷面對(duì)比(5000線)Fig.5 Comparisoll of resistivity from dual frequency IP method and CSAMT(line 5000)

    圖6 電阻率水平切片圖(海拔高度-600 m~+200 m)Fig.6 Horizontal resistance slices at elevation from +200 m to -600 m

    圖7 烏溪電阻率立體圖Fig.7 A perspective view of resistivity in Wuxi

    表1 烏溪測區(qū)巖石物理特征(萬芬,2014)

    3 電磁綜合勘探工作

    根據(jù)地質(zhì)資料,礦區(qū)內(nèi)斷裂構(gòu)造比較發(fā)育,并且在斷裂附近形成巖脈,成礦較為有利,因此將本次工作的重點(diǎn)圍繞淺表已知的斷裂和巖脈展開.為了了解斷裂在地下的走向及傾角,選擇垂直主斷裂的北西45°方向作為測線方向(圖2),為了有效的控制異常,將測網(wǎng)布設(shè)為200×25 m.在工作重點(diǎn)區(qū)域,同時(shí)進(jìn)行地面高精度磁法、雙頻激電法和CSAMT三種電磁方法,這樣不僅能夠控制區(qū)內(nèi)磁性物質(zhì)分布,還可以進(jìn)一步了解不同電阻率的巖石在地下的分布情況,從而尋找成礦有利區(qū)域.實(shí)際工作中為了進(jìn)一步控制磁性異常和激電異常的邊界范圍,將激電和磁法測線都進(jìn)行了外延和加密.在地面高精度磁測過程中,針對(duì)測區(qū)中磁場變化劇烈的區(qū)域進(jìn)行了測點(diǎn)加密,進(jìn)一步確認(rèn)異常的真實(shí)可靠性.

    3.1 地面高精度磁測結(jié)果及解釋

    圖3是烏溪測區(qū)地面磁測異常圖.磁異常結(jié)果顯示,在測區(qū)的北側(cè)主要分布團(tuán)塊狀正異常,異常范圍寬度500~1000 m,從北向南磁異常逐漸減弱,在測區(qū)東南側(cè)顯示有北東向分布的負(fù)磁異常.引起磁異常減弱可能的原因有兩種,磁性物質(zhì)減少或者磁性物質(zhì)分布深度增加.而在測區(qū)東南部相對(duì)負(fù)磁異常處,相應(yīng)地質(zhì)圖上顯示為第四紀(jì)覆蓋區(qū),由于第四紀(jì)地層不具有磁性,由此推測引起區(qū)內(nèi)磁異常減弱的原因可能是磁性物質(zhì)的埋深相對(duì)增加,同時(shí)不排除是在磁場變化區(qū)域存在斷裂,使得兩側(cè)巖石磁性差異.而區(qū)內(nèi)標(biāo)本測量結(jié)果顯示(表1),花崗巖類是唯一磁化率較高的巖石類型,從而初步推測區(qū)內(nèi)花崗巖巖體在礦區(qū)的分布范圍.測區(qū)中間分布大范圍的磁場不均勻區(qū)域,根據(jù)巖石標(biāo)本測量結(jié)果,蝕變和礦化都可以引起磁化率的減弱,因此推測該區(qū)域的不均勻磁場可能與巖石的蝕變相關(guān).

    3.2 激電觀測結(jié)果及解釋

    圖4是烏溪測區(qū)極化率異常平面圖.極化率結(jié)果顯示,高極化異常集中分布在測區(qū)中間部位,區(qū)內(nèi)存在一條和測線正交的高極化帶貫穿整個(gè)測區(qū),寬度500 m以上.根據(jù)前期地質(zhì)勘探和鉆探結(jié)果,區(qū)內(nèi)已圈定成礦帶2個(gè)和多條礦脈,與地表出露的花崗斑巖脈伴生.根據(jù)鉆探結(jié)果(劉琛琛和楊錢江,2012),花崗斑巖脈傾向南東,在巖脈之下為志留系唐家塢組粉砂巖,在唐家塢組內(nèi)陸續(xù)分布多處礦化破碎蝕變,初步判斷為地表礦脈在深部的延伸.已查明異常均位于高極化帶內(nèi),因此可以推測區(qū)內(nèi)高極化帶區(qū)域都存在很好的成礦條件.另外測區(qū)內(nèi)已查明多個(gè)金礦化點(diǎn),都位于高極化區(qū)域及邊界位置,推測是含礦熱液被圍巖阻擋,在斑巖-圍巖接觸帶產(chǎn)出一系列礦脈和金礦化點(diǎn).因此根據(jù)地表出露的斑巖脈、礦化蝕變帶、金礦化點(diǎn)和地表高極化異常帶,可以推測在高極化率異常之下,隱伏有一個(gè)含礦斑巖體,預(yù)測礦區(qū)具有很好的成礦遠(yuǎn)景.

    3.3 CSAMT觀測結(jié)果及解釋

    將可控源音頻大地電磁法的結(jié)果經(jīng)過近場校正、靜態(tài)校正和曲線平滑處理之后,再進(jìn)行一維電阻率深度反演計(jì)算,即可得到每條測線的電阻率隨深度變化的剖面圖.對(duì)比電阻率測深結(jié)果和激電測量的地表電阻率,兩者基本保持一致,進(jìn)一步驗(yàn)證兩種方法的可靠性(圖5).從電阻率測深結(jié)果來看,電阻率分層的傾向與測線方向相同,傾角30~40°,進(jìn)一步印證了前期的地質(zhì)工作得到的地層傾向南東的結(jié)果.從測深剖面結(jié)果來看,在測線800~1000 m之間存在縱向延伸較深的低阻異常帶,可推測為斷裂破碎帶,同時(shí)在破碎帶的邊界,對(duì)應(yīng)了兩個(gè)高極化率異常,電阻率低值與高極化率對(duì)應(yīng),較好的對(duì)應(yīng)了含礦巖石的低阻高極化特性,因此可以推測區(qū)內(nèi)蝕變礦化與斷裂的相關(guān)關(guān)系.

    根據(jù)測區(qū)的CSAMT電阻率測深結(jié)果進(jìn)行切片處理(圖6),電阻率相對(duì)較高的點(diǎn)在地表分布較分散,而隨著深度的增加逐漸向小號(hào)測線和測線小號(hào)點(diǎn)匯聚,整體電阻率呈條帶狀分布.從地表的電阻率水平切片圖上(H=200 m)看,高阻體的延伸方向由南北向逐漸轉(zhuǎn)為北東向,與測線方向垂直.整體切片圖的結(jié)果表明隨著深度的增加,測區(qū)深部整體電阻率較地表電阻率偏低,且又呈帶狀分布,推測是在構(gòu)造作用下圍巖受到強(qiáng)烈蝕變改造而形成,是可能的成礦通道和成礦空間.測區(qū)北半部分的高阻礦脈逐漸向測區(qū)邊界匯聚,測區(qū)南半部分高阻脈也向西南方向匯聚.而在測區(qū)中間位置形成大規(guī)模的低阻異常帶,深部如此大規(guī)模的低阻異??梢酝茰y深部存在較大的斷裂構(gòu)造.這也與區(qū)內(nèi)測量得到的大規(guī)模極化率異常相對(duì)應(yīng),從而推測區(qū)內(nèi)斷裂構(gòu)造為深部巖漿活動(dòng)提供了大范圍的活動(dòng)空間.

    從地質(zhì)圖上看,測區(qū)內(nèi)出露的花崗斑巖脈的走向基本與CSAMT測量結(jié)果的高阻延伸方向一致.測區(qū)中間出現(xiàn)的北東東向斷層有部分花崗斑巖出露,電阻率切片圖上顯示出高阻特性.而測區(qū)東南側(cè)出露的含礦化的花崗斑巖脈,電阻率測量結(jié)果相對(duì)較小,前期鉆探結(jié)果顯示該點(diǎn)的花崗斑巖脈出露地表,厚度只有50 m,并且傾向東南側(cè).根據(jù)測區(qū)電阻率立體圖(圖7)可以看到,電阻率分布情況也傾向南東,與測區(qū)內(nèi)地層的傾向相一致,說明巖體存在順層侵入.

    3.4 電磁法綜合解釋與成礦動(dòng)力來源

    通過地面高精度磁測結(jié)果,磁場有從南到北逐漸增大的趨勢,結(jié)合礦區(qū)地質(zhì)資料,我們可以推測得出巖漿來源于礦區(qū)東南側(cè)的深部巖基,沿著區(qū)內(nèi)巖層、斷裂、裂隙逐漸上侵到礦區(qū)西北側(cè)的淺部.而極化率和電阻率的結(jié)果都顯示礦區(qū)存在一條北東向的高極化率和低電阻異常帶,并且兩者在空間上基本可以完全吻合,因此可以基本得出區(qū)內(nèi)存在的低阻高極化區(qū)域?yàn)閰^(qū)內(nèi)最有利的成礦空間.結(jié)合三種電磁方法結(jié)果,依據(jù)巖石樣品的測量結(jié)果,測區(qū)西北角表現(xiàn)為高磁低極化特征,說明該區(qū)域內(nèi)巖石可能為花崗斑巖或者花崗巖巖體,同時(shí)電阻率表現(xiàn)為高阻特征,進(jìn)一步證明了此處可能的巖體存在.考慮到區(qū)域地層傾向,此處巖體可能是南東側(cè)深部巖體的延伸.測區(qū)中間測得的高極化帶與低阻構(gòu)造帶相對(duì)應(yīng),隨著深度的增加電阻率逐漸減小,推測表層巖石沒有受到蝕變作用的影響,隨著深度的增加深部巖石逐漸發(fā)生蝕變作用逐漸強(qiáng)烈.同時(shí)低阻條帶在深部范圍更大,說明了高極化帶不僅在淺部有較好的成礦可能性,在深部有更大的礦化蝕變空間.另外高極化帶區(qū)域磁場分布的不均勻性說明成礦過程中可能有磁性物質(zhì)的生成.電阻率測深結(jié)果顯示,高阻體在深部匯聚于測區(qū)西南角,而對(duì)應(yīng)的磁場和極化率情況顯示該區(qū)磁場較高極化率很低,因此推斷測區(qū)西南角部分深部可能沒有受到熱液活動(dòng)的影響,沒有成礦條件.

    綜合所得到的地球物理和前期的地質(zhì)、地球化學(xué)結(jié)果,將分別從以下三個(gè)不同的尺度分析礦區(qū)的成礦動(dòng)力來源:

    (1) 從烏溪礦區(qū)的綜合電磁法結(jié)果分析,可以初步推斷礦區(qū)的成礦動(dòng)力模型:磁異常表明巖漿可能來源于礦區(qū)東南部,而礦區(qū)中部的磁場不均勻區(qū)域推測為后期熱液蝕變?cè)斐傻?;極化率結(jié)果顯示礦區(qū)中間存在大面積的高極化率異常帶,推測為巖體和熱液活動(dòng)形成的礦化蝕變沿?cái)嗔研纬傻牡V化帶;電阻率測深結(jié)果進(jìn)一步描述了礦區(qū)深部存在的低電阻率異常帶,同時(shí)與高極化異常帶位置基本一致.因此推測巖漿熱液來源于深部,以區(qū)內(nèi)存在的大量斷裂裂隙為活動(dòng)空間上侵進(jìn)入淺層.單純的巖漿活動(dòng)不足以解釋磁異常與極化率/電阻率異常的分布特征,因此推測區(qū)內(nèi)存在巖漿期后熱液活動(dòng).熱液活動(dòng)范圍主要在極化率異常范圍附近,并沒有到達(dá)之前的巖漿活動(dòng)的邊界,因此區(qū)內(nèi)可能存在未蝕變巖體,磁異常結(jié)果可以驗(yàn)證這一推測.同時(shí)巖心鉆探結(jié)果表明在高極化區(qū)域深部在巖漿普遍存在隱爆角礫巖和巖石物性測試結(jié)果都能證明區(qū)內(nèi)存在熱液活動(dòng).熱液聚集過程中在深部形成壓力和溫度的增大,從而為引爆作用的發(fā)生提供了動(dòng)力來源,熱液上侵過程中對(duì)圍巖進(jìn)行蝕變改造和礦物富集,在成礦有利空間金屬礦物沉積成礦.

    (2) 對(duì)于礦區(qū)內(nèi)巖漿熱液來源及其動(dòng)力學(xué)背景,可以從區(qū)域地質(zhì)方面進(jìn)一步討論.對(duì)礦區(qū)周圍巖體的成礦專屬性分析、巖石學(xué)分析和地球化學(xué)特征分析表明,榔橋巖體為區(qū)域內(nèi)主要賦礦圍巖.對(duì)礦區(qū)內(nèi)出露花崗斑巖以及鉆孔樣品年代學(xué)測定結(jié)果表明,榔橋巖體成巖時(shí)代與烏溪礦區(qū)巖體的形成時(shí)代基本一致,分別為139.6±1.7 Ma (ZK7301)、137.3±1.6 Ma (ZK7001)、137.3±1.1 Ma (10WX-1),推測榔橋巖體與礦區(qū)斑巖脈屬于同期巖漿作用(李雙等,2015),這期巖漿活動(dòng)為皖南地區(qū)普遍的多金屬礦化期.區(qū)域內(nèi)分布多條深大斷裂,其中北北東向湯口斷裂束切割整個(gè)榔橋巖體并從礦區(qū)穿過,為巖漿和成礦流體提供了充分的運(yùn)移通道和賦存空間,有利于金屬礦的形成(李雙等,2012).

    (3) 榔橋巖體為中生代巖漿活動(dòng),而華南普遍存在大規(guī)模中生代巖漿作用和成礦作用.國內(nèi)對(duì)華南中生代巖漿作用已經(jīng)進(jìn)行了大量的研究,中生代發(fā)生的大規(guī)模巖漿活動(dòng),使得華南從南到北形成了三條形成時(shí)代漸新的成礦帶,而從成礦礦物的種類來看,華南成礦帶明顯類似于南美俯沖作用所引起的成礦分帶,因此對(duì)于榔橋巖體的動(dòng)力學(xué)背景可以用洋脊俯沖假說來解釋.

    綜上所述烏溪礦區(qū)的成礦動(dòng)力學(xué)模型基本可以概括為:太平洋板塊在中生代(125~140 Ma)時(shí)向南西方向俯沖,而依澤納吉板塊則向北北西方向俯沖,兩個(gè)板塊的剪切作用方向與同時(shí)期中國東部近南北向盆地拉張的力學(xué)性質(zhì)是耦合的.拉張?jiān)斐蓭r石圈減薄,軟流圈卸載上涌,發(fā)生減壓部分熔融.熔體上涌時(shí)先將古老巖石圈地幔中易熔的富集組分熔融(Ling et al.,2009;Sun et al.,2010),從而形成這一時(shí)期大規(guī)模巖漿活動(dòng)和巖漿的富集特征.榔橋巖體是華南中生代巖漿活動(dòng)中的其中一個(gè)巖漿事件.與榔橋巖體相關(guān)的巖漿活動(dòng)受到區(qū)內(nèi)湯口斷裂的控制,在巖漿上侵過程中到達(dá)礦區(qū)范圍,并在區(qū)內(nèi)可能存在的斷裂空間內(nèi)結(jié)晶分異和礦物富集成礦.通過綜合電磁法結(jié)果推測礦區(qū)范圍內(nèi)在巖漿作用下的成礦物質(zhì)的運(yùn)移方向.磁場分布反映了巖漿的運(yùn)移方向?yàn)閺牡V區(qū)南側(cè)深部沿著深部斷裂到達(dá)北側(cè)的淺部,在深部結(jié)晶分異并固結(jié)成巖.而區(qū)內(nèi)反映礦體賦存空間位置的極化率異常和電阻率異常表明,礦體可能分布的礦區(qū)正中,呈北東向分布.而這一結(jié)果用與磁場分布存在明顯差異,用巖漿活動(dòng)不能很好的解釋極化率和電阻率異常,因此推測礦區(qū)內(nèi)存在后期的熱液活動(dòng).而由于前期的巖漿活動(dòng)的結(jié)晶分異和固結(jié)成巖,造成深部斷裂裂隙減少,隨著熱液不斷上涌,在深部由于熱液所形成的壓力增大,在區(qū)內(nèi)構(gòu)造應(yīng)力的誘發(fā)下會(huì)發(fā)生大量的隱爆作用.而隱爆作用發(fā)生的的同時(shí),深部熱液的溫壓都急劇下降,這種作用對(duì)于成礦非常有利.

    4 深部成礦機(jī)制與鉆孔驗(yàn)證

    根據(jù)礦區(qū)的前期地質(zhì)工作,礦區(qū)內(nèi)主要出露地層為志留系舉坑組、泥盆系五通組、石炭系中下統(tǒng)和第四紀(jì)全新統(tǒng),地層傾角較大.區(qū)內(nèi)發(fā)育一組北東向斷裂構(gòu)造,控制著區(qū)內(nèi)花崗斑巖的分布,同時(shí)區(qū)內(nèi)礦化蝕變也受這組斷裂控制.從而推測區(qū)內(nèi)不僅存在熱液活動(dòng),前期還存在巖漿侵入.測區(qū)位于榔橋巖體北側(cè),區(qū)內(nèi)出露的花崗斑巖脈可以認(rèn)為是榔橋巖體的分支,區(qū)內(nèi)出現(xiàn)的礦化蝕變可以推測為巖漿期后熱液活動(dòng)對(duì)圍巖的改造.

    根據(jù)三種電磁方法的綜合結(jié)果分析,測區(qū)深部較低電阻率范圍較表層大,熱液從深部進(jìn)入,由于溫度壓力較高,可能發(fā)生強(qiáng)烈的隱爆作用,從而使深部發(fā)育較強(qiáng)礦化蝕變.由于蝕變使花崗巖類巖石磁性減弱,同時(shí)熱液活動(dòng)有可能形成磁黃鐵礦,因此出現(xiàn)磁場的不均勻分布.隨著熱液向上運(yùn)移,溫度壓力降低,使熱液開始沿區(qū)內(nèi)的斷裂構(gòu)造分布,在不同的位置不同的圍巖內(nèi)形成不同的蝕變和礦化,因此淺層電阻率異常分布基本與斷裂走向一致.而在遠(yuǎn)離巖體且沒有明顯斷裂活動(dòng)的位置(小號(hào)點(diǎn)),熱液的作用未能將圍巖蝕變和改造,因此還保持了早期巖性,沒有明顯的電磁異常.綜合結(jié)果表明,磁場從北向南逐漸減小的原因并不是由于巖體埋深增加,而是由于蝕變作用使相應(yīng)的磁性物質(zhì)的磁性減弱.

    根據(jù)地質(zhì)地球物理資料,選取成礦條件最有利的7000線(625點(diǎn))和7300線(650點(diǎn))進(jìn)行鉆探驗(yàn)證.7001孔深1250 m,7301孔深1000 m,深部巖層普遍礦化和強(qiáng)烈硅化,部分發(fā)育浸染狀黃鐵礦化,較好的證明了鉆孔高極化率的特征.從電阻率和鉆孔疊加圖上看(圖8),電阻率分布特征與區(qū)內(nèi)巖性分布基本一致,鉆探結(jié)果進(jìn)一步證明區(qū)內(nèi)地表礦化不明顯,整體電阻率較高,巖體相對(duì)其他巖石表現(xiàn)為高阻特征,礦化能夠明顯的降低巖層電阻率.

    從7301鉆孔結(jié)果(圖9)可以看到,花崗斑巖以脈體的形式侵入到上覆志留系地層的砂巖、粉砂巖中,巖脈普遍伴生黃鐵礦化,巖脈附近泥化破碎帶內(nèi)多金屬礦物更加富集.淺部發(fā)育青磐巖化、絹云母化,往深部各類巖石強(qiáng)烈蝕變,同時(shí)巖心節(jié)理發(fā)育,節(jié)理面普遍存在礦化,深部斷續(xù)可見隱爆角礫巖和碳酸鹽化,說明在深部存在強(qiáng)烈的熱液活動(dòng).從地下400 m深度開始,巖石普遍發(fā)育黃鐵礦化、鉛鋅礦化,隨著深度增加出現(xiàn)黃銅礦化,從900 m位置出現(xiàn)明顯的銅、鉬礦化,950 m開始發(fā)育含輝鉬礦的花崗斑巖(見圖10).

    地球化學(xué)分析結(jié)果顯示不管從巖石年代學(xué)還是巖石成因?qū)W分析都顯示烏溪礦區(qū)花崗斑巖和榔橋巖體屬于同一巖漿活動(dòng).烏溪成礦花崗斑巖鉆孔樣品以及地表出露巖體的鋯石定年結(jié)果一致,分別為139.6±1.7 Ma(ZK7301)、137.3±1.6 Ma(ZK7001)、137.3±1.1 Ma(10WX-1),表明該巖體形成時(shí)代為早白堊世.同時(shí)鉆孔結(jié)果顯示花崗斑巖脈明顯受到后期熱液活動(dòng)的蝕變,出現(xiàn)強(qiáng)烈的巖石破碎和泥化,泥化帶礦化明顯,因此推測研究區(qū)主要成礦期并不是巖漿活動(dòng)期,而是巖漿期后熱液活動(dòng)期.

    圖8 電阻率剖面&鉆孔疊加圖Fig.8 Resistivity profile and drilling holes superposition

    圖9 烏溪礦區(qū)7301鉆孔地質(zhì)圖(李雙等,2015)Fig.9 7301 Drilling geological colamns in Wuxi(Li et al.,2015)

    圖10 烏溪礦區(qū)7301孔蝕變巖及礦化照片(959~966 m)Fig.10 Rock alteration and mineralization of 7301 drilling in Wuxi(959~966 m)

    因此綜合各種地質(zhì)地球物理資料結(jié)果,可以推測礦區(qū)深部的成礦機(jī)制:

    (1) 測區(qū)早期存在巖漿活動(dòng),巖漿沿?cái)嗔焉锨郑趲r體周圍形成較弱的圍巖蝕變,并且有少量成礦物質(zhì)的代入.根據(jù)斷裂的走向,推測除了地表出露的花崗斑巖脈,在測區(qū)北側(cè)磁場值較高區(qū)域有大量巖漿侵入,在淺層形成隱伏巖體.

    (2) 區(qū)內(nèi)存在強(qiáng)烈的熱液活動(dòng),地表熱液活動(dòng)主要在斷裂及其附近,并且存在一定的順層侵入,推測深部低阻區(qū)域存在較大范圍的圍巖蝕變和礦化空間.同時(shí)隱爆角礫巖的出現(xiàn)表明深部存在強(qiáng)烈的隱爆作用,膠結(jié)物中發(fā)育網(wǎng)脈狀鉛鋅礦,說明隱爆作用有利于礦物沉積.

    (3) 烏溪金礦礦區(qū)內(nèi)發(fā)育的大量斷裂構(gòu)造為成礦流體提供了充分的運(yùn)移通道,巖漿熱液主要沿通道活動(dòng),雖然在深部可以通過隱爆作用形成熱液的通道,但是隨著熱液運(yùn)移,溫壓條件不斷降低,在沒有斷裂活動(dòng)的位置則保存較為原始的圍巖.推測小號(hào)線整體為沉積層,大號(hào)線較好的保存了前期巖漿活動(dòng)形成的巖體.

    (4) 鉆孔結(jié)果表明區(qū)內(nèi)淺部發(fā)育青磐巖化和絹云母化,深部發(fā)育強(qiáng)烈蝕變斑巖和Cu-Pb-Zn礦化,結(jié)合地球物理結(jié)果,在礦區(qū)深部應(yīng)大范圍的發(fā)育蝕變斑巖和礦化,因此初步判定烏溪礦區(qū)為斑巖型多金屬礦床.

    進(jìn)一步的野外勘測以及地球化學(xué)工作將對(duì)烏溪礦區(qū)及周邊探礦和找礦工作具有重要的指示意義.

    5 結(jié)論

    5.1 在研究區(qū)使用激發(fā)極化法、地面磁測和CSAMT電阻率測深法,能對(duì)地表礦化作用的分布、磁性異常地層分布和深部地層的電阻率分布有一個(gè)整體的和立體的描述.結(jié)合已知的地表地質(zhì)、地球化學(xué)、地表礦化點(diǎn)和地表斑巖脈的分布等信息,可以推斷地下深部隱伏有含礦的斑巖體.根據(jù)這一思路和地球物理探測結(jié)果布設(shè)的兩個(gè)鉆井,揭露了隱伏的斑巖體礦體的存在,驗(yàn)證了推斷正確性和地球物理方法的有效性.這一套思路和地球物理勘探方法,可以在皖南和長江中下游類似的條件下借鑒使用.

    5.2 通過綜合電磁方法的應(yīng)用,揭示了礦區(qū)內(nèi)巖體和蝕變的分布特征及相互關(guān)系,進(jìn)一步推測成礦物質(zhì)運(yùn)移方向和深部成礦機(jī)制.結(jié)合地球化學(xué)和鉆孔地質(zhì)資料,推測研究區(qū)主要成礦期并不是巖漿活動(dòng)期,而是巖漿期后熱液活動(dòng)期.因此遠(yuǎn)離榔橋巖體的高磁異常對(duì)應(yīng)的高阻體推測為前期巖漿活動(dòng)的產(chǎn)物.

    5.3 鉆井驗(yàn)證研究區(qū)深部成礦受斑巖控制,礦床類型為斑巖型銅(鉬金)礦床.后期熱液活動(dòng)將圍巖和前期巖漿活動(dòng)產(chǎn)物進(jìn)一步蝕變,使成礦物質(zhì)進(jìn)一步富集,從而深部形成大面積的低阻異常區(qū)域,同時(shí)與地面極化率異常較好對(duì)應(yīng).根據(jù)巖體的極化率和低阻體的規(guī)模推斷,烏溪礦區(qū)是一個(gè)極其有利的斑巖型銅(鉬金)成礦遠(yuǎn)景區(qū),這也進(jìn)一步擴(kuò)大了有俯沖作用形成的華南成礦分帶中斑巖型銅礦的可能范圍,有利于引導(dǎo)華南的區(qū)域成礦地質(zhì)研究.

    致謝 中國科技大學(xué)鄧江紅博士和中科院廣州地球化學(xué)研究所孫賽君博士,在論文寫作過程中提供了相關(guān)地質(zhì)和地球化學(xué)資料;中科院地質(zhì)與地球物理研究所王若老師和王妙月老師在文章修改過程中提出了許多寶貴的修改建議;審稿人提出了多處中肯的建議,對(duì)文章做出了有益的貢獻(xiàn),在此一并致以衷心的感謝.

    References

    Bastani M,Malehmir A,Ismail N,et al.2009.Delineating hydrothermal stockwork copper deposits using controlled-source and radio-magnetotelluric methods:A case study from northeast Iran.Geophysics,74(5):B167-B181.

    Chen Y L,Han K,Chen Y X,et al.2015.The application of CSAMT in Karst collapse investigation.Progress in Geophysics (in Chinese),30(6):2616-2622,doi:10.6038/pg20150620.

    Cooke D R,Hollings P,Walshe J L.2005.Giant porphyry deposits:characteristics,distribution,and tectonic controls.Economic Geology,100(5):801-818.

    Cui M L,Zhang B L,Liang G H,et al.2010.The technical combination of comprehensive geophysical prospecting in the Molybdenum mines with loess-covered:A case study at the shapoling Molybdenum mine.Progress in Geophysics (in Chinese),25(2):602-611,doi:10.3969/j.issn.1004-2903.2010.02.033.

    Deng J Z,Chen H,Yin C C,et al.2015.Three-dimensional electrical structures and significance for mineral exploration in the Jiujiang-Ruichang District.Chinese J.Geophys.(in Chinese),58(12):4465-4477,doi:10.6038/cjg20151211.

    Ding G M,Zhu Z Q,Chang R F,et al.2015.Application of integrative geophysical methods to predicting deep ore body in Yechangping molybdenum deposit,Henan province.Progress in Geophysics (in Chinese),30(1):325-331,doi:10.6038/pg20150147.

    Hao T Y,Jiang W W.1998.Application of comprehensive geophysical methods in looking for hidden gold mine in Bailidian region.Acta Geophysica Sinica (in Chinese),41(S1):404-413.

    Li D Q,Di Q Y,Wang G J,et al.2008.Fault detection by CSAMT and its application to new district planning in Beijing.Progress in Geophysics (in Chinese),23(6):1963-1969.

    Ling M X,Wang F Y,Ding X,et al.2009.Cretaceous ridge subduction along the Lower Yangtze River Belt,eastern China.Economic Geology,104(2):303-321.

    Liu S A,Li S G,He Y S,et al.2010.Geochemical contrasts between early Cretaceous ore-bearing and ore-barren high-Mg adakites in central-eastern China:Implications for petrogenesis and Cu-Au mineralization.Geochimica et Cosmochimica Acta,74(24):7160-7178.

    Lou D B,Song G X,Li N,et al.2008.The application of magnetic method in national mineral prediction.Progress in Geophysics (in Chinese),23(1):249-256.

    Lv Q T,Dong S W,Tang J T,et al.2015.Multi-scale and integrated geophysical data revealing mineral systems and exploring for mineral deposits at depth:A synthesis from SinoProbe-03.Chinese Journal of Geophysics (in Chinese),58(12):4319-4343,doi:10.6038/cjg20151201.

    Shao L S,Liu Z D,Lv Q T,et al.2015.Deep fine structure of Guichi Ore concentrated area:The understanding of the integrated geophysical detection results.Chinese Journal of Geophysics (in Chinese),58(12):4490-4504,doi:10.6038/cjg20151213.

    Smith R.2014.Electromagnetic induction methods in mining geophysics from 2008 to 2012.Surveys in Geophysics,35(1):123-156.

    Sun W D,Xie Z,Chen J F,et al.2003.Os-Os dating of copper and molybdenum deposits along the middle and lower reaches of the Yangtze River,China.Economic Geology,98(1):175-180.

    Sun W D,Ding X,Hu Y H,et al.2007.The golden transformation of the Cretaceous plate subduction in the west Pacific.Earth and Planetary Science Letters,262(3-4):533-542.

    Sun W D,Ling M X,Yang X Y,et al.2010.Ridge subduction and porphyry copper-gold mineralization:An overview.Science China-Earth Sciences,53(4):475-484.

    Sun W D,Liang H Y,Ling M X,et al.2013.The link between reduced porphyry copper deposits and oxidized magmas.Geochimica Et Cosmochimica Acta,103:263-275.

    Sun B,Li T L,Li H,et al.2015.Study on the sounding of CSAMT.Progress in Geophysics (in Chinese),30(2):836-839,doi:10.6038/pg20150247.

    Sun X G,Liu J M,Liu H T,et al.2007.The application of integrated geophysical prospecting method to the evaluation of Haolibao copper deposits.Progress in Geophysics (in Chinese),22(6):1910-1915.

    Tang J T,Liu Z J,Liu F Y,et al.2015.The denoising of the audio magnetotelluric data set with strong interferences.Chinese Journal of Geophysics (in Chinese),58(12):4636-4647,doi:10.6038/cjg20151225.

    Wan H P,Xu G L,Yu X,et al.2015.The geophysics evaluation on granite type of high-level radioactive waste disposal repository——taking Aqishan potential site as an example.Progress in Geophysics (in Chinese),30(6):2778-2784,doi:10.6038/pg20150642.

    Wang F Y,Liu S A,Li S G,et al.2013.Contrasting zircon Hf-O isotopes and trace elements between ore-bearing and ore-barren adakitic rocks in central-eastern China:implications for genetic relation to Cu-Au mineralization.Lithos,156-159:97-111.

    Wang R,Yin C C,Wang M Y,et al.2014.CSAMT sensitivity analysis for 1D models.Progress in Geophysics (in Chinese),29(3):1284-1291,doi:10.6038/pg20140339.

    Wang X Y,Tang J T,Zhang L C,et al.2015.Lithospheric electrical structure in the middle and lower reach of Yangtze River metallogenic belt inferred from magnetotelluric sounding.Chinese Journal of Geophysics (in Chinese),58(12):4403-4414,doi:10.6038/cjg20151206.

    Wu G X.2007.High accuracy magnetic prospecting on earth surface in gold mine exploration——Take the Shisangongli exploring area out of Wulaga gold deposit in Heilongjiang province as an example.Progress in Geophysics (in Chinese),22(5):1637-1641.

    Xu X Q,Su L J,Liang S Q.2015.A review of geophysical detection methods of landslide structure characteristics.Progress in Geophysics (in Chinese),30(3):1449-1458,doi:10.6038/pg20150361.

    Yang Z W,Yan J Y,Chen X B.2013.The application of spectral induced plarization in Shaxi porphyry copper in Anhui Province.Progress in Geophysics (in Chinese),28(4):2014-2023,doi:10.6038/pg20130445.

    Yu C M.1998.The application of CSAMT method in looking for hidden gold mine.Acta Geophysica Sinica (in Chinese),41(1):133-138.

    Zhang G Z,Zhou L G,Wang Y H.2015.Application of integrated electrical mathods to silver lead-zinc mine zone of Chaaobao in Inner Mongolia.Progress in Geophysics (in Chinese),30(2):867-871,doi:10.6038/pg20150252.

    Zhang K,Yan J Y,Lv Q T,et al.2015.The study of upper crust electrical structure in Ningwu and adjacent area.Chinese Journal of Geophysics (in Chinese),58(12):4505-4521,doi:10.6038/cjg20151214.

    Zhang Z,Liu J M,Yu C M,et al.2013.Application of integrated geophysical prospecting methods in the evaluation of BIF deposits—a case study in Inner Mongolia Aohanqi Sijiazi BIF deposits.Progress in Geophysics (in Chinese),28(4):2078-2084,doi:10.6038/pg20130452.

    Zhou C,Tang J T,Ren Z Y,et al.2015.Application of the Rhoplus method to audio magnetotelluric dead band distortion data.Chinese Journal of Geophysics (in Chinese),58(12):4648-4660,doi:10.6038/cjg20151226.

    附中文參考文獻(xiàn)

    安徽省地質(zhì)礦產(chǎn)局.1987.安徽省區(qū)域地質(zhì)志.北京:地質(zhì)出版社.

    陳偉軍,劉紅濤,劉建明等.2008.隱伏礦床定位預(yù)測常用的綜合物化探方法.有色礦冶,24(2):2-7.

    陳玉玲,韓凱,陳貽祥等.2015.可控源音頻大地電磁法在巖溶塌陷勘察中的應(yīng)用.地球物理學(xué)進(jìn)展,30(6):2616-2622,doi:10.6038/pg20150620.

    崔敏利,張寶林,梁光河等.2010.黃土覆蓋區(qū)鉬礦綜合地球物理找礦技術(shù)組合:以沙坡嶺鉬礦為例.地球物理學(xué)進(jìn)展,25(2):602-611,doi:10.3969/j.issn.1004-2903.2010.02.033.

    鄧居智,陳輝,殷長春等.2015.九瑞礦集區(qū)三維電性結(jié)構(gòu)研究及找礦意義.地球物理學(xué)報(bào),58(12):4465-4477,doi:10.6038/cjg20151211.

    丁高明,朱自強(qiáng),常榮鳳等.2015.綜合物探法在河南夜長坪鉬礦深部找礦預(yù)測中的應(yīng)用.地球物理學(xué)進(jìn)展,30(1):325-331,doi:10.6038/pg20150147.

    郝天珧,江為為.1998.綜合地球物理方法在山東百里店地區(qū)找尋隱伏金礦中的應(yīng)用.地球物理學(xué)報(bào),41(S1):404-413.

    何繼善.2006.雙頻激電法.北京:高等教育出版社.

    蘭學(xué)毅,周存亭,王建偉等.2012.安徽省重力異常特征分區(qū)與地質(zhì)構(gòu)造單元?jiǎng)澐?安徽地質(zhì),22(1):1-8.

    李帝銓,底青云,王光杰等.2008.CSAMT探測斷層在北京新區(qū)規(guī)劃中的應(yīng)用.地球物理學(xué)進(jìn)展,23(6):1963-1969.

    李雙,孫衛(wèi)東,楊曉勇等.2012.烏溪金礦蝕變分帶及元素地球化學(xué)特征.地球科學(xué)進(jìn)展,27(S1):34-35.

    李雙,楊曉勇,孫衛(wèi)東等.2014.皖南涇縣榔橋巖體鋯石U-Pb定年、Hf同位素和地球化學(xué)特征及其找礦指示意義.地質(zhì)學(xué)報(bào),88(8):1561-1578.

    李雙,孫賽軍,楊曉勇等.2015.皖南烏溪斑巖型金礦床賦礦侵入巖體的巖石地球化學(xué)及年代學(xué)研究.大地構(gòu)造與成礦學(xué),39(1):153-166.

    劉琛琛,楊錢江.2012.安徽省涇縣烏溪金礦地質(zhì)特征及控礦因素淺析.科技視界,(25):343-344,210.

    劉惠華.2003.烏溪金礦成礦后期構(gòu)造分析在坑道探礦中的作用.西部探礦工程,15(9):65-66.

    劉惠華,朱寧.2004.皖南涇縣烏溪地區(qū)金礦成礦條件分析.安徽地質(zhì),14(1):30-32.

    柳建新,劉春明,佟鐵鋼等.2004.雙頻激電法在西藏某銅多金屬礦帶的應(yīng)用.地質(zhì)與勘探,40(2):59-61.

    婁德波,宋國璽,李楠等.2008.磁法在我國礦產(chǎn)預(yù)測中的應(yīng)用.地球物理學(xué)進(jìn)展,23(1):249-256.

    盧焱,李健,白雪山等.2008.地面磁法在隱伏鐵礦勘查中的應(yīng)用—以河北灤平Ⅱ號(hào)鐵礦為例.吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),38(4):698-702.

    呂慶田,董樹文,湯井田等.2015.多尺度綜合地球物理探測:揭示成礦系統(tǒng)、助力深部找礦——長江中下游深部探測(SinoProbe-03)進(jìn)展.地球物理學(xué)報(bào),58(12):4319-4343,doi:10.6038/cjg20151201.

    邵陸森,劉振東,呂慶田等.2015.安徽貴池礦集區(qū)深部精細(xì)結(jié)構(gòu)——來自綜合地球物理探測結(jié)果的認(rèn)識(shí).地球物理學(xué)報(bào),58(12):4490-4504,doi:10.6038/cjg20151213.

    時(shí)彬.2012.CSAMT在深部礦產(chǎn)勘查中的研究與應(yīng)用[碩士論文].長春:吉林大學(xué).

    孫博,李桐林,李鶴等.2015.可控源音頻大地電磁測深法勘查深度研究.地球物理學(xué)進(jìn)展,30(2):836-839,doi:10.6038/pg20150247.

    孫興國,劉建明,劉洪濤等.2007.綜合物探方法在好力寶銅礦床的應(yīng)用.地球物理學(xué)進(jìn)展,22(6):1910-1915.

    譚章坤.2013.CSAMT在深部勘探中的效果研究[碩士論文].成都:成都理工大學(xué).

    湯井田,劉子杰,劉峰屹等.2015.音頻大地電磁法強(qiáng)干擾壓制試驗(yàn)研究.地球物理學(xué)報(bào),58(12):4636-4647,doi:10.6038/cjg20151225.

    萬芬,王光杰,楊曉勇等.2014.皖南云嶺金礦地球物理找礦模型.地球科學(xué)與環(huán)境學(xué)報(bào),36(1):185-192.

    萬芬.2014.皖南烏溪金礦綜合地球物理方法找礦模式研究[碩士論文].武漢:中國地質(zhì)大學(xué)(武漢).

    萬漢平,徐貴來,喻翔等.2015.花崗巖型高放廢物處置庫選址中的地球物理評(píng)價(jià)——以新疆阿奇山地段候選場址為例.地球物理學(xué)進(jìn)展,30(6):2778-2784,doi:10.6038/pg20150642.

    王建偉,李仁和,胡開勇等.2009.安徽省磁場基本特征初析.安徽地質(zhì),19(4):268-271.

    王若,殷長春,王妙月等.2014.CSAMT法一維層狀介質(zhì)靈敏度分析.地球物理學(xué)進(jìn)展,29(3):1284-1291,doi:10.6038/pg20140339.王顯瑩,湯井田,張林成等.2015.長江中下游成礦帶中段巖石圈電性結(jié)構(gòu)研究.地球物理學(xué)報(bào),58(12):4403-4414,doi:10.6038/cjg20151206.

    吳國學(xué).2007.金礦勘查中的地面高精度磁法測量—以黑龍江烏拉嘎金礦外圍十三公里勘查區(qū)為例.地球物理學(xué)進(jìn)展,22(5):1637-1641.

    武煒,張寶林,梁光河等.2009.雙頻激電法在我國西部兩類典型覆蓋區(qū)金屬礦體預(yù)測中的應(yīng)用.地質(zhì)與勘探,45(6):669-675.

    徐興倩,蘇立君,梁雙慶.2015.地球物理方法探測滑坡體結(jié)構(gòu)特征研究現(xiàn)狀綜述.地球物理學(xué)進(jìn)展,30(3):1449-1458,doi:10.6038/pg20150361.

    楊振威,嚴(yán)加永,陳向斌.2013.頻譜激電法在安徽沙溪斑巖銅礦中的應(yīng)用.地球物理學(xué)進(jìn)展,28(4):2014-2023,doi:10.6038/pg20130445.

    于昌明.1998.CSAMT方法在尋找隱伏金礦中的應(yīng)用.地球物理學(xué)報(bào),41(1):133-138.

    張光之,周立國,王延浩.2015.綜合電法在內(nèi)蒙古查敖包銀鉛鋅礦區(qū)的應(yīng)用.地球物理學(xué)進(jìn)展,30(2):867-871,doi:10.6038/pg20150252.

    張昆,嚴(yán)加永,呂慶田等.2015.寧蕪火山巖盆地及鄰區(qū)上地殼電性結(jié)構(gòu)研究.地球物理學(xué)報(bào),58(12):4505-4521,doi:10.6038/cjg20151214.

    張壯,劉建明,于昌明等.2013.綜合地球物理方法在鞍山式鐵礦勘查中的應(yīng)用—以內(nèi)蒙古敖漢旗四家子鐵礦為例.地球物理學(xué)進(jìn)展,28(4):2078-2084,doi:10.6038/pg20130452.

    趙永利,史春旺,劉惠華.2013.安徽省涇縣烏溪金礦地質(zhì)特征及控礦因素.西部探礦工程,25(8):106-108.

    周聰,湯井田,任政勇等.2015.音頻大地電磁法“死頻帶”畸變數(shù)據(jù)的Rhoplus校正.地球物理學(xué)報(bào),58(12):4648-4660,doi:10.6038/cjg20151226.

    (本文編輯 劉少華)

    Two-dimensional joint inversion of seismic velocity and electrical resistivity using seismic travel times and full channel electrical measurements based on alternating cross-gradient structural constraint

    GAO Ji1,2,ZHANG Hai-Jiang1*

    1 Laboratory of Seismology and Physics of Earth′s Interior,School of Earth and Space Sciences, University of Science and Technology of China,Hefei 230026,China2 Anhui Wantai Geophysical Technology Co.Ltd.,Hefei 230026,China

    When imaging complex near surface structures using different geophysical exploration methods,because of the incompleteness of the observation system and insensitivity of geophysical data to some geophysical properties of media,there exist relatively large uncertainties and significant inconsistency in imaging results from separate inversions.For seismic travel time tomography and electrical resistivity tomography,both methods are subject to the issue of imaging shadow zones.For seismic travel time tomography,seismic rays tend to travel through high velocity zones,causing poor ray path coverage and thus low resolution of low velocity zones.For electrical resistivity tomography,the distribution of electric field lines is relatively sparse in high resistivity zones,leading to lower resolution.In order to reliably image underground structure,Gallardo and Meju (2003) proposed a joint geophysical inversion method based on cross-gradient structure constraint,which requires separate models fitting the data and at the same time structurally consistent by having approximately zero cross gradients.To more effectively realize the cross-gradient based structure constraint,we propose a new joint inversion scheme that alternatively inverts one type of data but structurally constrained by another model.With this new joint inversion scheme,we test two-dimensional cross-hole seismic travel time tomography and full channel resistivity tomography.Synthetic results show that the new scheme can efficiently realize the cross-gradient based joint inversion.Compared to separate inversions,the joint inversion can more reliably determine seismic velocity and electrical resistivity anomalies.

    Seismic travel time tomography;Electrical resistivity tomography;Structure similarity;Cross-gradient;Joint inversion

    Application of comprehensive electromagnetic study in deep mineralization mechanism—A case study of the Wuxi polymetallic ore deposit,south Anhui

    LIN Fang-Li1,2,WANG Guang-Jie1,YANG Xiao-Yong3

    1 Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing 100029,China2 University of Chinese Academy of Sciences,Beijing 100049,China3 School of Earth and Space Sciences,University of Science and Technology of China,Hefei 230036,China

    高級(jí),張海江.2016.基于交叉梯度交替結(jié)構(gòu)約束的二維地震走時(shí)與全通道直流電阻率聯(lián)合反演.地球物理學(xué)報(bào),59(11):4310-4322,

    10.6038/cjg20161131.

    Gao J,Zhang H J.2016.Two-dimensional joint inversion of seismic velocity and electrical resistivity using seismic travel times and full channel electrical measurements based on alternating cross-gradient structural constraint.Chinese J.Geophys.(in Chinese),59(11):4310-4322,doi:10.6038/cjg20161131.

    國家自然科學(xué)基金項(xiàng)目(41474039)資助.

    高級(jí),男,1983年生,博士,主要從事地震資料處理、彈性波層析成像、電阻率層析成像、多屬性聯(lián)合反演研究.E-mail:gaoji617@mail.ustc.edu.cn

    *通訊作者 張海江,男,1973年生,教授、博士生導(dǎo)師,主要從事發(fā)展新的地震成像算法和聯(lián)合地球物理成像算法以及對(duì)斷層、火山、俯沖帶結(jié)構(gòu)的成像,油氣礦產(chǎn)開發(fā)過程誘發(fā)的微地震等研究.E-mail:zhang11@ustc.edu.cn

    10.6038/cjg20161131

    P631 收稿日期2015-11-23,2016-09-22收修定稿

    ?? 中圖分類號(hào) P631

    2015-11-18,2016-05-20收修定稿

    猜你喜歡
    斑巖
    玲瓏金礦田煌斑巖與礦脈關(guān)系的探索及應(yīng)用
    東天山赤湖地區(qū)原生暈異常結(jié)構(gòu)特征對(duì)尋找斑巖型銅鉬礦床的指示意義
    遼寧調(diào)兵山西調(diào)斑巖型鉬礦床特征及找礦標(biāo)志
    柴北緣阿木尼克山地區(qū)斑巖系Cu、Mo-Pb、Zn、Ag-Au成礦模型初步研究
    斑巖型礦床含礦斑巖與非含礦斑巖鑒定特征綜述
    黑云母和鋯石化學(xué)組分對(duì)巖漿結(jié)晶條件的約束:以滇西北衙超大型金礦床為例*
    巖型礦床含礦斑巖與非含礦斑巖鑒定特征綜述
    煌斑巖的研究進(jìn)展
    斑巖型銅礦床研究現(xiàn)狀與進(jìn)展
    山東煙臺(tái)牟乳地區(qū)煌斑巖與金礦成礦關(guān)系
    地球(2016年7期)2016-04-14 22:00:20
    日本爱情动作片www.在线观看| 一级毛片黄色毛片免费观看视频| 熟妇人妻不卡中文字幕| 成人无遮挡网站| 久久青草综合色| 人人妻人人添人人爽欧美一区卜| 色视频在线一区二区三区| 日本爱情动作片www.在线观看| 美女福利国产在线| 丝袜美足系列| 欧美日韩国产mv在线观看视频| 国产精品嫩草影院av在线观看| 日本欧美国产在线视频| 日韩 亚洲 欧美在线| 国产精品一区www在线观看| av不卡在线播放| 在线 av 中文字幕| 亚洲成人手机| 国产色婷婷99| 亚洲丝袜综合中文字幕| 国产精品久久久久久久久免| 亚洲精品乱码久久久v下载方式| 免费人成在线观看视频色| 中文字幕人妻丝袜制服| 精品人妻一区二区三区麻豆| 丝袜脚勾引网站| 久久久精品区二区三区| 老司机亚洲免费影院| 亚洲av不卡在线观看| 久久这里有精品视频免费| 91精品三级在线观看| 亚洲经典国产精华液单| 制服诱惑二区| 久久久久久伊人网av| 啦啦啦在线观看免费高清www| 制服人妻中文乱码| 成人18禁高潮啪啪吃奶动态图 | 大片电影免费在线观看免费| 一本—道久久a久久精品蜜桃钙片| av免费观看日本| 亚洲欧美日韩另类电影网站| 国产老妇伦熟女老妇高清| 男人操女人黄网站| 欧美人与善性xxx| 极品少妇高潮喷水抽搐| 亚洲婷婷狠狠爱综合网| 亚洲综合色惰| 看非洲黑人一级黄片| 亚洲精品日韩av片在线观看| 日韩三级伦理在线观看| 大片电影免费在线观看免费| 天美传媒精品一区二区| 嫩草影院入口| 免费日韩欧美在线观看| 亚洲国产精品成人久久小说| 久久久久国产精品人妻一区二区| 青春草国产在线视频| 美女视频免费永久观看网站| 高清毛片免费看| 一级爰片在线观看| 免费看不卡的av| 九九在线视频观看精品| 日韩一本色道免费dvd| 91aial.com中文字幕在线观看| av福利片在线| 99久久中文字幕三级久久日本| 一级毛片电影观看| 一级毛片电影观看| 男女国产视频网站| 伦理电影免费视频| 99九九线精品视频在线观看视频| 亚洲色图综合在线观看| 99久久综合免费| 国产亚洲精品久久久com| 国产片特级美女逼逼视频| 91精品国产九色| 久久狼人影院| 国产不卡av网站在线观看| 波野结衣二区三区在线| 日韩一本色道免费dvd| 国产极品粉嫩免费观看在线 | 久久人人爽人人爽人人片va| 多毛熟女@视频| 亚洲精品av麻豆狂野| 精品卡一卡二卡四卡免费| 97精品久久久久久久久久精品| 超碰97精品在线观看| 丁香六月天网| 国产不卡av网站在线观看| 免费观看在线日韩| av福利片在线| 日韩av免费高清视频| 国产免费视频播放在线视频| 国产精品一区www在线观看| av黄色大香蕉| 亚洲,欧美,日韩| 日本黄大片高清| 日本av免费视频播放| 国产男人的电影天堂91| 免费大片18禁| 老司机影院成人| 妹子高潮喷水视频| 日本91视频免费播放| 精品一品国产午夜福利视频| 精品人妻熟女毛片av久久网站| 男女免费视频国产| 少妇被粗大猛烈的视频| 亚洲av福利一区| 91精品国产国语对白视频| 亚洲av免费高清在线观看| 美女视频免费永久观看网站| 欧美一级a爱片免费观看看| 亚洲久久久国产精品| 国产av精品麻豆| 国产永久视频网站| 久久精品久久久久久噜噜老黄| 高清av免费在线| 久久久久久久精品精品| 中文乱码字字幕精品一区二区三区| 中文字幕亚洲精品专区| 久久久精品免费免费高清| 妹子高潮喷水视频| 大香蕉97超碰在线| 国产精品国产三级国产专区5o| 美女视频免费永久观看网站| 秋霞伦理黄片| 久久久久久久久久久免费av| 亚洲精品自拍成人| 美女国产视频在线观看| 国产亚洲最大av| 国产成人免费无遮挡视频| 老司机亚洲免费影院| 欧美日韩亚洲高清精品| 国产色婷婷99| 成人二区视频| 熟妇人妻不卡中文字幕| 涩涩av久久男人的天堂| 夜夜骑夜夜射夜夜干| 又大又黄又爽视频免费| 日韩一本色道免费dvd| 成人毛片60女人毛片免费| 男的添女的下面高潮视频| 满18在线观看网站| 免费观看a级毛片全部| 亚洲综合色网址| 最近中文字幕高清免费大全6| 国产老妇伦熟女老妇高清| 欧美日韩亚洲高清精品| 国产淫语在线视频| videosex国产| 欧美日韩成人在线一区二区| 国产男人的电影天堂91| 中文字幕av电影在线播放| 最近的中文字幕免费完整| 久久精品国产鲁丝片午夜精品| 久久久久久久久久成人| 亚洲av欧美aⅴ国产| 国产成人91sexporn| 在线观看美女被高潮喷水网站| 曰老女人黄片| 国产亚洲精品第一综合不卡 | 欧美人与善性xxx| 日韩av在线免费看完整版不卡| 精品久久蜜臀av无| 99久久综合免费| 热re99久久精品国产66热6| 亚洲第一av免费看| .国产精品久久| 看非洲黑人一级黄片| 9色porny在线观看| 大香蕉97超碰在线| 免费黄网站久久成人精品| 亚洲性久久影院| 男的添女的下面高潮视频| 成人国语在线视频| 中文字幕av电影在线播放| 国产精品国产三级国产专区5o| 一级爰片在线观看| 日韩熟女老妇一区二区性免费视频| 国产极品天堂在线| 亚洲精品国产色婷婷电影| 成年人午夜在线观看视频| 免费黄频网站在线观看国产| 成人亚洲精品一区在线观看| 人人妻人人澡人人看| 亚洲色图综合在线观看| 欧美丝袜亚洲另类| 精品人妻在线不人妻| 亚洲第一区二区三区不卡| 我的女老师完整版在线观看| 高清不卡的av网站| 女的被弄到高潮叫床怎么办| 中国美白少妇内射xxxbb| 亚洲欧美精品自产自拍| 久久久久视频综合| 丰满少妇做爰视频| 久久精品久久久久久久性| 午夜视频国产福利| 三级国产精品欧美在线观看| 亚洲国产av新网站| 日日摸夜夜添夜夜爱| 免费看av在线观看网站| 热99久久久久精品小说推荐| 母亲3免费完整高清在线观看 | 欧美日韩亚洲高清精品| 国产毛片在线视频| 久久久久久久久久久久大奶| 精品久久蜜臀av无| 国产在视频线精品| 纯流量卡能插随身wifi吗| kizo精华| a级毛片在线看网站| 亚洲欧洲日产国产| 精品卡一卡二卡四卡免费| 成人二区视频| 亚洲欧美成人精品一区二区| 日韩av不卡免费在线播放| 又粗又硬又长又爽又黄的视频| 夜夜骑夜夜射夜夜干| 国产欧美日韩一区二区三区在线 | 色哟哟·www| 春色校园在线视频观看| av在线老鸭窝| 欧美日韩综合久久久久久| 国产伦精品一区二区三区视频9| 黄色怎么调成土黄色| 欧美丝袜亚洲另类| 欧美老熟妇乱子伦牲交| 男女啪啪激烈高潮av片| 国产欧美另类精品又又久久亚洲欧美| 好男人视频免费观看在线| 国产欧美日韩一区二区三区在线 | 久久影院123| 午夜激情av网站| 精品视频人人做人人爽| 国产免费福利视频在线观看| 欧美+日韩+精品| av网站免费在线观看视频| 国产一区二区在线观看日韩| 少妇的逼好多水| 中文乱码字字幕精品一区二区三区| 精品亚洲成国产av| 国产在线免费精品| www.色视频.com| av不卡在线播放| 日韩欧美一区视频在线观看| 亚洲精品视频女| 男女边吃奶边做爰视频| 欧美精品亚洲一区二区| 中文字幕久久专区| 亚洲av日韩在线播放| 中文精品一卡2卡3卡4更新| a 毛片基地| 午夜免费男女啪啪视频观看| 国产精品成人在线| 一级毛片电影观看| 午夜激情福利司机影院| 国产一区有黄有色的免费视频| 大陆偷拍与自拍| 一区二区av电影网| 亚洲精品国产av成人精品| 一级爰片在线观看| 成人毛片a级毛片在线播放| 2018国产大陆天天弄谢| 欧美 亚洲 国产 日韩一| 久久精品国产鲁丝片午夜精品| 久久国产精品大桥未久av| 国产一区二区在线观看av| 好男人视频免费观看在线| 亚洲少妇的诱惑av| 夜夜看夜夜爽夜夜摸| 成人国语在线视频| 男人爽女人下面视频在线观看| 久久久a久久爽久久v久久| 国产成人精品无人区| 多毛熟女@视频| av女优亚洲男人天堂| 国产日韩欧美在线精品| 天天躁夜夜躁狠狠久久av| 亚洲精品国产色婷婷电影| 亚洲精品成人av观看孕妇| 午夜视频国产福利| 亚洲精品亚洲一区二区| 综合色丁香网| 在线亚洲精品国产二区图片欧美 | 国精品久久久久久国模美| 亚洲国产精品国产精品| 夫妻性生交免费视频一级片| 女人精品久久久久毛片| 精品一区二区三区视频在线| 国产精品国产av在线观看| 蜜臀久久99精品久久宅男| 午夜免费鲁丝| 久久久久久久国产电影| 大陆偷拍与自拍| 亚洲精品色激情综合| 少妇人妻精品综合一区二区| 亚洲精品美女久久av网站| 亚洲国产av新网站| 99久久精品国产国产毛片| 中文字幕久久专区| 久久99精品国语久久久| 日日爽夜夜爽网站| 乱码一卡2卡4卡精品| 国产成人精品久久久久久| 国产一区亚洲一区在线观看| 国产黄色视频一区二区在线观看| 少妇被粗大的猛进出69影院 | 高清视频免费观看一区二区| 亚洲国产精品专区欧美| 一区二区三区乱码不卡18| 国产日韩欧美亚洲二区| 欧美xxⅹ黑人| √禁漫天堂资源中文www| 啦啦啦啦在线视频资源| 久久鲁丝午夜福利片| 搡女人真爽免费视频火全软件| 春色校园在线视频观看| 秋霞伦理黄片| 国产精品成人在线| 大片免费播放器 马上看| 日本91视频免费播放| 亚洲一区二区三区欧美精品| 亚洲av电影在线观看一区二区三区| 亚洲精品中文字幕在线视频| 亚洲成人av在线免费| 日本vs欧美在线观看视频| 男人爽女人下面视频在线观看| 亚洲性久久影院| 国产精品成人在线| 免费观看无遮挡的男女| 老女人水多毛片| 亚洲精品自拍成人| 国产黄片视频在线免费观看| 亚洲精品乱码久久久v下载方式| 边亲边吃奶的免费视频| 美女xxoo啪啪120秒动态图| 亚洲精品国产av蜜桃| 天堂中文最新版在线下载| 黑人猛操日本美女一级片| 国产视频内射| 亚洲av二区三区四区| 久久久国产一区二区| 人妻系列 视频| 日韩熟女老妇一区二区性免费视频| 桃花免费在线播放| 狂野欧美白嫩少妇大欣赏| 欧美日本中文国产一区发布| av免费在线看不卡| 亚洲国产精品成人久久小说| 国产不卡av网站在线观看| 国产一区二区在线观看av| 性色avwww在线观看| 一区二区日韩欧美中文字幕 | 欧美日韩综合久久久久久| 色哟哟·www| 狠狠精品人妻久久久久久综合| 大话2 男鬼变身卡| 一级二级三级毛片免费看| 十八禁高潮呻吟视频| 少妇 在线观看| 久久精品人人爽人人爽视色| 亚洲内射少妇av| 特大巨黑吊av在线直播| 男女免费视频国产| 校园人妻丝袜中文字幕| 国产日韩欧美亚洲二区| 亚洲精品亚洲一区二区| 久久久久久久精品精品| 久久久久久久大尺度免费视频| 精品熟女少妇av免费看| 美女xxoo啪啪120秒动态图| 久久精品国产鲁丝片午夜精品| 日韩av不卡免费在线播放| 日日摸夜夜添夜夜爱| 免费人成在线观看视频色| 又大又黄又爽视频免费| 最近中文字幕高清免费大全6| 十八禁高潮呻吟视频| 国产永久视频网站| 黄色怎么调成土黄色| 在线观看一区二区三区激情| 69精品国产乱码久久久| 大香蕉97超碰在线| 亚洲欧美成人精品一区二区| 日本wwww免费看| 97精品久久久久久久久久精品| 十八禁网站网址无遮挡| 亚洲精品色激情综合| 欧美3d第一页| 永久网站在线| 人人澡人人妻人| 十八禁高潮呻吟视频| 欧美日韩成人在线一区二区| 亚洲国产最新在线播放| 中国美白少妇内射xxxbb| 国产精品免费大片| 中文字幕免费在线视频6| 欧美精品人与动牲交sv欧美| 国产成人免费无遮挡视频| 蜜桃在线观看..| 欧美精品国产亚洲| 亚洲美女搞黄在线观看| 亚洲精品乱码久久久久久按摩| 99九九在线精品视频| 十八禁网站网址无遮挡| 日韩三级伦理在线观看| 日本黄色片子视频| 成年美女黄网站色视频大全免费 | 两个人免费观看高清视频| 最新中文字幕久久久久| 九色亚洲精品在线播放| 久久久久久久久久人人人人人人| 欧美性感艳星| 美女中出高潮动态图| 婷婷色综合大香蕉| 51国产日韩欧美| 不卡视频在线观看欧美| 自拍欧美九色日韩亚洲蝌蚪91| 国产爽快片一区二区三区| 男女啪啪激烈高潮av片| 日韩欧美精品免费久久| 中国美白少妇内射xxxbb| 五月玫瑰六月丁香| 中文字幕免费在线视频6| 这个男人来自地球电影免费观看 | 亚洲精品456在线播放app| 欧美精品一区二区大全| 久久午夜福利片| 人人妻人人爽人人添夜夜欢视频| 久久久国产一区二区| 亚洲精品乱码久久久v下载方式| 在线观看免费日韩欧美大片 | 国产熟女欧美一区二区| 综合色丁香网| 午夜免费观看性视频| 国产伦精品一区二区三区视频9| 精品一区二区三区视频在线| 成年人免费黄色播放视频| 久久久久人妻精品一区果冻| 国产日韩欧美在线精品| 亚洲欧洲国产日韩| 另类精品久久| 欧美xxⅹ黑人| 大片免费播放器 马上看| 在线精品无人区一区二区三| 午夜久久久在线观看| 欧美精品高潮呻吟av久久| 国产极品粉嫩免费观看在线 | 日韩在线高清观看一区二区三区| 母亲3免费完整高清在线观看 | a级片在线免费高清观看视频| 亚洲国产欧美在线一区| 内地一区二区视频在线| 国产黄色免费在线视频| 国产老妇伦熟女老妇高清| 国产成人精品在线电影| 久久青草综合色| 国产午夜精品一二区理论片| 在线观看三级黄色| 亚洲少妇的诱惑av| 久久久久精品久久久久真实原创| 97精品久久久久久久久久精品| 国产免费福利视频在线观看| 国模一区二区三区四区视频| 黄色配什么色好看| 黑人欧美特级aaaaaa片| 久久毛片免费看一区二区三区| 777米奇影视久久| 国产视频首页在线观看| 国产乱人偷精品视频| 欧美日韩视频精品一区| 日韩 亚洲 欧美在线| 免费观看在线日韩| 日韩欧美精品免费久久| 国产在线免费精品| 99九九线精品视频在线观看视频| 久久这里有精品视频免费| 亚州av有码| 精品亚洲乱码少妇综合久久| 欧美日韩在线观看h| 精品人妻熟女av久视频| 精品熟女少妇av免费看| 大片免费播放器 马上看| 国产在线免费精品| 久久99热这里只频精品6学生| 亚洲精华国产精华液的使用体验| 多毛熟女@视频| av免费在线看不卡| 国产精品蜜桃在线观看| 91国产中文字幕| 日韩一区二区三区影片| 交换朋友夫妻互换小说| av不卡在线播放| 全区人妻精品视频| 亚洲国产最新在线播放| 夫妻性生交免费视频一级片| 少妇猛男粗大的猛烈进出视频| 国产成人freesex在线| 18禁在线播放成人免费| 日韩欧美一区视频在线观看| 亚洲国产欧美在线一区| 午夜福利视频在线观看免费| 99热网站在线观看| 老熟女久久久| 国产欧美日韩综合在线一区二区| 嘟嘟电影网在线观看| 欧美精品一区二区免费开放| 青春草国产在线视频| freevideosex欧美| 最近的中文字幕免费完整| 中文天堂在线官网| 曰老女人黄片| 成人毛片a级毛片在线播放| 一本一本综合久久| 欧美日韩一区二区视频在线观看视频在线| 毛片一级片免费看久久久久| 精品酒店卫生间| 久久99热这里只频精品6学生| 在线精品无人区一区二区三| 国产成人av激情在线播放 | 视频在线观看一区二区三区| 水蜜桃什么品种好| 男的添女的下面高潮视频| av女优亚洲男人天堂| 亚洲av欧美aⅴ国产| 国产精品一区www在线观看| 久久久久久久久久久久大奶| 熟女人妻精品中文字幕| 久久热精品热| 午夜激情福利司机影院| 18在线观看网站| 亚洲精品成人av观看孕妇| 国产黄色免费在线视频| 少妇 在线观看| 亚洲美女黄色视频免费看| 免费看av在线观看网站| 国产一区二区三区av在线| 如何舔出高潮| 久久久久精品久久久久真实原创| 18禁观看日本| 亚洲综合色惰| 国产成人精品一,二区| 国产成人精品婷婷| 丝袜喷水一区| 亚洲精品色激情综合| 国产免费又黄又爽又色| 2018国产大陆天天弄谢| 欧美变态另类bdsm刘玥| 制服人妻中文乱码| 高清视频免费观看一区二区| 日本av手机在线免费观看| www.av在线官网国产| 热re99久久国产66热| 欧美日韩综合久久久久久| 女性被躁到高潮视频| 永久免费av网站大全| 狂野欧美白嫩少妇大欣赏| 久久99精品国语久久久| 啦啦啦啦在线视频资源| 亚洲天堂av无毛| 三上悠亚av全集在线观看| 极品人妻少妇av视频| 极品少妇高潮喷水抽搐| 五月玫瑰六月丁香| 美女中出高潮动态图| 狂野欧美激情性xxxx在线观看| 尾随美女入室| 校园人妻丝袜中文字幕| 岛国毛片在线播放| 99九九在线精品视频| 国产免费福利视频在线观看| 丝袜脚勾引网站| 久久99精品国语久久久| 国产精品蜜桃在线观看| www.色视频.com| 国产免费又黄又爽又色| 久久久精品94久久精品| 中国美白少妇内射xxxbb| 免费观看性生交大片5| 亚洲av二区三区四区| 午夜福利影视在线免费观看| 天堂中文最新版在线下载| 国产黄色视频一区二区在线观看| 久久99精品国语久久久| 欧美亚洲日本最大视频资源| 亚洲美女视频黄频| 免费少妇av软件| 欧美另类一区| 在线精品无人区一区二区三| 少妇猛男粗大的猛烈进出视频| 在线天堂最新版资源| a 毛片基地| 18在线观看网站| 丝袜在线中文字幕| 国产一区二区在线观看av| 国产一区有黄有色的免费视频| 2018国产大陆天天弄谢| 久久毛片免费看一区二区三区| 免费观看性生交大片5| 欧美日韩综合久久久久久| 曰老女人黄片| 日本与韩国留学比较| 蜜桃久久精品国产亚洲av| 成人国语在线视频| 成人午夜精彩视频在线观看| 黄色怎么调成土黄色| 九九在线视频观看精品| 午夜激情av网站| 婷婷色综合大香蕉| videos熟女内射| 男女国产视频网站| 人妻制服诱惑在线中文字幕| 成人无遮挡网站| 永久免费av网站大全|