孫伯軒，侯振隆，周文月，鞏恩普，鄭玉君，程浩

(1.東北大學 資源與土木工程學院，遼寧 沈陽 110819；2.中國自然資源航空物探遙感中心，北京 100083)

0 引言

重力梯度數(shù)據(jù)是引力位數(shù)據(jù)沿坐標軸方向的二階導數(shù)，具有分辨率高、抗干擾能力強的特點。近年來基于重力梯度數(shù)據(jù)的密度反演[1]、邊界識別[2-3]和成像方法[4]等研究發(fā)展迅速。

在上述方法中，歐拉反褶積是一種能自動計算場源位置的方法，Thompson于1982年將歐拉反褶積拓展到重力異常數(shù)據(jù)中[5]，但是在歐拉反褶積的應(yīng)用過程中會產(chǎn)生發(fā)散解，降低計算精度，因此有必要改進方法，提高分辨率。許多國內(nèi)外學者圍繞該方法開展了研究并取得良好的效果：Cooper通過計算單個數(shù)據(jù)點的解和選擇場源位置，提高了結(jié)果的準確性[6]；Petar等拓展了Thompson[5]提出的構(gòu)造指數(shù)的范圍[7]；Wang等[8]在歐拉反褶積計算中引入了奇異值分解技術(shù)，降低了噪聲對計算結(jié)果的影響;Eric等[9]嘗試將重力數(shù)據(jù)歐拉反褶積與震源邊緣檢測等方法結(jié)合，可以更清晰地顯示細節(jié)；Zhang等[10]提出了利用重力梯度數(shù)據(jù)進行歐拉反褶積，提高了場源數(shù)據(jù)的解釋精度；Majid等[11]進一步討論了結(jié)構(gòu)指數(shù)的使用，該參數(shù)能夠減少地質(zhì)體定位時的不確定性；周文月等[12]提出了不同高度觀測數(shù)據(jù)的聯(lián)合歐拉反褶積方法，拓展了歐拉反褶積的應(yīng)用范圍；侯振隆等[13]對重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積又做了進一步改進，避免了構(gòu)造指數(shù)的選擇和分量轉(zhuǎn)換，提高了計算的精度。馬國慶等[14]利用重力梯度張量及其導數(shù)改進了聯(lián)合歐拉反褶積算法矩陣，降低背景異常差異的干擾。相對于重力異常數(shù)據(jù)的歐拉反褶積，使用重力梯度數(shù)據(jù)的方法對目標地質(zhì)體具有更高的識別分辨率。

GeoProbe[15]、Geoquest[16]、GeoEast[17]等專業(yè)軟件近年來發(fā)展迅速。許多國內(nèi)學者在重磁處理解釋軟件架構(gòu)與開發(fā)技術(shù)的自主研發(fā)方面做了很多研究：鄭元滿等[18]研發(fā)了針對金屬礦勘探的重磁綜合軟件，利用協(xié)同交互操作提高了軟件工作效率。陳靖等[19]設(shè)計了一種基于組件加插件的分層體系架構(gòu)的重磁軟件，實現(xiàn)多元數(shù)據(jù)集成分析解釋。王浩然[20]基于重磁數(shù)據(jù)物性反演方法設(shè)計了一種卷積神經(jīng)網(wǎng)絡(luò)，降低了重力密度成像所需的訓練量。張浩平等[21]用Visual C++6.0調(diào)用Matlab Engine，以Matlab腳本為核心開發(fā)了一種磁性數(shù)據(jù)的三維歐拉反褶積計算軟件，為相關(guān)的軟件開發(fā)提供了思路。從重磁數(shù)據(jù)處理解釋軟件的研發(fā)現(xiàn)狀可知，當前已開發(fā)的商業(yè)軟件較多，多數(shù)的功能較為完備和復雜，但不利于數(shù)據(jù)的高效解釋和軟件的二次開發(fā)。同時，在數(shù)據(jù)處理解釋方法的研究中，絕大多數(shù)算法是通過編程語言單獨開發(fā)腳本和函數(shù)實現(xiàn)的，缺少一個集成的軟件平臺。

