朱益民 余 騰 王 鑫 張 鵬

1 中國地質大學(北京)土地科學技術學院，北京市學院路29號，100083 2 宿遷學院建筑工程學院，江蘇省宿遷市黃河南路399號，223800 3 宿遷學院地質環(huán)境與工程健康監(jiān)測重點實驗室, 江蘇省宿遷市黃河南路399號，223800 4 應急管理部國家自然災害防治研究院，北京市安寧莊路1號，100085 5 南京工業(yè)大學交通運輸工程學院，南京市浦珠南路30號，211816

通過分析重力異常的梯度變化能夠確定異常體的邊緣與規(guī)模，進而研究地殼內密度異常變化及其與地震活動性的聯(lián)系[1]。部分學者通過垂線偏差反演構造應力場，探討應力狀態(tài)與地質體運動的基本特征[2-3]。國內外學者利用重力梯度數(shù)據(jù)識別地質體邊界、劃分構造單元[4-5]，借助重力總水平梯度和垂線偏差研究區(qū)域構造應力場[6-7]，基于重力場模型計算不同深度下的地殼密度異常與構造應力場[8]。

針對郯廬帶的構造特征研究多基于GPS位移分析[9]、地震反射剖面[10]、地應力測量[11]、地殼速度結構及泊松比[12]、地震重定位與震源機制解以及地震活動性[13-15]等數(shù)據(jù)，針對郯廬帶及鄰區(qū)的重力異常數(shù)據(jù)處理多采用小波多尺度分解方法[16-17]以及Tilt導數(shù)法[18]，鮮有基于重力全張量梯度數(shù)據(jù)進行的郯廬帶中南段構造特征研究?；谥亓鲆浑A導數(shù)定義的傳統(tǒng)方法僅適用于淺層構造和礦產勘探目標界定[19]，且會因為同時存在正負異常而產生錯誤邊界；而重力全張量梯度數(shù)據(jù)已在南海斷裂圈定[20]、四川盆地地質構造水平位置界定[21]等研究與應用中展現(xiàn)出較好的適用性。因此，基于重力全張量研究郯廬帶中南段地質體交接關系與構造單元劃分具有一定的實用意義。

本文通過重力全張量數(shù)據(jù)重新定義水平導數(shù)THDR和解析信號振幅ASM，并提出改進的傾斜導數(shù)NTDM。經過模型驗證，該方法可均衡識別不同埋深下的地質體邊界，具有較好的穩(wěn)定性與收斂性。同時，利用垂線偏差數(shù)據(jù)反演研究區(qū)構造應力場與殼幔密度差發(fā)現(xiàn)，二者在空間上具有較強的相關性。研究區(qū)最大主應力方向總體呈EW向，與前人研究結果一致。

1 研究區(qū)概況

1.1 郯廬帶中南段地形及構造特征

郯廬帶是我國東部規(guī)模最大的NNE走向活動斷裂，其中新沂以北至渤海為郯廬帶中段，新沂以南為郯廬帶南段。郯廬帶中南段貫穿本文研究區(qū)(32.7°～35.2°N、117.0°～119.5°E)，研究區(qū)內以淺源地震為主，震群較為集中。區(qū)域內多條斜列的NWW、EW向斷裂與郯廬帶中南段交接，組成復雜的斷陷帶(圖1)。第四紀以來，在西太平洋板塊弧后擴張與印度板塊向北碰撞的共同作用下，郯廬帶表現(xiàn)為逆沖兼右旋走滑的活動性質。

圖1 研究區(qū)地形與主要構造特征[22]Fig.1 Topography and main structural features of the study area[22]

1.2 郯廬帶中南段布格重力場、莫霍面埋深及地殼厚度特征

本文收集由BGI官方網(wǎng)站發(fā)布的2′×2′格網(wǎng)大小的WGM2012模型計算布格重力異常(圖2(a))，此外還收集了Sandwell V31的1′×1′垂線偏差數(shù)據(jù)，用于反演研究區(qū)構造應力場、殼幔密度差和重力總水平梯度。地殼厚度與莫霍面埋深數(shù)據(jù)來自Crust1.0(圖2(b))。

