趙文舉，趙荔，楊戰(zhàn)軍，陶德強

(中國石油東方地球物理勘探有限責任公司 綜合物化探處，河北 涿州 072751)

0 引言

隨著計算機的廣泛應(yīng)用，重磁異常分離方法早已經(jīng)由手動發(fā)展為自動，且方法眾多。廣義來說，改變量綱的位場變換方法也屬于位場分離，這里我們討論的主要是不改變量綱的分離方法。主要有波數(shù)域[1]和空間域兩種，波數(shù)域方法如溫納濾波[2]、匹配濾波[3-4]、小波逼近[5-6]、正則化濾波[7]、巴特沃斯濾波等，空間域方法如最小二乘法[8]、窗口滑動平均法、趨勢分析法、圓周平均法、插值切割法等眾多方法，不勝枚舉。

插值切割法[9]由文百紅等人提出，該方法采用圓周上的4個正交點值以及半徑來構(gòu)建切割算子，由于其引入的半二階差分量與負二階水平導數(shù)成正比，因此，其對位場的不同非線性部分有不同的作用，從而提高了對不同特征重磁異常的分辨率。該方法近年來在位場數(shù)據(jù)的去噪、半定量分離等方面被應(yīng)用，甚至用于估計場源體的深度和三維成像[10-17]。

筆者首先介紹插值切割的原理。其次，就該方法在模型和實際應(yīng)用時所碰到的問題進行了部分改進，并對改進前后的效果進行了對比分析。最后，給出了該方法在深層構(gòu)造異常提取和火成巖三維成像方面的應(yīng)用實例，取得了較好的應(yīng)用效果。

1 方法原理

1.1 基本原理

圖1中g(shù)0是某點的重力異常值，g1，g2，g3和g4是半徑為r的周圍4個點的異常值，用下列公式計算g0點第一次切割的區(qū)域場：

圖1 插值切割法示意Fig.1 Interpolating cut method

(1)

(2)

α=f(g0,g1,g2,g3,g4),

(3)

其中：α是與g0，g1，g2，g3，g4有關(guān)的函數(shù)。

α=b+c，

δgx=g1-g3, δgy=g2-g4。

由此得出全區(qū)的第一次切割的區(qū)域場R1。用R1取代原來的異常場，再用式(1)～(3)計算得到第二次切割的區(qū)域場R2。如此迭代，最后有

實際切割次數(shù)約50～100次。最后得到局部場L，L=g-R。

由于切割算子中的α系數(shù)引入了半二階差分量，它們與負二階水平導數(shù)即曲率成正比，因此插值切割對不同非線性異常部分有不同的作用，從而提高了對不同特征異常場分離的分辨率。

1.2 分離效果對比

為了說明插值切割方法的分離效果，對理論模型數(shù)據(jù)和實際數(shù)據(jù)進行了異常分離，并與上延法和趨勢分析等常用異常分離方法進行了對比。圖2a為理論局部磁異常，圖2b圖為理論區(qū)域磁異常，二者相加得到理論的局部磁異常與區(qū)域磁異常的疊加異常(圖2c)。

a—理論局部磁異常；b—理論區(qū)域磁異常；c—理論疊加磁異常a—theoretical local magnetic anomaly；b—theoretical regional magnetic anomaly；c—theoretical superposition magnetic anomaly

為了從疊加異常中分離出局部異常與區(qū)域異常，利用插值切割方法進行了不同切割半徑下的分離實驗。當切割半徑為150 m、切割次數(shù)100次時，分離出來的局部磁異常與理論局部磁異常無論在異常幅值以及形態(tài)上都具有極大的一致性，分離效果極佳。區(qū)域異常的分離效果也很好，形態(tài)和幅值上均相當，只是分離處理的區(qū)域磁異常的梯度較理論模型小一些。圖3展示了利用插值切割(圖3b)、向上延拓(圖3c)、趨勢分析(圖3d)等方法分離的局部和區(qū)域磁異常與理論異常(圖3a)的對比情況，其中向上延拓分離時的區(qū)域場為上延250 m確定，趨勢分析分離的區(qū)域場為1階多項式確定。實驗表明，插值切割方法有對所設(shè)計模型的局部和區(qū)域場進行分離的能力，較一般方法有優(yōu)越性。