可見，使用具有豐富函數(shù)庫、可移植性高的語言，開發(fā)易用的可視化軟件是十分必要的。這有利于降低方法二次開發(fā)與軟件移植的難度，使研究人員減少對數(shù)據(jù)管理與可視化的開發(fā)成本，擴展解釋方法的適用范圍。因此，本文基于簡潔實用的原則，利用Python語言及其函數(shù)庫開發(fā)了一個支持數(shù)據(jù)/文件管理、二/三維可視化、邊界識別、重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積等功能的軟件系統(tǒng)，為今后的重磁數(shù)據(jù)解釋研究提供技術(shù)支持。

1 方法原理

1.1 重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積

根據(jù)歐拉齊次方程的定義[13]，歐拉反褶積公式如下：

(1)

式中：(x，y，z)與(x0，y0，z0)分別為笛卡爾坐標系中觀測點和地質(zhì)體中心坐標；f和B分別表示觀測異常場和背景場；N為構(gòu)造指數(shù)。將式(1)轉(zhuǎn)化為向量積形式，即：

(2)

用引力位V在x、y、z方向的一階導數(shù)Vx、Vy、Vz代替式(2)中各方向的f。因此，?f/?x、?f/?y、?f/?z在x方向分別轉(zhuǎn)化為Vxx、Vxy、Vxz形式，y軸與z軸方向分量同理。由于求高階導數(shù)影響計算精度，因此，利用式(3)計算引力位的三階垂向?qū)?shù)：

(3)

最后，整理得出重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積公式：

(4)

重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積公式聯(lián)合多個梯度分量建立方程組，將構(gòu)造指數(shù)N與場源坐標一同視為未知數(shù)，有效地避免了選擇構(gòu)造指數(shù)所產(chǎn)生的系統(tǒng)誤差，可以更準確地描述地質(zhì)體空間位置。

1.2 基于方向總水平導數(shù)相關(guān)系數(shù)的邊界識別方法

本文提出一種利用方向總水平導數(shù)生成的相關(guān)系數(shù)作為邊界識別的方法：對方向總水平導數(shù)法(edge detector of total horizontal derivative，EDT)[22]歸一化可得到歸一化的方向總水平導數(shù)法(normalized edge detector of total horizontal derivative，NEDT)[3]，計算EDT與NEDT的相關(guān)系數(shù)R，用此相關(guān)系數(shù)作為邊界識別結(jié)果。公式如下：

(5)

其中:cov(EDT,NEDT)表示兩種方法的協(xié)方差;D(EDT)與D(NEDT)分別表示兩種方法的方差，故式(5)可改寫為：

(6)

其中:EDTi與NEDTi分別表示窗口內(nèi)第i個EDT與NEDT的值;N表示窗口內(nèi)的測點數(shù);R為EDT與NEDT結(jié)果的相關(guān)程度，無量綱，R∈[-1,1]，R越接近-1，該處為地質(zhì)體的邊界的可能性越大，理想情況地質(zhì)體邊界處皆為極小值-1。本文稱上述方法為基于方向總水平導數(shù)相關(guān)系數(shù)的邊界識別方法(edge recognition based on correlation coefficient ofEDT，ERCC)。

1.3 基于相關(guān)系數(shù)邊界識別約束的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積

為降低重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算結(jié)果的發(fā)散性，可應(yīng)用傳統(tǒng)的水平梯度濾波、聚散度準則等方法[22]進行結(jié)果篩選，消除發(fā)散解。為了進一步提高篩選的精度，本文提出基于相關(guān)系數(shù)邊界識別約束的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積，首先通過ERCC預先劃定地質(zhì)體水平分布范圍，然后使用重力梯度數(shù)據(jù)進行聯(lián)合歐拉反褶積計算得到初步結(jié)果，之后刪除劃定范圍外的發(fā)散解，最后聯(lián)合傳統(tǒng)篩選方法獲取最終的可靠解。

除ERCC外，軟件中內(nèi)置的邊界識別方法有：傾斜角法[23]、解析信號法[24]、總水平導數(shù)法[25]、歸一化方向總水平導數(shù)法[3,26]和改進的方向總水平導數(shù)法[27]、方向總水平導數(shù)法[28]。使用模型試驗來驗證內(nèi)置的邊界識別方法的準確性，立方體模型在x軸與y軸水平方向分布范圍為-400～400 m，埋深為200 m。邊界識別方法效果如圖1所示。

a—總水平導數(shù)法；b—方向總水平導數(shù)法；c—解析信號法；d—改進的方向總水平導數(shù)法；e—歸一化方向總水平導數(shù)法；f—傾斜角法a—THDR; b—EDT; c—ASM; d—EEDT; e—NEDT; f—TILT圖1 邊界識別結(jié)果可視化(圖中紅色線框表示立方體邊界)Fig.1 Visualization of edge detection results(the red wireframe in the figure indicates the boundary of the cube)