郯廬帶是一條重要的地球物理分隔帶。由圖2(a)可見，NNE走向的重力異常分隔帶呈現(xiàn)西低東高的特征，反映出華北板塊和揚子板塊基底性質的不同，深部構造復雜；西側的華北板塊表現(xiàn)為數(shù)值偏低的布格異常區(qū)，變化趨勢為EW向。F10(淮陰-響水斷裂)和F11(洪澤-溝墩斷裂)共同形成弧形串珠狀低重力異常區(qū)，分隔下?lián)P子板塊。

郯廬帶也是一條NNE走向的莫霍面突變帶。由圖2(b)可見，總體上東側莫霍面埋深較淺，西側莫霍面埋深較深。莫霍面埋深由北至南逐漸增大，局部還存在較為明顯的分段性。

2 數(shù)據(jù)與方法

2.1 重力場方向導數(shù)法與全張量梯度

對布格重力異常網(wǎng)格數(shù)據(jù)進行x、y、z三個坐標軸方向的求導，得到重力矢量g在笛卡爾坐標系中沿坐標軸方向的分量gx、gy、gz，可反映淺層地質體的交接關系，初步揭示區(qū)域構造格局。重力全張量梯度為重力向量沿坐標軸方向分量的一階導數(shù)，其在三維直角坐標系中可表示為：

(1)

總水平導數(shù)THDR表達式為[1]：

(2)

解析信號振幅ASM表達式為[23]：

(3)

圖2 研究區(qū)布格重力異常和地殼厚度與莫霍面埋深Fig.2 Bouguer gravity anomaly and crustal thickness and Moho face buried depth of the study area

由于式(2)、(3)是基于一階導數(shù)定義的邊界識別方法，在圈定較深地質構造的邊界時存在較大誤差，因此為使淺層與深部地質體邊界圈定達到均衡狀態(tài)，本文利用重力全張量數(shù)據(jù)重新定義式(2)、(3)：

(4)

(5)

對THDRM的總水平導數(shù)與ASMM的垂向導數(shù)進行融合處理，定義基于重力全張量的傾斜角法NTDM：

(6)

THDRM和ASMM的局部極值和最值處可以界定地質體邊界，傾斜導數(shù)NTDM在±π/2的振幅范圍內隨傾斜角發(fā)生顯著變化。零交叉點接近異常體結構邊緣，局部極值能有效圈定其邊界，同時也能判定地質體在平面位置的延伸方向。經過3種不同實驗環(huán)境下的模型驗證后認為，該方法在均衡探測深淺斷裂的同時，對邊界識別也具有較好的穩(wěn)定性與收斂性，能清晰區(qū)分異常體區(qū)域與邊界，便于邊界解譯。圖3為式(2)～(5)的數(shù)值模擬對比圖，圖中黑色虛線為地質體水平投影所對應的邊界。

2.2 構造應力場與殼幔密度差

構造應力是在正壓力(靜巖應力)上附加的一種應力，文獻[2]在Turcotte公式的基礎上，基于垂線偏差推導其與構造應力場的轉換公式?；诖?，本文推導出如下公式：

(7)

(8)

(9)

由垂線偏差定義的最大主應力方向為：

(10)

圖3 數(shù)值模擬對比Fig.3 Comparison of numerical simulation

3 地質體交接關系識別與分析

3.1 重力異常方向導數(shù)

對布格重力數(shù)據(jù)進行方向導數(shù)處理，可突出其局部梯度變化信息，位場解析延拓也可使目標場源特征更加明顯。

圖4(a)～(c)為經方向導數(shù)處理后的重力異常與研究區(qū)斷裂位置的對應情況。由圖4(a)可見，重力場x方向導數(shù)gx表現(xiàn)出近NS向構造邊界，橫向上具有顯著差異，西側為正梯度帶，東側為負梯度帶，變化幅度為0.4×10-3mGal/km。由圖4(b)可見，重力場y方向導數(shù)gy表現(xiàn)出近EW向地質體邊界，縱向上的正負異常梯度交替出現(xiàn)。以F6(無錫-宿遷斷裂)為界，研究區(qū)中部以低異常梯度為主，分隔南北正異常梯度帶。由圖4(c)可見，重力場z方向導數(shù)gz表現(xiàn)出垂向變化特征，與區(qū)域布格重力異常的相似性較高，西北側的華北板塊轉變?yōu)槠骄彽牡彤惓Ｌ荻?，無明顯變化。淺層地殼的剛體性質較為明顯，郯廬帶中南段處在北寬南窄的正異常梯度帶上。

