位場數據處理算法的數據挖掘試驗與應用

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

(1. 南京大學 地球科學與工程學院；能源科學研究院，南京 210023；2. 中國海洋大學 海洋地球科學學院； 海底科學與探測技術教育部重點實驗室，山東 青島 266100)

0 引 言

數據挖掘是從海量數據中挖掘知識的技術，隨著大數據時代的到來，數據挖掘發展迅猛。作為一種通用的技術，數據挖掘在各行業都扮演著舉足輕重的作用[1]。地球物理學是通過定量的物理方法研究地球以及尋找地球內部礦產資源的學科，重力勘探是地球物理的方法之一，通過地表測得的重力異常確定地表以下的構造或地質體位置的方法[2]。地球物理數據作為一種可挖掘的數據，讓地質與地球物理學界的專家學者獲益匪淺。

位場是一種場的大小與位置有關的場，重力場是位場的一種，通過對位場的變換和處理可以挖掘出場源相關的信息。然而位場數據直接反映場源的能力較差，需要對位場數據進行適當的變換。國內外學者對位場數據的處理方法有一定的研究，Cooper等[3]提出的THDR、Theta Map、HTA等算法提高了邊界增強后圖像識別的效果；張超等[4]提出的Sigmoid算法實現異常值網格數據的拉升和灰度級像素的壓低，凸顯了地質體的邊界；張沖等[5]提出向下延拓3階Adams-Bashforth公式法，相比起傳統的延拓方法更穩定，不容易產生邊界效應，使延拓結果更準確。

本文設計了位場數據處理的Matlab算法[6]：梯度算法，解析延拓算法，邊界提取算法，并對規則形體的重力異常使用上述3種方法進行處理，討論了它們各自的特點。

1 位場數據處理算法

在重力位場中，引力位V的1階垂向導數為重力異常，則不同的位場數據處理算法。

梯度算法

VZ垂向導數：

(1)

VZ水平0°導數：

(2)

VZ水平45°導數：

Vz45(x,y)≈

(3)

VZ水平90°導數：

(4)

VZ水平135°導數：

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 模型試算

為了驗證算法較位場原始數據的優越性以及對場源體的識別效果，建立了若干模型進行計算分析。本文選定地面觀測數據區域為51×51點距的平面網格面積，并通過規則形體的重力異常正演公式計算得到重力異常數據[7]。對模型進行以下幾種算法的處理，探討算法各自的特點(見圖1)。圖中：Vz表示重力異常；Vzz表示重力異常的1階垂向導數；Vzzz表示重力異常的2階垂向導數；Vz0,Vz45,Vz90,Vz135分別表示重力異常方位角為0°，45°，90°，135°的水平導數；VzUP和VzDOWN分別表示重力異常向上延拓和向下延拓的計算結果。圖中比例尺中的單位: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)所示，幾何形態截然不同的兩個模型具有相似的重力異常特征(圖1(c)、(d))，模型的垂向1階和2階梯度 (圖1(e)～(h))表現出了明顯的差異，梯度階次的增加使異常平面圖表現出更明顯的低頻特性，異常的平面特征與場源體的幾何形態趨于吻合。模型1和2在0°、45°、90°、135°4個方向的水平1階梯度如圖2所示。場源體處水平梯度產生明顯異常條帶(圖2(a)、(c)、(e)、(g))。

模型3(圖3(a))是兩個水平位置靠近的條帶狀場源體，其重力異常(圖3(c))和單個條帶狀場源體模型4(圖3(b))產生的異常(圖3(d))的平面特征相近，基本不可直接區分。垂向1階梯度和2階梯度(圖3(e)、(g))能將模型3的相互靠近的兩個異常場分離，且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個埋深互不相同的條帶狀場源體組成，模型的垂向1階梯度和2階梯度 (圖6(b),(c))隨著階次的增加，異常的低頻成分愈發明顯，這一規律有利于分離不同埋深的異常體。

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

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

(a) Vz(b) Vzz