2 軟件開發(fā)與功能

2.1 軟件系統(tǒng)開發(fā)

軟件是基于Windows10系統(tǒng)，其硬件配置為Intel Core i7-9750H處理器、8GB內(nèi)存，利用Python語言的PyQt5.15.4、Numpy1.20.2和Matplotlib3.4.1等函數(shù)庫開發(fā)的，軟件結(jié)構(gòu)圖如圖2所示。利用NumPy函數(shù)庫實現(xiàn)邊界識別與重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積的計算；利用Matplotlib函數(shù)庫實現(xiàn)可視化功能，包括觀測數(shù)據(jù)二維可視化、邊界識別結(jié)果可視化、歐拉反褶積結(jié)果可視化等；軟件界面與文件管理是通過PyQt5函數(shù)庫實現(xiàn)的(主界面見圖3)，通過PyQt5還實現(xiàn)了重力梯度數(shù)據(jù)文件的選擇、預覽、輸入、修改與儲存等功能。

圖2 軟件結(jié)構(gòu)Fig.2 Software architecture diagram

圖3 軟件主頁面Fig.3 Main interface of software

在二次開發(fā)方面，為了方便算法的改進與軟件功能的拓展，軟件使用“信號槽”機制將界面與計算模塊、可視化模塊相連接，這種松散的耦合機制為二次開發(fā)提供便利條件?？衫谩靶盘柌邸睓C制將其他位場數(shù)據(jù)處理解釋方法加載到軟件中，增加軟件的功能；也可將后續(xù)改進的算法直接寫入軟件計算模塊替換當前算法；同時，由于對重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積算法、邊界識別等計算模塊進行了封裝，可將其直接用于其他軟件平臺的后續(xù)開發(fā)。二次開發(fā)時只需重點關(guān)注算法代碼的編寫，降低了軟件開發(fā)難度。

2.2 軟件系統(tǒng)設(shè)計

軟件系統(tǒng)以重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算為主要功能，使用直觀的操作界面，得到可視化結(jié)果?；谶吔缱R別結(jié)果約束的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積方法的算法配準如下(見圖4)：①邊界識別，劃定地質(zhì)體水平邊界；②重力梯度聯(lián)合歐拉反褶積計算，得到初步的計算結(jié)果；③最后對初步結(jié)果進行篩選，得出最終結(jié)果。根據(jù)算法流程和使用要求，軟件系統(tǒng)應(yīng)設(shè)計如下功能：①數(shù)據(jù)讀取與管理功能；②邊界識別功能；③重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算功能；④初步計算結(jié)果篩選功能；⑤結(jié)果可視化功能。

圖4 算法流程Fig.4 Flow chart of algorithm

1)數(shù)據(jù)讀取與管理功能：對數(shù)據(jù)進行操作時可能出現(xiàn)讀取文件錯誤、數(shù)據(jù)可操作性低、數(shù)據(jù)不易查看的問題。為了降低誤操作的可能性，設(shè)計數(shù)據(jù)預覽窗口(圖5)，可在發(fā)現(xiàn)讀取錯誤文件后及時退回文件選取窗口，重新選擇數(shù)據(jù)文件。為了增加對讀入數(shù)據(jù)的可操作性，在數(shù)據(jù)管理可視化窗口可以實現(xiàn)復制、粘貼、修改等操作(圖6)，文件導出功能可將修改后的數(shù)據(jù)以文本文檔形式保存到計算機中(圖7)。為了方便查看數(shù)據(jù)，根據(jù)重力梯度數(shù)據(jù)的獲取方式，數(shù)據(jù)管理可視化窗口設(shè)計了分線功能(圖8)，以平行x軸方向作為測線布置方向，左上角下拉菜單可選擇想要查看的測線。

圖5 文件預覽窗口Fig.5 File preview window

圖6 數(shù)據(jù)操作窗口Fig.6 Data operation window

圖7 文件儲存路徑選擇窗口Fig.7 Window of file storage path selection

圖8 測線切換窗口Fig.8 Line switching window