(a)、(b)中的斷裂分別為篩選后的NS向、EW向斷裂圖4 郯廬帶中南段布格重力場方向導數(shù)與向上延拓情況Fig.4 Directional dervatives and upward continuation of Bouguer gravity anomaly in the middle-southern segment of Tan-Lu fault zone

3.2 基于全張量梯度的NTDM邊界識別方法

3.2.1 模型驗證

為驗證NTDM方法的可行性，建立4個不同體積與底部埋深的棱柱體模型(圖5)，對THDRM、ASMM、NTDM以及重力全張量gxx、gxy、gxz、gyy、gyz、gzz進行模擬計算與驗證。 模型基本參數(shù)如表1所示。

表1 棱柱體模型參數(shù)

圖6為模型的重力梯度異常(圖中紅、橙、黑、綠實線位置表示模型M1、M2、M3、M4在x軸上的投影，圖7、8同理)，其中gxx、gyy、gzz可識別特定方向的邊界。埋深較深的模型M2、M4的異常振幅與周圍存在明顯差異，能界定異常體的規(guī)模。gxx對邊界位置識別較為發(fā)散，gxy可反映異常體角點信息,gxz、gyz可對應邊界信息。與THDRM相比，ASMM對不同埋深的異常振幅大小的均衡探測效果較好，二者均對埋深較深的模型異常體規(guī)模具有較好的揭示效果。對于埋深最淺的模型M1，THDRM未能對其進行有效界定。剖面圖能更直觀地展示異常振幅的差異和不同方法的分辨率特征。如圖6(g)中黑色虛線和紅色箭頭所示，THDRM的局部極小值可界定異常體邊界；而ASMM主要以最大值圈定邊界，均衡圈定的能力較好，且ASMM邊界識別的收斂性強于THDRM。由圖6(i)可見，NTDM具有邊界增強作用，邊界識別的分辨率高于圖6(a)～(h)，NTDM剖面圖的極大值可對應異常體邊界。

圖5 模型剖面圖、三維圖、重力異常Fig.5 Section, 3D graph, gravity anomaly of model

圖6 模型重力異常的全張量及經全張量組合后的邊界識別結果Fig.6 Full tensors of gravity anomaly and the boundary identification results after the combination of the full tensors

本文在重力梯度分量中加入5%的高斯隨機噪聲，以驗證NTDM的抗噪能力(圖7)。從圖7(g)～(i)可以看出，NTDM具有較強的抑噪能力，剖面圖中虛線位置與模型的水平投影位置對應較好，NTDM在均衡探測異常體邊界與抗噪能力上均最優(yōu)。但由于NTDM的表達式中包含三階導數(shù)，會放大異常體周邊的噪聲影響，因此剖面線的異常擾動較為明顯。

通過向上延拓的方法對數(shù)據(jù)進行濾波處理，延拓后的數(shù)據(jù)可抑制較淺的異常體。由圖8(g)～(i)所示，模型M1底部埋深最淺，THDRM、ASMM、NTDM對其邊界信息的識別總體較弱，其中NTDM的識別效果優(yōu)于THDRM、ASMM，這是因為THDRM和ASMM的識別易受到周邊異常體的干擾，結果較為發(fā)散。3個邊界識別方法對埋深較大異常體的規(guī)模圈定效果較好，NTDM在地質體周邊的正負異常對比使其具有更穩(wěn)定的邊界識別能力。

圖7 加入5%高斯噪聲重力異常的全張量及經全張量組合后的邊界識別結果Fig.7 Full tensors of gravity anomalies with 5% Gaussian noise and the boundary identification results after the combination of the full tensors

圖8 加入5%高斯噪聲并向上延拓1 km后重力異常的全張量及經全張量組合后的邊界識別結果Fig.8 Full tensors of gravity anomaly with 5% Gaussian noise and 1 km upward continuation, and the boundary identification results after the combination of the full tensors

經3種實驗環(huán)境驗證，THDRM、ASMM、NTDM的識別效果總體上優(yōu)于單個重力全張量梯度。THDRM主要以局部極小值確定邊界，ASMM、NTDM主要以局部極大值和最大值圈定邊界。3種方法均具有較好的抗噪性與收斂性，不會產生錯誤邊界，其中NTDM的邊界識別結果更清晰。

3.2.2 實際數(shù)據(jù)應用