(c) Vzzz

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

2.2 解析延拓算法

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

(a) VzUP(b) Vz

(c) VzDOWN

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

2.3 邊界提取算法

邊界提取算法是增強位場邊界識別和檢測效果的一種方法[8]。位場垂向梯度在場源邊界處發生急劇變化，通過tanh、sign、ReLU、Leaky ReLU等激活函數(見圖8)的作用能夠一定程度反映場源體的邊界。激活函數在深度學習與神經網絡當中十分常見，能夠一定程度地將線性模型轉化為非線性，豐富模型內容[9]。

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

tanh(Vzz(x,y)/n)

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

圖8 激活函數

圖9 模型5邊界提取

3 位場數據挖掘實例

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

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

3.1 數據清理、集成與選擇

BGI將實測重力數據通過基點聯測或平滑處理消除地表附近密度不均因素引起的異常。BGI將不同來源的重力數據(如地面重力測量，航空重力測量等)組合形成數據庫。實測重力數據經過地形校正、中間層校正、自由空氣校正、均衡校正等得到布格異常、自由空氣異常、均衡異常，均存放在數據庫中，學者根據研究目的選出相應的數據。本文選擇能夠反映地球內部剩余密度分布的布格重力異常。考慮到對位場數據做處理時的邊界效應，通常使用數據的經緯度范圍比研究區域大2°。

3.2 數據變換

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

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

3.3 數據挖掘

本例中數據挖掘步驟是對算法處理過的數據進行相應的解釋。刪除盆地邊界以外的數據點，僅留下與盆地形狀大小相近的區域作為處理結果并可視化。根據3.2中提到判別標準可以確定出斷層的平面展布的情況。

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

(a) 垂向1階梯度Vzz

(c) 向上延拓14.8 kmVzUP

(d) 向下延拓14.8 kmVzZOWN

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

(a) 水平梯度Vz0

(b) Vz45

(c) Vz90

(d) Vz135

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

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

(a) tanh

(b) sign

(c) ReLU

(d) Leaky ReLU

3.4 模式評估

區域主要斷裂示意圖(圖10)可識別挖掘所得結論的可靠性：盆地主要由五大斷裂系統組成，它們分別是天山山前斷裂系統、庫魯克塔格斷裂系統、阿爾金山前斷裂系統、昆山山前斷裂系統和巴楚凸起斷裂系統。其中巴楚凸起斷裂系統位于盆地內部，另外4個斷裂系統位于盆地的邊緣，走向為NE、NEE和NW，且基本與盆地邊緣平行。本例挖掘出的斷裂系統的分布模式與地質事實吻合度高。

運用不同的挖掘方式挖掘同樣的數據能挖掘出不同的知識[14]，使用不同的算法處理區域異常數據既能得到構造水平方向的位置信息也能得到垂向的深度信息。對不同的數據用同樣的挖掘方式也能得到不同結論。除本文介紹的應用之外，位場數據處理技術還廣泛應用于磁力探測、遙感、計算機視覺、特征學習等領域當中[15]。

4 結 語

本文通過模型試算比較了幾種位場數據處理算法的特點：梯度算法能夠區別不同形狀的場源，突出淺部短波段信號，并分離靠近的場源異常，且梯度階數越高效果越明顯；解析延拓算法能夠突出不同深度的場源異常，向上延拓能突出深部場源，向下則突出淺部；邊界提取算法以梯度帶，階躍等形式刻畫出場源邊界。通過挖掘地球物理數據可得到地表以下的情況，塔里木盆地斷裂系統尤其是盆地邊緣的深大斷裂在梯度算法的作用下主要表現為梯度帶和條帶狀異常，斷裂的深淺不一在解析延拓算法處理結果中表現為長短波段，深大斷裂在向上延拓圖中表現為低頻信號，向下延拓能夠識別和圈定一些細小的異常特征。塔里木盆地的主要斷裂系統走向為NE, NEE, NW且基本與盆地邊緣平行分布，與研究區域實際情況相符。