圖3展示了插值切割方法對理論模型異常的分離能力，下面再簡單對比一下實際資料的處理效果。東北某區(qū)磁異常比較發(fā)育，火成巖的磁異常疊加于區(qū)域場上，分別利用插值切割、向上延拓和趨勢分析3種方法進行火成巖引起的局部磁異常的分離。圖4展示了該區(qū)化極磁異常(圖4a)以及插值切割(圖4b)、向上延拓(圖4c)、趨勢分析(圖4d)方法分離的區(qū)域場和局部場。其中插值切割選取的切割半徑為1 km、切割次數(shù)為100次，向上延拓分離的區(qū)域場為上延1 km，趨勢分析分離的區(qū)域場的多項式階次為2階。對比局部場的分離效果，插值切割分離的局部磁異常在異常細節(jié)表現(xiàn)力和分辨率等方面都優(yōu)于其他兩種常用方法。

a—某區(qū)化極磁異常；b—差值切割分離的區(qū)域場(上)及局部場(下)；c—向上延拓分離的區(qū)域場(上)及局部場(下)；d—趨勢分析分離的區(qū)域場(上)及局部場(下)a—abnormal magnetic polarization of a certain region；b—regional field (upper) and local field (lower) separated by interpolation cutting；c—regional field (upper) and local field (lower) separated by upward extension；d—trend analysis separate regional field (upper) and local field (lower)

a—理論區(qū)域場(上)及局部區(qū)域場(下)；b—差值切割分離的區(qū)域場(上)及局部場(下)；c—向上延拓分離的區(qū)域場(上)及局部場(下)；d—趨勢分析分離的區(qū)域場(上)及局部場(下)a—theoretical regional field (upper) and local field (lower);b—the regional field (upper) and local field (lower) separated by interpotation cutting；c—continue upward to separate the regional field (upper) and local (lower)；d—regional field (upper) and partial (lower) separated by trend analysis

2 改進方案及效果

插值切割方法雖然在異常分離上較一般方法具有一定的優(yōu)勢，但仍存在一些問題。筆者采用模型的計算來說明并給出改進前后的對比。圖5為深淺不同磁性磁源的疊合模型，3個綠色方塊(磁性0.05SI)為淺部磁源，產(chǎn)生局部磁異常，一個紫色橢圓(磁性0.1 SI)為深部磁源，產(chǎn)生區(qū)域磁異常。疊加異常如圖6所示。

圖5 深淺源不同磁性疊合模型Fig.5 Different magnetic superposition models of deep and shallow sources

圖6 疊合模型產(chǎn)生的磁異常Fig.6 Magnetic anomaly generated by superposition model

2.1 區(qū)域場圓滑

模型和實際數(shù)據(jù)應(yīng)用中發(fā)現(xiàn)，切割算子獲得的區(qū)域場往往會帶有高頻的抖動，這與區(qū)域場的光滑性相矛盾。針對這一現(xiàn)象，在獲得區(qū)域場后，采用切割窗口大小的滑動平均法進行光滑獲得區(qū)域場，進而獲得局部場。由圖7可看出圓滑前后對比情況，區(qū)域場的光滑性得到改善，進而分離的局部場也不會引入額外的高頻噪聲，區(qū)域場圓滑后獲得的局部場明顯得到改善。

a—圓滑前區(qū)域場(上)及局部場(下)；b—圓滑后區(qū)域場(上)及局部場(下)a—smooth front regional field (upper) and local field (lower)；b—smooth posterior regional field (upper) and local field (lower)

2.2 八方位切割