2)邊界識別功能：引入邊界識別方法并設(shè)置為聯(lián)合歐拉反褶積計算的一個步驟，可以對地質(zhì)體的平面分布具有初步的認識，為后續(xù)的計算結(jié)果篩選提供邊界依據(jù)，減少誤差。系統(tǒng)設(shè)置了包括總水平導數(shù)法、解析信號法、傾斜角法、方向總水平導數(shù)法、歸一化的方向總水平導數(shù)法、改進的方向總水平導數(shù)法、傾斜角方向總水平導數(shù)法、相關(guān)系數(shù)法等邊界識別方法，可以依據(jù)識別效果擇優(yōu)選用。篩選計算時只保留邊界內(nèi)的數(shù)據(jù)，降低后續(xù)篩選工作的計算量，提高計算效率。邊界識別功能除了可以劃定地質(zhì)體邊界范圍外，也可將可視化結(jié)果單獨以圖片形式輸出(圖9)。

圖9 軟件輸出的邊界識別結(jié)果Fig.9 Edge detection result diagram from software

3)重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算功能：計算前要先確定地質(zhì)體水平分布范圍，為了方便對照邊界識別結(jié)果，設(shè)計了如下的參數(shù)輸入界面(圖10)。只需要輸入x與y軸方向的起止坐標作參數(shù)，即可劃出形狀為矩形的邊界范圍，簡化了操作步驟。確定邊界參數(shù)后，矩形線框會在邊界識別可視化中顯示(圖11)。界面設(shè)置為交互式界面，邊界參數(shù)可返回上一步修改，并實時更新圖中線框位置，提高了便捷性。在輸入計算窗尺寸與選擇聯(lián)合歐拉反褶積計算所需的梯度數(shù)據(jù)后，即可開始計算。參與計算的3個梯度數(shù)據(jù)分量可通過界面中(圖11)的下拉框選擇，軟件依據(jù)式(4)預先設(shè)定選擇的梯度數(shù)據(jù)分量為Vxz、Vyz與Vzz。軟件的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算公式與周文月等[12]提出的3個不同高度數(shù)據(jù)聯(lián)合歐拉求解矩陣表達式格式一致，只是選取的數(shù)據(jù)類型不同，所以軟件不僅能聯(lián)合重力梯度數(shù)據(jù)進行聯(lián)合歐拉反褶積計算，也可聯(lián)合3種不同高度的重力異常數(shù)據(jù)進行聯(lián)合歐拉反褶積計算。

圖10 劃定地質(zhì)體范圍輸入界面Fig.10 Delimit the geological body range input interface

圖11 重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算參數(shù)輸入界面Fig.11 Joint Euler deconvolution of multiple gravity gradiometry tensors calculation parameters input interface

4)結(jié)果篩選功能：當聯(lián)合歐拉反褶積計算后的結(jié)果(圖12)與地質(zhì)體實際情況不相符時，需要通過有效的篩選方法篩選計算結(jié)果，使結(jié)果更加準確地描繪地質(zhì)體位置與形狀。軟件內(nèi)置的篩選方法有水平梯度濾波法[22]、聚散度準則[22]等(圖13)。主要包括：①水平梯度濾波法：計算窗口內(nèi)區(qū)域的水平梯度的模，并與用戶輸入的系數(shù)相乘得到篩選標準，刪除水平梯度的模小于篩選標準的結(jié)果。用戶輸入的3個系數(shù)分別對應(yīng)著先前重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算時選取的重力梯度分量。②主體異常距離準則：根據(jù)重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算的原理，計算結(jié)果的水平位置應(yīng)該位于當前計算窗口范圍內(nèi)。對于超出計算窗口范圍的點坐標，認定該計算結(jié)果誤差較大，予以刪除。③聚散度準則：根據(jù)重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積在地質(zhì)體上方相鄰多個滑動窗口所得到的計算結(jié)果相關(guān)性較強的特點，正確的計算結(jié)果應(yīng)該比較密集的出現(xiàn)。過于分散的結(jié)果視為誤差較大的結(jié)果，予以刪除。用戶可在篩選標準輸入界面輸入聚散度準則所需的作用半徑與聚散度指數(shù)。④邊界識別：重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算時依據(jù)邊界識別結(jié)果可視化劃定的范圍作為篩選邊界，刪除邊界外的計算結(jié)果，避免誤差和屏蔽范圍外的無效結(jié)果。篩選后的結(jié)果如圖14所示。

圖12 重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積初步計算結(jié)果Fig.12 Preliminary calculation results of joint Euler deconvolution of multiple gravity gradiometry tensors

