位場(chǎng)數(shù)據(jù)處理算法的數(shù)據(jù)挖掘試驗(yàn)與應(yīng)用

魏華敬， 尹宏偉， 姜素華， 汪 剛， 趙斐宇

(1. 南京大學(xué) 地球科學(xué)與工程學(xué)院；能源科學(xué)研究院，南京 210023；2. 中國(guó)海洋大學(xué) 海洋地球科學(xué)學(xué)院； 海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東 青島 266100)

0 引 言

數(shù)據(jù)挖掘是從海量數(shù)據(jù)中挖掘知識(shí)的技術(shù)，隨著大數(shù)據(jù)時(shí)代的到來(lái)，數(shù)據(jù)挖掘發(fā)展迅猛。作為一種通用的技術(shù)，數(shù)據(jù)挖掘在各行業(yè)都扮演著舉足輕重的作用[1]。地球物理學(xué)是通過(guò)定量的物理方法研究地球以及尋找地球內(nèi)部礦產(chǎn)資源的學(xué)科，重力勘探是地球物理的方法之一，通過(guò)地表測(cè)得的重力異常確定地表以下的構(gòu)造或地質(zhì)體位置的方法[2]。地球物理數(shù)據(jù)作為一種可挖掘的數(shù)據(jù)，讓地質(zhì)與地球物理學(xué)界的專家學(xué)者獲益匪淺。

位場(chǎng)是一種場(chǎng)的大小與位置有關(guān)的場(chǎng)，重力場(chǎng)是位場(chǎng)的一種，通過(guò)對(duì)位場(chǎng)的變換和處理可以挖掘出場(chǎng)源相關(guān)的信息。然而位場(chǎng)數(shù)據(jù)直接反映場(chǎng)源的能力較差，需要對(duì)位場(chǎng)數(shù)據(jù)進(jìn)行適當(dāng)?shù)淖儞Q。國(guó)內(nèi)外學(xué)者對(duì)位場(chǎng)數(shù)據(jù)的處理方法有一定的研究，Cooper等[3]提出的THDR、Theta Map、HTA等算法提高了邊界增強(qiáng)后圖像識(shí)別的效果；張超等[4]提出的Sigmoid算法實(shí)現(xiàn)異常值網(wǎng)格數(shù)據(jù)的拉升和灰度級(jí)像素的壓低，凸顯了地質(zhì)體的邊界；張沖等[5]提出向下延拓3階Adams-Bashforth公式法，相比起傳統(tǒng)的延拓方法更穩(wěn)定，不容易產(chǎn)生邊界效應(yīng)，使延拓結(jié)果更準(zhǔn)確。

本文設(shè)計(jì)了位場(chǎng)數(shù)據(jù)處理的Matlab算法[6]：梯度算法，解析延拓算法，邊界提取算法，并對(duì)規(guī)則形體的重力異常使用上述3種方法進(jìn)行處理，討論了它們各自的特點(diǎn)。

1 位場(chǎng)數(shù)據(jù)處理算法

在重力位場(chǎng)中，引力位V的1階垂向?qū)?shù)為重力異常，則不同的位場(chǎng)數(shù)據(jù)處理算法。

梯度算法

VZ垂向?qū)?shù)：

(1)

VZ水平0°導(dǎo)數(shù)：

(2)

VZ水平45°導(dǎo)數(shù)：

Vz45(x,y)≈

(3)

VZ水平90°導(dǎo)數(shù)：

(4)

VZ水平135°導(dǎo)數(shù)：

Vz135(x,y)≈

(5)

解析延拓算法

向上延拓：

(6)

向下延拓：

(7)

邊界提取算法

ReLU：

max(Vzz(x,y),0)

(8)

Leaky ReLU：

max(Vzz(x,y),0.01Vzz(x,y))

(9)

tanh：

tanh(Vzz(x,y)/n)

(10)

sign：

sign(Vzz(x,y))

(11)

2 模型試算