插值切割的算法是四方位(EW—SN)切割，實際應(yīng)用中發(fā)現(xiàn)四方位切割會引起區(qū)域場的整體形態(tài)畸變，針對這個問題，采用八方位切割來解決(圖8)，即在原四方位(EW-SN)切割的基礎(chǔ)上，增加NE—NW正交方向的切割運算，半徑上的異常值(g1′、g2′、g3′、g4′)采用雙線性插值方法確定，區(qū)域場值g0′采用與四方位切割相同的公式計算，這里不進行贅述。最終的區(qū)域場值采用兩組方向(EW—SN和NE—NW)切割算子的平均值(g0′/2+g0/2)代替。

圖8 八方位切割原理示意Fig.8 Eight-position interpolation cutting method

圖9展示了八方位切割前后對比情況。a上圖為四方位切割的區(qū)域場，區(qū)域場形態(tài)趨近方形，b上圖為八方位切割的區(qū)域場，區(qū)域場的形態(tài)更接近于橢圓形，與源的形狀相同。a下圖和b下圖分別為四方位、八方位切割的局部場。

a—四方位切割的區(qū)域場(上)及局部場(下)；b—八方位切割的區(qū)域場(上)及局部場(下)a—regional field (upper) and local field (lower) cut by the square；b—regional field (upper) and local field (lower) cut by octagonal position

2.3 任意實數(shù)半徑

原插值切割算法的半徑一般按點距的整數(shù)倍(如1dx、2dx、…)進行選擇，這種設(shè)置方法在去噪處理中不適用，過大的切割半徑不僅會消除噪聲，也會對有效信號進行過度濾除。針對這個問題，采用實數(shù)倍點距(如0.1dx、1.3dx、…)作為半徑，相應(yīng)半徑上的異常值采用雙線性插值方法確定。圖10a為加入噪聲(均值0，標準差1 nT)的疊加磁異常圖。從圖10b、c中局部異常形態(tài)和幅值的對比可以明顯看出，利用0.5倍點距半徑濾波后的磁異常不僅去除了噪聲，而且對有效信號的幅值最小程度地削弱。

a—加入噪聲的疊加磁異常；b—1倍點距半徑濾波后的磁異常；c—0.5倍點距半徑濾波后的磁異常a—superposition magnetic anomaly with added noise；b—magnetic anomaly after 1 radius of point distance filtering；c—magnetic anomaly after 0.5 times of point distance radius filtering

3 實際應(yīng)用效果

3.1 重力異常半定量分離

中國西部某區(qū)域地表復雜，地震地質(zhì)條件極差，地震資料的品質(zhì)較差，難以解決問題。為了加快該區(qū)油氣勘探的步伐，于2007年部署高精度重磁勘探，測網(wǎng)500 m×200 m，重力異常精度0.027 mGal，磁異常精度1.26 nT。2011年經(jīng)過鉆井和地震聯(lián)合攻關(guān)，確定了1 500 m左右深度的淺部構(gòu)造，構(gòu)造形態(tài)見圖11中的黑色虛線。

從圖11a可見，地震構(gòu)造沒有與布格重力異常的任何高異常點吻合，并且偏差很大。為了研究插值切割方法對不同深度異常的提取能力，利用不同半徑提取的局部重力異常獲得了代表不同深度段的重力異常圖，并通過向下延拓方法將不同深度異常延拓至對應(yīng)深度頂面。圖11中b、c、d、e、f、g分別代表200～600 m、600～1 000 m、1 000～1 400 m、1 400～1 800 m、1 800～2 200 m、2 200～2 600 m深度段的重力異常。從這些異常與地震構(gòu)造的對比來看，圖11e是1 400～1 800 m深度的重力異常與地震構(gòu)造吻合最好，而且深度也基本吻合。這反映了插值切割算法的異常半定量分離能力。

a—布格重力異常；b—200～600 m深度重力異常；c—600～1 000 m深度重力異常；d—1 000～1 400 m深度重力異常；e—1 400～1 800 m深度重力異常；f—1 800～2 200 m深度重力異常；g—2 200～2 600 m深度重力異常a—bouguer gravity anomaly; b—gravity anomaly at depths of 200～600 m; c—gravity anomaly at depths of 600～1 000 m; d—gravity anomaly at depths of 1 000～1 400 m; e—gravity anomaly at depths of 1 400～1 800 m; f—gravity anomaly at depths of 1 800～2 200 m; g—gravity anomaly at depths of 2 200～2 600 m