圖13 結(jié)果篩選標準輸入界面Fig.13 Results filtering standard input interface

圖14 重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算結(jié)果可視化Fig.14 Visualization of joint Euler deconvolution of multiple gravity gradiometry tensors results

5)可視化：為了可以直觀地觀察梯度數(shù)據(jù)的情況和邊界識別結(jié)果，利用等勢線圖能夠清晰地描繪數(shù)值的分布情況(圖15)。圖15顯示的是x軸與y軸方向分布為-400～400 m、埋深為200 m的單立方體模型的重力異常數(shù)據(jù)與重力梯度數(shù)據(jù)等值線圖，圖中紅色線框表示立方體模型的邊界。在立方體邊界處重力梯度數(shù)據(jù)的數(shù)值變化明顯，重力異常數(shù)據(jù)可視化與模型實際位置相符。

a—Vxx分量；b—Vxy分量；c—Vxz分量；d—Vyy分量；e—Vyz分量；f—Vzz分量；g—Vz圖15 重力異常數(shù)據(jù)與重力梯度數(shù)據(jù)等值線(圖中紅色線框表示立方體邊界)Fig.15 Contour figure of gravity data and gravity gradient data(the red wireframe in the figure indicates the boundary of the cube)

利用三維散點的形式表現(xiàn)重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算結(jié)果的空間位置關(guān)系，簡單且高效地將與x、y、z坐標相關(guān)的空間位置轉(zhuǎn)化為易于觀察的分布在三維空間的散點。同時，邊界識別的結(jié)果也可置于三維散點圖的上方(圖14)，使得結(jié)果展示得更為直觀、全面。為了更準確地觀察計算結(jié)果與正演模型相對位置關(guān)系，判斷算法的分辨率，可視化模塊中還具有添加輔助線框功能。軟件可根據(jù)用戶輸入的坐標數(shù)據(jù)(圖16)繪制立方體線框并顯示在計算結(jié)果三維散點可視化中，輔助用戶更好地對比正演模型與計算結(jié)果之間的區(qū)別，或輔助劃定實測數(shù)據(jù)解釋的地質(zhì)體范圍。

圖16 輔助線框坐標輸入界面Fig.16 Auxiliary wireframe coordinate input interface

3 數(shù)據(jù)試驗

3.1理論模型數(shù)據(jù)試驗

為了驗證ERCC對不同深度地質(zhì)體邊界的探測精度，使用加入均值為0，標準差為0.01 mGal的高斯隨機噪聲的立方體模型數(shù)據(jù)進行試驗。取x和y軸方向范圍均為0～2 000 m、深度方向范圍為0～1 000 m的空間，存在3個立方體模型，具體分布情況見表1。邊界識別結(jié)果可視化如圖17所示。

表1 立方體模型參數(shù)表Toble 1 Cube model parameters

圖17 邊界識別結(jié)果可視化(圖中紅色線框表示立方體邊界)Fig.17 Visualization of edge detection(the red wireframe in the figure indicates the boundary of the cube)

地質(zhì)體頂面埋深分別為100 m(A)、400 m(B)和700 m(C)，由圖17可知，隨著地質(zhì)體埋深的增加，相關(guān)系數(shù)變得越來越大。無地質(zhì)體邊界分布的區(qū)域相關(guān)系數(shù)范圍為0.75～1，3個地質(zhì)體邊界的相關(guān)系數(shù)范圍由A至C依次為：-0.75～-0.25、-0.25～0.25與0～0.5，埋深淺的地質(zhì)體的相關(guān)系數(shù)小于深的，但都與0.75有一定差距。根據(jù)相關(guān)系數(shù)的極小值表示地質(zhì)體邊界的特性，埋深最淺的地質(zhì)體A劃定的范圍為x軸向1200～1 850 m，y軸向150～850 m，最深的地質(zhì)體C劃定的x與y軸向的范圍都是200～800 m，可見ERCC可以準確地識別出同一區(qū)域內(nèi)不同深度的地質(zhì)體的邊界。同時ERCC也良好地壓制了噪聲的影響，結(jié)果中并未產(chǎn)生與實際情況不符的邊界。