為了驗(yàn)證算法較位場(chǎng)原始數(shù)據(jù)的優(yōu)越性以及對(duì)場(chǎng)源體的識(shí)別效果，建立了若干模型進(jìn)行計(jì)算分析。本文選定地面觀測(cè)數(shù)據(jù)區(qū)域?yàn)?1×51點(diǎn)距的平面網(wǎng)格面積，并通過(guò)規(guī)則形體的重力異常正演公式計(jì)算得到重力異常數(shù)據(jù)[7]。對(duì)模型進(jìn)行以下幾種算法的處理，探討算法各自的特點(diǎn)(見(jiàn)圖1)。圖中：Vz表示重力異常；Vzz表示重力異常的1階垂向?qū)?shù)；Vzzz表示重力異常的2階垂向?qū)?shù)；Vz0,Vz45,Vz90,Vz135分別表示重力異常方位角為0°，45°，90°，135°的水平導(dǎo)數(shù)；VzUP和VzDOWN分別表示重力異常向上延拓和向下延拓的計(jì)算結(jié)果。圖中比例尺中的單位:Vz,VzUP,VzDOWN的單位為mGal (1 mGal=10-5m/s2)；Vz0,Vz45,Vz90,Vz135和Vzz的單位為E(1E=10-9/s2)；Vzzz的單位為nMKS (1nMKS=10-12/(m·s2))

2.1 梯度算法

模型1和2形狀如圖1(a)、(b)所示，幾何形態(tài)截然不同的兩個(gè)模型具有相似的重力異常特征(圖1(c)、(d))，模型的垂向1階和2階梯度 (圖1(e)～(h))表現(xiàn)出了明顯的差異，梯度階次的增加使異常平面圖表現(xiàn)出更明顯的低頻特性，異常的平面特征與場(chǎng)源體的幾何形態(tài)趨于吻合。模型1和2在0°、45°、90°、135°4個(gè)方向的水平1階梯度如圖2所示。場(chǎng)源體處水平梯度產(chǎn)生明顯異常條帶(圖2(a)、(c)、(e)、(g))。

模型3(圖3(a))是兩個(gè)水平位置靠近的條帶狀場(chǎng)源體，其重力異常(圖3(c))和單個(gè)條帶狀場(chǎng)源體模型4(圖3(b))產(chǎn)生的異常(圖3(d))的平面特征相近，基本不可直接區(qū)分。垂向1階梯度和2階梯度(圖3(e)、(g))能將模型3的相互靠近的兩個(gè)異常場(chǎng)分離，且2階分離效果好于1階。

圖1 模型1

圖2 模型2

圖3 模型3

(a)及其Vz(c),Vzz(e),Vzzz(g);(b)及其Vz(d),Vzz(f),Vzzz(h)

模型5(圖4)由3個(gè)埋深互不相同的條帶狀場(chǎng)源體組成，模型的垂向1階梯度和2階梯度 (圖6(b),(c))隨著階次的增加，異常的低頻成分愈發(fā)明顯，這一規(guī)律有利于分離不同埋深的異常體。

圖4 模型5平面圖和立體圖

從梯度算法模型試算的結(jié)果可以看出，不同形態(tài)的場(chǎng)源體的梯度具有不同的平面特征。梯度算法不僅能將相互靠近的場(chǎng)源產(chǎn)生的疊加異常分離，還能突出淺部場(chǎng)源壓制深部，且梯度階次越高效果越明顯。實(shí)測(cè)重力異常往往是多場(chǎng)源疊加的結(jié)果，梯度算法能夠分離不同性質(zhì)的場(chǎng)源產(chǎn)生的異常，有利于分類和解釋。

(a) Vz(b) Vzz

(c) Vzzz

圖6 模型6平面圖和立體圖

2.2 解析延拓算法