3.2 磁異常三維成像

南祁連盆地位于祁連山南部，南部隔宗務(wù)隆山與柴達木盆地相鄰，西北部緊鄰阿爾金山，北東通過中祁連地塊與河西走廊相望(圖12)。該盆地是在加里東褶皺基底上發(fā)育起來的山間盆地，石炭系—三疊系為地臺期一套相對穩(wěn)定的淺海相和海陸交互相沉積，侏羅系—白堊系缺失，上新近系—第四系為陸相沉積。研究區(qū)的磁性體主要發(fā)育在志留系，巖性為千枚狀凝灰?guī)r、千枚巖和安山巖等，磁化率達到 (1 000～1 600)×10-5SI。磁異常主要是志留系基底內(nèi)部磁性體的反映。ΔT磁異常經(jīng)過化極處理后得到了化極磁異常圖，見圖13，它主要反映了志留系磁性體的分布。

圖13 化極磁異常Fig.13 RTP magnetic anomaly

圖12 研究區(qū)構(gòu)造位置Fig.12 Schematic location of study area

常規(guī)的重磁資料處理，可以獲得較為可靠的構(gòu)造平面展布信息，深度計算一般采用剖面正反演或基于規(guī)則模型的特征點法，它們的精度或依賴于已知控制點成果，或依賴于異常與典型模型的符合程度，而且由于計算和反演周期較長，往往僅作數(shù)條剖面。為了解決上述問題，我們引入了插值切割法，通過插值切割三維成像研究區(qū)內(nèi)火成巖的三維空間分布。

根據(jù)前人關(guān)于插值切割數(shù)據(jù)與地質(zhì)體深度密切相關(guān)的研究成果，我們把不同切割半徑的插值切割數(shù)據(jù)賦予相應(yīng)的深度概念，構(gòu)建三維數(shù)據(jù)體，在一定程度上反映了地質(zhì)體的三維空間特征，即地質(zhì)體三維磁成像。三維成像數(shù)據(jù)體構(gòu)建采用以下步驟：

1) 以異常數(shù)據(jù)網(wǎng)格距及其整倍數(shù)作為切割半徑,采用插值切割法分別對異常場求取局部場和區(qū)域場；

2) 以對應(yīng)各切割半徑的局部場(也可用減去上一切割半徑的局部場的插值)作為異常，以切割半徑作為深度，建立各深度的異常數(shù)據(jù)體；

3) 將各深度的異常數(shù)據(jù)體組成三維數(shù)據(jù)體；

4) 采用三維顯示軟件對異常進行顯示，就獲得異常三維成像。

圖14是插值切割三維磁成像，圖中紅色表示磁異常強度大，一定程度上代表了磁化率高，異?？臻g展布表示了磁性體的三維空間分布?？梢钥闯?，三維成像很好地反映了該區(qū)火成巖的空間分布特征。

圖14 插值切割三維磁成像效果Fig.14 Three-dimensional magnetic imaging effect of interpolation cutting

4 結(jié)論

從區(qū)域場圓滑、八方位切割和任意實數(shù)半徑等幾個方面對插值切割算法進行了完善，較有效地解決了應(yīng)用中的部分問題，模型數(shù)據(jù)計算的對比分析表明了改進的有效性。

從重力異常半定量分離和磁異常三維成像兩個應(yīng)用實例來看，插值切割在重磁異常精細分離和三維快速成像方面應(yīng)用效果明顯，在缺少其他定量分離和反演手段時，不失為一種替代方法。

需要指出的是，該方法依然是一種半定量方法，主要用于研究地質(zhì)體空間展布特征，具體的深度會存在一定的誤差，這在應(yīng)用中需要注意。