為驗證基于相關(guān)系數(shù)邊界識別約束的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積的準確性和開發(fā)的軟件的有效性，在笛卡爾坐標系下建立理論模型并進行模型試驗。選取x和y軸方向范圍均為-1 000～1 000 m、深度方向范圍為0～1 000 m的空間，設(shè)地下有一個立方體形狀地質(zhì)體，長和寬為800 m，高為200 m，頂面埋深為200 m，測線間距為20 m，測線上測點間距為20 m，地質(zhì)體剩余密度為1.0 g/cm3，并加入均值為0，標準差為0.01 mGal的高斯隨機噪聲。聯(lián)合歐拉反褶積計算滑動窗口為19×19，水平梯度濾波法系數(shù)設(shè)為1，聚散度指數(shù)設(shè)為5、半徑設(shè)為1.5倍測線間距。計算結(jié)果如圖18所示，圖中黑色線框為添加的輔助線框，表示模型所在位置。

圖18 篩選后的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積計算結(jié)果Fig.18 Screened results of joint Euler deconvolution of multiple gravity gradiometry tensors

由圖18可知，篩選后的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積的計算結(jié)果主要分布在200～280 m的深度范圍，水平方向投影在立方體邊界的位置上，在立方體頂點附近聚集著深度大于280 m的解。計算結(jié)果總體上分布在地質(zhì)體的邊界處，指示出地質(zhì)體的分布范圍與形狀。通過模型數(shù)據(jù)試驗，證明計算結(jié)果可以準確地標示出地質(zhì)體的水平位置與埋深，驗證了該方法的準確性與有效性。

3.2 實測數(shù)據(jù)試驗

為了進一步驗證方法和軟件的適用性與準確性，采用文頓鹽丘(Vinton Dome)地區(qū)的實測重力梯度數(shù)據(jù)進行試驗。文頓鹽丘位于美國路易斯安那州西南部與得克薩斯州交界處，地層以沉積巖為主[29]。該地區(qū)大陸架盆地的形成時間為新生代，地層以侵入巖為主。鹽丘巖性主要為頁巖和砂巖，剩余密度為2.2 g/cm3，蓋巖的成分為石膏和硬石膏，該區(qū)域的重力梯度異常主要是其蓋巖引起的。實測數(shù)據(jù)位于WGS84坐標系下，x軸代表EW方向，范圍為440.5～444.5 km；y軸代表SN方向，范圍為3 332.8～3 336.8 km。測區(qū)平均分布100條測線，每條測線平均分布100個測點。選取實測數(shù)據(jù)中Vxz、Vyz和Vzz這3個分量進行聯(lián)合歐拉反褶積計算,計算滑動窗口為11×11，水平梯度濾波法系數(shù)設(shè)為1.5，聚散度指數(shù)設(shè)為6、半徑設(shè)為1.5倍測線間距,結(jié)果見圖19。

a—三維顯示；b—z方向視圖a—3D display; b—z-direction view圖19 文頓鹽丘數(shù)據(jù)聯(lián)合歐拉反褶積計算結(jié)果Fig.19 Calculation results of joint Euler deconvolution of multiple gravity gradiometry tensors of Vinton Dome data

從圖19a可以看出蓋巖的埋深約為160～550 m，西側(cè)蓋巖埋藏較淺，埋深約為160 m；南側(cè)蓋巖埋藏較深，埋深約為550 m。從圖19b可以看出巖蓋整體呈現(xiàn)弧形，x與y軸向上分布范圍分別為44 1700～442 900 m和333 380 0～333 480 0 m，蓋巖的埋深由東南側(cè)向東側(cè)、南側(cè)逐漸變淺。綜上，軟件計算結(jié)果與其他學者研究結(jié)果[30-31]一致，驗證了軟件計算的準確性。

4 結(jié)論

本文提出了基于相關(guān)系數(shù)邊界識別約束的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積，并開發(fā)了一種基于Python及其函數(shù)庫的重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積可視化軟件。提出的約束方法可以有效地減少歐拉反褶積的發(fā)散解，提高準確性；可視化軟件具有數(shù)據(jù)/文件管理、二/三維可視化、邊界識別、重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積等功能。經(jīng)模型試驗和實測數(shù)據(jù)試驗，驗證了計算的準確性和軟件的有效性，且軟件操作方法較為簡單，可視化效果較好，為算法研究與應(yīng)用提供了可靠的平臺。同時，重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積算法已被模塊化封裝，方便后續(xù)研究進行二次開發(fā)。

致謝：感謝Bell Geospace Inc.提供文頓鹽丘地區(qū)的實測數(shù)據(jù)。