模型4由4個(gè)不同埋深的條帶狀場(chǎng)源組成(見(jiàn)圖6)，模型的重力異常(圖7(b))表現(xiàn)出與埋深呈現(xiàn)負(fù)相關(guān)的態(tài)勢(shì)：由淺到深的場(chǎng)源產(chǎn)生由大到小的異常，對(duì)應(yīng)云圖上的顏色為由深到淺。云圖7(a)由異常值引起的顏色差異較云圖7(b)不明顯，模型重力異常的向上延拓結(jié)果(圖7(a))表現(xiàn)出異常值隨深度的差異較延拓前(圖7(b))不明顯，相對(duì)地突出了深部場(chǎng)源的產(chǎn)生的異常，淺部異常受壓制。反之，圖7(c)中深淺場(chǎng)源引起的顏色差異顯著，表示淺部場(chǎng)被突出，深部場(chǎng)被壓制。

(a) VzUP(b) Vz

(c) VzDOWN

從解析延拓算法模型試算結(jié)果可以得到結(jié)論：向下延拓突出淺部異常體，向上延拓突出深部異常體。

2.3 邊界提取算法

邊界提取算法是增強(qiáng)位場(chǎng)邊界識(shí)別和檢測(cè)效果的一種方法[8]。位場(chǎng)垂向梯度在場(chǎng)源邊界處發(fā)生急劇變化，通過(guò)tanh、sign、ReLU、Leaky ReLU等激活函數(shù)(見(jiàn)圖8)的作用能夠一定程度反映場(chǎng)源體的邊界。激活函數(shù)在深度學(xué)習(xí)與神經(jīng)網(wǎng)絡(luò)當(dāng)中十分常見(jiàn)，能夠一定程度地將線性模型轉(zhuǎn)化為非線性，豐富模型內(nèi)容[9]。

對(duì)于tanh函數(shù)，若自變量值絕對(duì)值過(guò)大，則函數(shù)變換后的值通常為1或-1，tanh函數(shù)邊界提取的效果與sign函數(shù)近似，需要通過(guò)引入合適的參數(shù)n使自變量絕對(duì)值縮小。自變量落在0附近函數(shù)變化率大的范圍內(nèi)能夠更好地由線性信息得到非線性信息，從而使算法更加實(shí)用。引入?yún)?shù)后的改良tanh函數(shù)為：

tanh(Vzz(x,y)/n)

對(duì)模型5的重力異常作邊界提取處理，從圖9對(duì)比可以看出，激活函數(shù)對(duì)邊界的提取效果各有千秋。改良tanh函數(shù)用梯度帶來(lái)反映邊界，sign函數(shù)在邊界附近表現(xiàn)出符號(hào)階躍，ReLU和Laeky ReLU用接近0的極值條帶反映了場(chǎng)源覆蓋的范圍，0值與正值的交界能一定程度反映邊界。4種激活函數(shù)提取出的模型邊界略大于實(shí)際邊界，邊界的寬度能反映場(chǎng)源的深淺。

圖8 激活函數(shù)

圖9 模型5邊界提取

3 位場(chǎng)數(shù)據(jù)挖掘?qū)嵗?/h2>

位于我國(guó)西北部新疆的塔里木盆地蘊(yùn)藏著豐厚的石油及天然氣資源[10]。塔里木盆地經(jīng)緯度范圍為75°E～9E°，36°E～42°，盆地?cái)嗔严到y(tǒng)發(fā)育十分良好(圖10)，且大多數(shù)油氣藏均與斷裂系統(tǒng)有著密不可分的關(guān)聯(lián)。本文結(jié)合數(shù)據(jù)挖掘的基本思路將位場(chǎng)處理算法應(yīng)用于塔里木盆地的布格重力異常數(shù)據(jù)，得到塔里木盆地?cái)嗔严到y(tǒng)的相關(guān)信息。實(shí)驗(yàn)數(shù)據(jù)來(lái)源于BGI(Bureau Gravimetric International)[11]，數(shù)據(jù)精度為2′，數(shù)據(jù)經(jīng)緯度為72°E～92°E，34°E～44°。

圖10 塔里木盆地主要斷裂示意圖(據(jù)文獻(xiàn)[13]修改)

