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