姜 杰,張 濤,蔡曉仙,吳招才
(1.山東科技大學(xué)測繪與空間信息學(xué)院,山東 青島266000;2.自然資源部第二海洋研究所,浙江 杭州 310012;3.自然資源部海底科學(xué)重點實驗室,浙江 杭州 310012)
大地電磁法[1-2]采用天然電磁場作為場源,能反映不同深度范圍的電阻率信息,是海底深部結(jié)構(gòu)探測的主要手段。由于國外的技術(shù)封鎖,我國海洋大地電磁的探測近幾年才逐步興起,成為海底探測的一個主要熱點。國際上對海底火山的大地電磁探測與研究主要集中在洋中脊軸部火山區(qū),并取得了對其深部電阻率結(jié)構(gòu)的認(rèn)識:高電導(dǎo)率地殼與軟流層(5~50 Ω·m)由低阻巖石圈(50~100 Ω·m)連接,表明一小部分熔體從地幔的熔體源遷移到洋中脊軸部下的地殼巖漿室[3-4]。由于海洋板塊內(nèi)部火山(簡稱板內(nèi)海山)形成于較老的巖石圈上,其巖漿深部侵入和噴出的方式及相應(yīng)的結(jié)構(gòu)都與洋中脊軸部火山不同[5]。目前對遠(yuǎn)離洋中脊的板內(nèi)海山的電阻率結(jié)構(gòu)還未見報道,需要進(jìn)一步開展大地電磁探測工作以揭示其結(jié)構(gòu)特征。
蘇達(dá)海山位于馬爾庫斯-威克海山鏈,是典型的西太平洋板內(nèi)海山。依托于“大洋號”科考船2020年的西太平洋海底精細(xì)地球物理調(diào)查航次,自然資源部第二海洋研究所用我國最新研制的大地電磁接收機(jī)對蘇達(dá)海山進(jìn)行了首次大地電磁測量。本次調(diào)查揭示了蘇達(dá)海山的深部電阻率結(jié)構(gòu),對了解離軸板內(nèi)海山的構(gòu)建過程具有重要的意義。同時,由于蘇達(dá)海山是我國錳結(jié)核的礦區(qū),對其深部結(jié)構(gòu)的了解也具有一定的資源意義。
2020年9月—11月在西太平洋蘇達(dá)海山區(qū)域進(jìn)行了海洋大地電磁的測量,實驗區(qū)域及測點位置如 圖1 所示。西太平洋區(qū)域發(fā)育了一系列的板內(nèi)巖漿成因的海底山脈,包括馬爾庫斯-威克海山鏈、麥哲倫海山鏈和馬紹爾海山鏈等。蘇達(dá)海山位于馬爾庫斯-威克海山鏈北部,屬于遠(yuǎn)離洋中脊的板內(nèi)海山。
圖1 測區(qū)地形圖
近年來,我國大洋系列航次對蘇達(dá)海山的部分區(qū)域進(jìn)行了不規(guī)則的隨船多波束測量,并進(jìn)行了重力、磁力等地球物理調(diào)查和地質(zhì)采樣工作,但資料未見發(fā)表。KANEDA et al[5]利用OBS和多道地震資料揭示了位于蘇達(dá)海山西側(cè)的大型海山的地殼厚度可以超過15 km,巖漿侵入內(nèi)核寬度小于30~40 km,內(nèi)核與整體的體積比為15%~20%,表明大量的巖漿都通過噴發(fā)的形式釋放,形成了內(nèi)核周邊的火山碎屑沉積物。馬爾庫斯-威克海山區(qū)內(nèi)的海山沿NW—SE向散亂分布,整體形態(tài)并不呈現(xiàn)典型鏈狀,且沒有明確的年齡遞進(jìn)方向。前人使用40Ar/39Ar 定年法對馬爾庫斯-威克海山進(jìn)行了定年[6-9],馬爾庫斯-威克海山的年齡范圍在74~120 Ma之間,其中蘇達(dá)海山南側(cè)的拉蒙特海山的年齡為81.5~87.2 Ma。古地磁的研究也表明海山的形成過程沒有明確規(guī)律[10]。這與夏威夷-皇帝海山鏈成因為典型地幔柱存在很大的不同。NATLAND et al[11]認(rèn)為,包括巖石圈冷收縮、俯沖帶拖曳、海山負(fù)載等5種應(yīng)力的共同作用可能導(dǎo)致了古老巖石圈板塊內(nèi)的破裂,并使得地幔物質(zhì)發(fā)生熔融并沿裂隙形成火山活動。
本次實驗所用的儀器為OBEM-Ⅲ代海底大地電磁儀,是中國地質(zhì)大學(xué)(北京)與廣州海洋地質(zhì)調(diào)查局聯(lián)合研制的地球物理探測裝備,裝配2個電道、3個磁道,具有高可靠性、低噪聲、大動態(tài)范圍、低功耗、寬頻帶等優(yōu)點[12-16],儀器參數(shù)見表1。
表1 OBEM Ⅲ儀器參數(shù)
大地電磁數(shù)據(jù)處理主要是將采集到的時間域序列數(shù)據(jù)轉(zhuǎn)換為頻率域序列數(shù)據(jù),得到頻率域序列數(shù)據(jù)之間的傳遞函數(shù),最終獲得高質(zhì)量的阻抗張量、視電阻率和相位等參數(shù)的頻率響應(yīng)[17-19]。
海洋大地電磁數(shù)據(jù)處理的優(yōu)勢在于人為干擾所引起的電磁噪聲影響較小,對人為噪音的處理可以忽略[20];但由于海底環(huán)境的復(fù)雜性、海水的“低通濾波”作用、海水運動等產(chǎn)生的電磁噪音等[21],也增加了海洋大地電磁數(shù)據(jù)處理的難度。
本文按照圖2所示的數(shù)據(jù)處理流程,對海洋大地電磁接收機(jī)采集到的時間域原始數(shù)據(jù)進(jìn)行處理。首先為時間域數(shù)據(jù)處理,包括時間域數(shù)據(jù)預(yù)處理以及時間域到頻率域的轉(zhuǎn)換。時間域數(shù)據(jù)預(yù)處理是指時間域數(shù)據(jù)的挑選、噪音的去除以及轉(zhuǎn)換系數(shù)的校正;時間域到頻率域的轉(zhuǎn)換包括時間域數(shù)據(jù)時窗選擇(根據(jù)目標(biāo)頻率)和時間域到頻率域的轉(zhuǎn)換。其次為頻率域序列數(shù)據(jù)處理,包括選擇目標(biāo)頻點,以目標(biāo)頻點為中心完成時窗內(nèi)的數(shù)據(jù)疊加并進(jìn)行時窗內(nèi)的阻抗張量矩陣計算;進(jìn)一步通過最小二乘或魯棒(Robust)估計進(jìn)行不同時窗阻抗張量矩陣的疊加得到阻抗張量矩陣[22-23];通過阻抗張量可以得到初步的視電阻率和相位曲線;最后,根據(jù)姿態(tài)棒記錄的姿態(tài)數(shù)據(jù)對數(shù)據(jù)進(jìn)行方位校正[20],得到最終的視電阻率和相位數(shù)據(jù)。
圖2 海底大地電磁數(shù)據(jù)處理流程
圖3a為測點的原始時間域數(shù)據(jù),可看出Hx、Hy兩個通道存在兩種規(guī)律的周期性干擾信號,其分布特征分別為周期1 h、時長4 s以及周期6 h、時長2 s。這兩種規(guī)律噪音是由于儀器讀寫數(shù)據(jù)產(chǎn)生的。對于這種周期和時長固定的干擾信號,在時間序列文件中手動定位首個噪音的初始時間,根據(jù)其周期和時長即可對干擾信號進(jìn)行提取,進(jìn)而實現(xiàn)信噪分離。按照這種信噪分離的思路,去除Hx、Hy兩個通道中的規(guī)律噪音。隨后按照下述原則選擇數(shù)據(jù):(1)曲線應(yīng)具有很好的連續(xù)性,不連續(xù)的曲線,應(yīng)去除凌亂、斷檔、飛點等數(shù)據(jù);(2)Ex、Ey、Hx、Hy四個通道之間的數(shù)據(jù)曲線應(yīng)具有良好的一致性。最終選擇2020-10-01T18:58:07至2020-10-03T19:59:58區(qū)間(49 h)數(shù)據(jù)作為待處理數(shù)據(jù)。對選定數(shù)據(jù)進(jìn)行衰減系數(shù)、極距和增益校正等預(yù)處理后的結(jié)果如圖3b所示,數(shù)據(jù)曲線具有很好的連續(xù)性,各個通道間的曲線一致性良好。
圖3 時間域數(shù)據(jù)序列
采用SSMT2000方法對預(yù)處理后的數(shù)據(jù)進(jìn)行下一步處理。SSMT2000是加拿大鳳凰地球物理公司生產(chǎn)的商用大地電磁處理軟件,其時間域到頻率域的轉(zhuǎn)換采用傳統(tǒng)的傅里葉變換方法[24],不同窗口的疊加采用魯棒估算。魯棒估算方法是根據(jù)觀測誤差的剩余功率譜大小對數(shù)據(jù)進(jìn)行加權(quán),突出未被干擾的數(shù)據(jù),降低突變點數(shù)據(jù)的權(quán),使突變點對阻抗估算值的影響最小,從而明顯改善受電磁噪聲污染的單站大地電磁資料[22-23,25-26]。此次處理具體采用的魯棒估計形式為:比較電場和磁場的相干性,將相干性弱的去掉,去除各個通道的非相干噪音;權(quán)重方式為賦予磁場通道和電場通道之間一致性好的數(shù)據(jù)點更多的權(quán)重,可對數(shù)據(jù)進(jìn)行更好的處理,降低誤差棒的大小,使數(shù)據(jù)曲線更加光滑。
數(shù)據(jù)處理結(jié)果見圖4。海水是良導(dǎo)體,具備“低通濾波器”的效應(yīng),能夠壓制測量數(shù)據(jù)的高頻部分信號。假設(shè)海水的電阻率為0.3 Ω·m,根據(jù)趨膚定律,當(dāng)站位水深為1 434 m,小于30 s周期的數(shù)據(jù)會被海水衰減掉,因此在后續(xù)數(shù)據(jù)處理和反演中去除小于30 s周期的短周期數(shù)據(jù)。
圖4 SSMT2000方法處理結(jié)果
海洋大地電磁測量在海底進(jìn)行,儀器在海面釋放,以自由姿態(tài)到達(dá)海底,因此測量方位是任意的。需要通過姿態(tài)棒記錄的方位角θ對數(shù)據(jù)進(jìn)行方位校正。按照公式(1),將阻抗張量Z旋轉(zhuǎn)到磁北方向,旋轉(zhuǎn)后的數(shù)據(jù)如圖5所示。
圖5 方位校正后的視電阻率(a)和相位(b)數(shù)據(jù)
(1)
式中:Zxx、Zxy、Zyx、Zyy為任意方位的阻抗張量;Zx′x′、Zx′y′、Zy′x′、Zy′y′為旋轉(zhuǎn)后磁北方向的阻抗張量;θ為任意方位與磁北方向之間的夾角。
大地電磁阻抗張量具有旋轉(zhuǎn)不變的特性。目前研究表明,阻抗張量旋轉(zhuǎn)不變量不受電性主軸的旋轉(zhuǎn)和地形畸變等的影響,能夠很好地解釋三維構(gòu)造特征[27]。本文采用兩種旋轉(zhuǎn)不變量,第一種為BERDICHEVSKY et al[28]提出的方案[公式(2)],通過平均有效電阻率(等效于由阻抗張量行列式導(dǎo)出的視電阻率),從畸變數(shù)據(jù)估計平均一維剖面模型,它能有效地協(xié)調(diào)復(fù)雜地點結(jié)構(gòu)的方向特征。第二種為RUNG-ARUNWAN et al[29]提出的估算平均一維剖面模型的方法[公式(3)],該方法類似于使用BERDICHEVSKY et al[28]的平均值方法,但將旋轉(zhuǎn)不變量定義為阻抗張量平方元素的和(ssq阻抗)。如圖6所示,Zdet、Zssq的數(shù)據(jù)位于XY方向和YX方向數(shù)據(jù)范圍內(nèi),并且兩者基本一致,表明Zdet、Zssq數(shù)據(jù)是兩個方向數(shù)據(jù)的綜合。
圖6 XY方向、YX方向、Zdet和Zssq的視電阻率(a)和相位(b)比較
(2)
(3)
在一維反演前,需要先檢驗數(shù)據(jù)是否適合一維反演。采用Rho+算法對數(shù)據(jù)進(jìn)行檢驗,該算法是對一維結(jié)構(gòu)假設(shè)下MT響應(yīng)函數(shù)一致性的測試[30],并通過設(shè)置視電阻率和相位的上下界限清除溢界值[31]。
對測點的XY方向、YX方向、Zdet、Zssq數(shù)據(jù)進(jìn)行Rho+算法一致性測試,其結(jié)果如圖7所示。除去高頻段干擾后,4組視電阻率數(shù)據(jù)均與Rho+算法一致。相位數(shù)據(jù)方面,YX方向的一致性最好(圖7b);Zdet和Zssq在80~140 s區(qū)間(圖7c和7d),數(shù)據(jù)溢出較嚴(yán)重,其余數(shù)據(jù)一致性較好;XY方向(圖7a)整體溢出較嚴(yán)重,與Rho+算法的一致性差。因此,采用對一維結(jié)構(gòu)響應(yīng)最好的YX方向數(shù)據(jù)進(jìn)行一維反演。
圖7 Rho+算法測試結(jié)果
反演采用OCCAM1DCSEM算法[31-34],針對反演的多解性,該算法通過正則化反演[30,34]確定與給定數(shù)據(jù)集兼容的電導(dǎo)率模型。在模型中,對應(yīng)數(shù)據(jù)約束良好的特征被保留,而不受數(shù)據(jù)約束的特征將在模型中被平滑或者去掉[35-36]。
一維反演采用的模型如圖8a所示:模型為各向同性;采用右手坐標(biāo)系,Z軸垂直向下;從上到下分別為空氣層、水層、海底以下的自由電阻率層??諝鈱訛楣潭▽樱o定層寬10 km,電阻率1012Ω·m;水層為30層,層寬均勻分布,共1 434 m,每層海水電阻率值給定實際深度的電阻率值,避免海水層電阻率突變造成反演結(jié)果畸變[35],使用的海水電阻率值為太平洋圣選戈海槽區(qū)域的海水電阻率數(shù)據(jù),為 OCCAM1DCSEM 算法程序包中自帶的海水電阻率文件;海底以下的自由電阻率層為50層,從上到下層寬按對數(shù)分布,共100 km,允許數(shù)據(jù)的自由反演。對測點的YX方向數(shù)據(jù)進(jìn)行一維反演的結(jié)果如圖8b所示。海底下約2~3 km處存在電阻率突尖,可能反映出蘇達(dá)海山表層沉積物較薄,電阻率較大的巖層接近海底面;海底向下15 km處電阻率最小,表明電阻率較大巖層下存在電阻率較小層;15 km以下電阻率逐漸變大,在 40~60 km處形成一電阻率高值,隨后電阻率隨著深度的增加而逐漸減小。這種變化特征符合地殼、地幔的電阻率分布特征:15 km以下的電阻率升高可以解釋為地殼輝長巖向高電阻率的地幔橄欖巖的轉(zhuǎn)換;60 km以下逐步降低的電阻率反映了橄欖巖電阻率隨深度和溫度的增加而逐漸變小。
圖8 100 km反演模型及反演結(jié)果
實測數(shù)據(jù)周期最長可達(dá)到10 000 s(圖5),對應(yīng)趨膚定律計算的深度超過100 km。為進(jìn)一步研究該地區(qū)的深部構(gòu)造與巖石圈厚度,我們在保持其它參數(shù)不變的情況下,將反演模型加深到200 km(圖9a)。結(jié)合圖11可看出,在90 km以上的深度,2個模型的反演結(jié)果整體較為一致;90~110 km,電導(dǎo)率與深度之間的曲線斜率發(fā)生了顯著轉(zhuǎn)折,之后趨于穩(wěn)定,我們由此推測此處為巖石圈-軟流圈邊界區(qū)域,即蘇達(dá)海山的巖石圈厚度約為100 km。
圖9 200 km反演模型的反演結(jié)果
為進(jìn)一步確定海山的電阻率結(jié)構(gòu),本文結(jié)合YX方向反演結(jié)果與相關(guān)文獻(xiàn)[5,37-39],建立理論電阻率結(jié)構(gòu)模型,正演合成其數(shù)據(jù)后進(jìn)行反演,并與YX方向反演結(jié)果對比,找出更相符的地下電阻率結(jié)構(gòu)模型。
結(jié)合西北太平洋馬爾庫斯海山鏈的地震速度結(jié)構(gòu)特征[5]、西南太平洋Tūranganui海山的電性結(jié)構(gòu)特征[37],通過多次試驗,得到最佳模型如下:海底至 5 km 深度為電阻率100 Ω·m的玄武巖層;5~13 km為火山碎屑巖層,電阻率為20 Ω·m; 13~18 km為火山碎屑巖層,電阻率為3 Ω·m; 18~23 km為輝長巖層,電阻率為300 Ω·m;結(jié)合MATSUNO et al[38-39]對馬里亞納區(qū)域一維深部地下電性結(jié)構(gòu)的研究,給定23~100 km為橄欖巖層,電阻率為 1 000 Ω·m;100 km以下電阻率設(shè)為1 Ω·m(圖10a)。模型的反演結(jié)果與觀測值較為符合,可以復(fù)現(xiàn)海底下的電阻率突尖、15 km 左右的電阻率最小值、40~60 km間形成的電阻率高值(圖10b和圖11)。因此,3.5 km 厚的玄武巖層、13 km厚的火山碎屑巖以及5 km厚的輝長巖,可以解釋蘇達(dá)海山的觀測電阻率。在15 km處電阻率值大于實測數(shù)據(jù)反演值以及在40~60 km處電阻率值小于實測數(shù)據(jù)反演值,說明實測數(shù)據(jù)反演結(jié)果可能存在過度擬合的情況。
圖10 最佳模型(a)與其對應(yīng)的反演響應(yīng)曲線(b)
圖11 反演結(jié)果對比圖
蘇達(dá)海山的平均水深為1 500 m,周邊深海盆為 5 000 m,假定周邊深海盆為洋殼,平均厚度為6 km,洋殼密度為2.8 g/cm3,地幔密度為3.3 g/cm3?;诎锏貧ぞ饽P?,蘇達(dá)海山的地殼厚度約為 22 km,這與本文用電磁數(shù)據(jù)推斷的21.5 km的地殼厚度較為一致。蘇達(dá)海山較厚的火山碎屑巖說明其噴發(fā)部分占了主體,這可能與蘇達(dá)海山噴發(fā)時的水深較淺、噴發(fā)壓力較小有關(guān),也可能是由于蘇達(dá)海山的巖漿揮發(fā)性比較強(qiáng),更易于噴發(fā)。同時,蘇達(dá)海山的結(jié)構(gòu)與馬爾庫斯-威克海山鏈上小規(guī)模火山的地震推斷的地殼結(jié)構(gòu)類似[5],表明這種小型海山的形成可能都具有以噴發(fā)性為主、侵入性較弱的模式。
本文首次對板內(nèi)海底火山——蘇達(dá)海山的海洋大地電磁調(diào)查數(shù)據(jù)進(jìn)行了處理,對旋轉(zhuǎn)后的數(shù)據(jù)以及旋轉(zhuǎn)不變量進(jìn)行了一維反演,進(jìn)一步通過一維正演來推斷蘇達(dá)海山的巖石圈結(jié)構(gòu)。蘇達(dá)海山的電阻率模型表明:海山表層覆蓋著薄玄武巖層,厚度在 3.5 km 左右;薄玄武巖層下是低電阻率的火山碎屑巖,厚度為13 km左右;其下為5 km厚的輝長巖層。通過大地電磁數(shù)據(jù)推斷的蘇達(dá)海山地殼厚度約為 21.5 km,與艾里地殼均衡模型推斷的22 km的地殼厚度相符。蘇達(dá)海山較厚的火山碎屑巖層表明其噴發(fā)部分為主體,與地震推斷的馬爾庫斯-威克海山鏈上小規(guī)?;鹕降牡貧そY(jié)構(gòu)類似。
致謝感謝自然資源部第二海洋研究所2020年西太平洋海底精細(xì)地球物理調(diào)查航次的船長楊江志在數(shù)據(jù)采集過程中提供的幫助;感謝中國地質(zhì)大學(xué)(北京)陳凱博士及其團(tuán)隊在儀器及數(shù)據(jù)預(yù)處理方面的指導(dǎo);感謝自然資源部第二海洋研究所秦林江博士、南方科技大學(xué)趙云生博士在數(shù)據(jù)處理及正、反演方面提供的幫助;感謝山東科技大學(xué)趙俐紅博士提供的浪潮TS10K集群的算力支持;感謝山東科技大學(xué)王飛燕博士在反演算法上的指導(dǎo)。