為驗證本文方法在實際重力場數(shù)據(jù)中的應用效果，利用上述方法對郯廬帶中南段及鄰區(qū)布格重力異常數(shù)據(jù)進行處理，確定地質體的交接關系與構造邊界(圖9)。由圖可見，重力梯度gxx中央的低異常區(qū)與郯廬帶中南段對應較好，但異常區(qū)規(guī)模更大、東西向更寬；gyy、gyz能識別EW向斷裂，同時斜切的EW和NW向斷裂也使得郯廬帶中南段具有分段性；受到F6(無錫-宿遷斷裂)、F7(銅石-甘霖斷裂)的阻隔作用，gzz可反映郯廬帶中南段深部構造的不連通性。

相較于ASMM，THDRM邊界識別的分辨率更高，南部的串珠狀弧形分隔帶有利于揭示邊界信息，但ASMM在圈定地質體規(guī)模上具有一定的優(yōu)勢。新方法NTDM在均衡深部和淺部構造特征上更具優(yōu)勢，反映出的信息與已有成果對應較好，還可以展示周邊更豐富的細節(jié)信息，不受噪聲干擾，可為地質構造單元的劃分提供依據(jù)。

圖9(i)邊界識別結果顯示，以郯廬帶為界，研究區(qū)中部的地質體東西兩側受到郯廬帶的阻擋作用，延伸角度的弧度多為負值，方向呈NW、EW向，與郯廬帶交接；研究區(qū)南部地質體延伸方向發(fā)生明顯扭轉，在水平向與垂直向上均發(fā)生變化。F11(洪澤-溝墩斷裂)的東南側存在一個由EW向轉為NE向的異常地質體邊界，可作為構造單元劃分的依據(jù)。

圖9 郯廬帶中南段重力異常邊界識別結果Fig.9 Identification results of gravity anomaly boundary in the middle-southern Tan-Lu fault zone

綜合圖9(g)～(i)的識別結果，結合前人研究資料[24-26]，對郯廬帶中南段及其周邊地區(qū)進行劃分，共得到14個構造單元(圖10)。由圖可見，區(qū)域Ⅰ郯廬帶作為二級構造單元，分隔東西不同的地塊，區(qū)域Ⅱ可能為膠遼地塊南沿部分。F8(海泗斷裂)具有北延跡象，其北界與F9(邵店-桑墟斷裂)和郯廬帶交接位置可作為區(qū)域Ⅲ蘇北膠南斷塊的劃分依據(jù)。 F10(淮陰-響水口斷裂)不僅是華北陸塊區(qū)和秦祁昆造山系的分界線，同時也是魯西陸塊和大別-蘇魯造山帶的分界線。根據(jù)邊界識別結果可知，北延的次級斷裂可劃分區(qū)域Ⅳ連云港-泗洪斷塊。區(qū)域Ⅴ鹽阜坳陷、區(qū)域Ⅵ建湖隆起和區(qū)域Ⅶ東臺坳陷可由研究區(qū)東南方位的弧形異常邊界進行劃分，且異常邊界與F11(洪澤-溝墩斷裂)的對應較好。研究區(qū)西側呈現(xiàn)出弧形的NW向異常帶(方位角弧度值為負)，分隔了區(qū)域Ⅷ魯中隆起和區(qū)域Ⅸ魯西斷塊，以及區(qū)域Ⅹ邳縣-睢寧斷塊和區(qū)域Ⅺ徐淮坳褶帶。區(qū)域ⅩⅢ靈璧-泗縣凹陷屬于區(qū)域Ⅻ淮北斷陷帶，可由NW向陡變?yōu)镋W向的異常邊界進行劃分。ⅩⅢ靈璧-泗縣凹陷和區(qū)域ⅩⅣ蚌埠隆起可由EW向的弧形異常邊界進行劃分。

圖10 郯廬帶中南段構造單元劃分Fig.10 Division of tectonic units in the middle-southern segment of Tan-Lu fault zone

4 基于垂線偏差的構造應力場反演

根據(jù)§2.2理論基礎，借助重力垂線偏差數(shù)據(jù)反演研究區(qū)構造應力場、殼幔密度差和重力總水平梯度(圖11)。

藍色實線為沿郯廬帶剖面線作垂直于剖面線左右各30 km的剖面(圖中等間隔的白色實線)，繪制成上圖中藍色曲線圖11 構造應力場、殼幔密度差、重力總水平梯度Fig.11 Tectonic stress field, crust-mantle density difference, total horizontal gradient of gravity