3.1 數(shù)據(jù)清理、集成與選擇

BGI將實(shí)測(cè)重力數(shù)據(jù)通過(guò)基點(diǎn)聯(lián)測(cè)或平滑處理消除地表附近密度不均因素引起的異常。BGI將不同來(lái)源的重力數(shù)據(jù)(如地面重力測(cè)量，航空重力測(cè)量等)組合形成數(shù)據(jù)庫(kù)。實(shí)測(cè)重力數(shù)據(jù)經(jīng)過(guò)地形校正、中間層校正、自由空氣校正、均衡校正等得到布格異常、自由空氣異常、均衡異常，均存放在數(shù)據(jù)庫(kù)中，學(xué)者根據(jù)研究目的選出相應(yīng)的數(shù)據(jù)。本文選擇能夠反映地球內(nèi)部剩余密度分布的布格重力異常?？紤]到對(duì)位場(chǎng)數(shù)據(jù)做處理時(shí)的邊界效應(yīng)，通常使用數(shù)據(jù)的經(jīng)緯度范圍比研究區(qū)域大2°。

3.2 數(shù)據(jù)變換

數(shù)據(jù)變換是為了讓數(shù)據(jù)變換成適合挖掘的形式，本例中使用梯度算法、解析延拓算法、邊界提取算法將布格重力異常變換成適合解釋的形式。

通常，在重力異常和垂向梯度平面圖中梯級(jí)帶、不同特征異常區(qū)的分界線、線性分布的高低異常過(guò)渡帶能反映深大斷裂，有時(shí)也可以反映大范圍不同巖性的接觸帶[12]。在水平梯度平面圖中，異常條帶能夠反映深大斷裂的位置。向上和向下解析延拓的方法可以分別突出淺部或深部的異常特征。進(jìn)行延拓處理后再根據(jù)斷裂構(gòu)造在重力異常平面等值線圖上的特征對(duì)斷裂構(gòu)造加以識(shí)別。

3.3 數(shù)據(jù)挖掘

本例中數(shù)據(jù)挖掘步驟是對(duì)算法處理過(guò)的數(shù)據(jù)進(jìn)行相應(yīng)的解釋。刪除盆地邊界以外的數(shù)據(jù)點(diǎn)，僅留下與盆地形狀大小相近的區(qū)域作為處理結(jié)果并可視化。根據(jù)3.2中提到判別標(biāo)準(zhǔn)可以確定出斷層的平面展布的情況。

垂向1階(圖11(a))和2階梯度(圖11(b))與研究區(qū)域的布格重力異常場(chǎng)(圖12)相比，高頻成分豐富，表現(xiàn)出的斷裂的數(shù)量多且細(xì)致，突出了埋深較淺的斷裂。水平梯度 (圖13)以極值條帶的方式體現(xiàn)了斷裂的展布，圖中只確定了條帶較寬在圖中較為明顯的異常條帶，它們和深且大的斷裂相對(duì)應(yīng)。研究區(qū)域中不同斷裂的走向和規(guī)模不同，不同方向的水平梯度突出的深大斷裂也不相同，因此重力異常數(shù)據(jù)處理當(dāng)中常對(duì)多個(gè)方向進(jìn)行水平梯度計(jì)算，綜合考慮各方向處理結(jié)果才能得到合適的解釋。

(a) 垂向1階梯度Vzz

(c) 向上延拓14.8 kmVzUP

(d) 向下延拓14.8 kmVzZOWN

圖12 塔里木盆地布格重力異常

(a) 水平梯度Vz0

(b) Vz45

(c) Vz90

(d) Vz135

解析延拓中的單位長(zhǎng)度就是原始數(shù)據(jù)的精度，本文使用的原始數(shù)據(jù)的精度是2′，對(duì)應(yīng)的長(zhǎng)度為3.7 km。分別對(duì)研究區(qū)域布格重力異常場(chǎng)(圖12)向上延拓4個(gè)單位長(zhǎng)度(14.8 km)(圖11(c))和向下延拓4個(gè)單位長(zhǎng)度(14.8 km)(圖11(d))。向上延拓所得結(jié)果中長(zhǎng)波段信號(hào)豐富，反映大型斷裂的展布。向下延拓極大程度地突出了高頻成分，它們是由淺部小規(guī)模的異常源引起，淺部斷裂數(shù)量多展布方式復(fù)雜，圖中沒(méi)有一一描繪出斷裂形態(tài)。

用tanh, sign, ReLU, Leaky ReLU 4種激活函數(shù)處理研究區(qū)域的布格重力異常數(shù)據(jù)。由邊界提取算法模型試算的結(jié)果可知，tanh函數(shù)(圖14(a))梯度帶能表征斷裂邊界，邊界在sign函數(shù)(圖14(b))表現(xiàn)為階躍，ReLU(圖14(c))和Leaky ReLU(圖14(d))0值附近的極小值條帶即圖中黃色的條帶是深大斷裂存在的區(qū)域。

(a) tanh

(b) sign

(c) ReLU

(d) Leaky ReLU

3.4 模式評(píng)估

區(qū)域主要斷裂示意圖(圖10)可識(shí)別挖掘所得結(jié)論的可靠性：盆地主要由五大斷裂系統(tǒng)組成，它們分別是天山山前斷裂系統(tǒng)、庫(kù)魯克塔格斷裂系統(tǒng)、阿爾金山前斷裂系統(tǒng)、昆山山前斷裂系統(tǒng)和巴楚凸起斷裂系統(tǒng)。其中巴楚凸起斷裂系統(tǒng)位于盆地內(nèi)部，另外4個(gè)斷裂系統(tǒng)位于盆地的邊緣，走向?yàn)镹E、NEE和NW，且基本與盆地邊緣平行。本例挖掘出的斷裂系統(tǒng)的分布模式與地質(zhì)事實(shí)吻合度高。

運(yùn)用不同的挖掘方式挖掘同樣的數(shù)據(jù)能挖掘出不同的知識(shí)[14]，使用不同的算法處理區(qū)域異常數(shù)據(jù)既能得到構(gòu)造水平方向的位置信息也能得到垂向的深度信息。對(duì)不同的數(shù)據(jù)用同樣的挖掘方式也能得到不同結(jié)論。除本文介紹的應(yīng)用之外，位場(chǎng)數(shù)據(jù)處理技術(shù)還廣泛應(yīng)用于磁力探測(cè)、遙感、計(jì)算機(jī)視覺(jué)、特征學(xué)習(xí)等領(lǐng)域當(dāng)中[15]。

4 結(jié) 語(yǔ)

本文通過(guò)模型試算比較了幾種位場(chǎng)數(shù)據(jù)處理算法的特點(diǎn)：梯度算法能夠區(qū)別不同形狀的場(chǎng)源，突出淺部短波段信號(hào)，并分離靠近的場(chǎng)源異常，且梯度階數(shù)越高效果越明顯；解析延拓算法能夠突出不同深度的場(chǎng)源異常，向上延拓能突出深部場(chǎng)源，向下則突出淺部；邊界提取算法以梯度帶，階躍等形式刻畫(huà)出場(chǎng)源邊界。通過(guò)挖掘地球物理數(shù)據(jù)可得到地表以下的情況，塔里木盆地?cái)嗔严到y(tǒng)尤其是盆地邊緣的深大斷裂在梯度算法的作用下主要表現(xiàn)為梯度帶和條帶狀異常，斷裂的深淺不一在解析延拓算法處理結(jié)果中表現(xiàn)為長(zhǎng)短波段，深大斷裂在向上延拓圖中表現(xiàn)為低頻信號(hào)，向下延拓能夠識(shí)別和圈定一些細(xì)小的異常特征。塔里木盆地的主要斷裂系統(tǒng)走向?yàn)镹E, NEE, NW且基本與盆地邊緣平行分布，與研究區(qū)域?qū)嶋H情況相符。