由圖11可見，構造應力場、殼幔密度差和重力總水平梯度東西向差異明顯，水平方向上很不均勻，由圖中曲線可以看出，郯廬帶是一條構造應力差異變化分隔帶。在構造應力變化明顯的區(qū)域，殼幔密度和重力總水平梯度變化也較大。斷裂帶兩側構造應力等值線密集，多呈規(guī)則的環(huán)形，地震大多發(fā)生在2個環(huán)形等值線的交接處。同時，環(huán)形構造應力等值線或許可以作為局部構造單元的劃分依據(jù)，如多個環(huán)形等值線與區(qū)域ⅩⅢ靈璧-泗縣凹陷、區(qū)域ⅩⅣ蚌埠隆起的對應較好，研究區(qū)東南側的弧形變化區(qū)與F11(洪澤-溝墩斷裂)的延伸方向也具有相似性。

圖12為郯廬帶中南段構造應力場矢量圖，由圖可見，最大主應力方向整體以EW向為主。中段東北角主應力方向由NW轉為SW，該區(qū)域構造應力值較大、等值線密集，受到明顯的擠壓作用，最大主應力方向為NEE-EW；南段以F10(淮陰-響水口斷裂)為分界線，東南角由NW向一直延伸至蚌埠一帶，淮安-洪澤一帶由NW和SW向轉為EW向并延伸至郯廬帶及以西地區(qū)，最大主應力變化區(qū)也與高構造應力值對應較好。郯廬帶中南段應力場整體方向呈EW向，與受到太平洋板塊向歐亞大陸俯沖和印度洋板塊共同作用形成的中國大陸東部的應力場一致。表2為前人對最大主應力方向的研究成果。

實線為郯廬帶中段與南段分界線圖12 郯廬帶中南段最大主應力方向Fig.12 The direction of the maximum stress in middle-southern segment of Tan-Lu fault zone

表2 最大主應力方向對比

5 結 語

本文根據(jù)重力全張量定義的總水平導數(shù)與解析信號振幅，提出新的傾斜導數(shù)邊界識別方法NTDM，并將其應用到郯廬帶中南段布格重力異常數(shù)據(jù)的處理與解釋中。借助垂線偏差可反演得到郯廬帶中南段的構造應力場特征，得出以下結論：

1)當正負異常同時出現(xiàn)并含有噪聲時，NTDM能表現(xiàn)出良好的抑噪性與收斂性，能清晰識別異常體邊界、均衡不同深部的異常體振幅。相比于單一的全張量數(shù)據(jù)，THDRM、ASMM在界定局部邊界時具有一定的優(yōu)勢，NTDM能更準確地圈定邊界。

2)依據(jù)邊界識別結果并結合相關資料，以區(qū)域Ⅰ郯廬帶中南段為界，大致可將該區(qū)域劃分為區(qū)域Ⅱ膠遼地體、區(qū)域Ⅲ蘇北膠南斷塊、區(qū)域Ⅳ連云港-泗洪斷塊、區(qū)域Ⅴ鹽阜坳陷、區(qū)域Ⅵ建湖隆起、區(qū)域Ⅶ東臺坳陷、區(qū)域Ⅷ魯中隆起、區(qū)域Ⅸ魯西斷塊、區(qū)域Ⅹ邳縣-睢寧斷塊、區(qū)域Ⅺ徐淮坳褶帶、區(qū)域Ⅻ淮北斷陷帶、區(qū)域ⅩⅢ靈璧-泗縣凹陷、區(qū)域ⅩⅣ蚌埠隆起等14個構造單元。

3)構造應力場、殼幔密度差與重力總水平梯度均存在明顯的以郯廬帶中南段為界的水平差異，三者與區(qū)域重力異常的相關性較大，構造應力值范圍為-28～-12 Mpa，東西向最大差異約為16 Mpa，殼幔密度差水平向最大差異約為0.2 g/cm3。構造應力等值線密集區(qū)的構造運動強烈，地震頻發(fā)。

4)研究區(qū)最大主應力方向總體呈EW向，其中郯廬帶中段為NEE-EW，郯廬帶南段為EW-NWW，與中國大陸東部應力場保持較高的一致性。

致謝：應急管理部國家自然災害防治研究院和宿遷市地震服務中心提供相關資料；開源Python庫GEOIST(https:∥cea2020.gitee.io/geoistdoc/)對本文模型驗證提供支持；文中部分圖件由GMT繪制，并參考有關實例，在此一并表示感謝。