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

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

0 引言

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

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

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

可見，使用具有豐富函數(shù)庫(kù)、可移植性高的語(yǔ)言，開發(fā)易用的可視化軟件是十分必要的。這有利于降低方法二次開發(fā)與軟件移植的難度，使研究人員減少對(duì)數(shù)據(jù)管理與可視化的開發(fā)成本，擴(kuò)展解釋方法的適用范圍。因此，本文基于簡(jiǎn)潔實(shí)用的原則，利用Python語(yǔ)言及其函數(shù)庫(kù)開發(fā)了一個(gè)支持?jǐn)?shù)據(jù)/文件管理、二/三維可視化、邊界識(shí)別、重力梯度數(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)分別為笛卡爾坐標(biāo)系中觀測(cè)點(diǎn)和地質(zhì)體中心坐標(biāo)；f和B分別表示觀測(cè)異常場(chǎng)和背景場(chǎng)；N為構(gòu)造指數(shù)。將式(1)轉(zhuǎn)化為向量積形式，即：

(2)

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

(3)

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

(4)

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

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

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

(5)

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

(6)

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

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

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

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

a—總水平導(dǎo)數(shù)法；b—方向總水平導(dǎo)數(shù)法；c—解析信號(hào)法；d—改進(jìn)的方向總水平導(dǎo)數(shù)法；e—?dú)w一化方向總水平導(dǎo)數(shù)法；f—傾斜角法a—THDR; b—EDT; c—ASM; d—EEDT; e—NEDT; f—TILT圖1 邊界識(shí)別結(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語(yǔ)言的PyQt5.15.4、Numpy1.20.2和Matplotlib3.4.1等函數(shù)庫(kù)開發(fā)的，軟件結(jié)構(gòu)圖如圖2所示。利用NumPy函數(shù)庫(kù)實(shí)現(xiàn)邊界識(shí)別與重力梯度數(shù)據(jù)聯(lián)合歐拉反褶積的計(jì)算；利用Matplotlib函數(shù)庫(kù)實(shí)現(xiàn)可視化功能，包括觀測(cè)數(shù)據(jù)二維可視化、邊界識(shí)別結(jié)果可視化、歐拉反褶積結(jié)果可視化等；軟件界面與文件管理是通過PyQt5函數(shù)庫(kù)實(shí)現(xiàn)的(主界面見圖3)，通過PyQt5還實(shí)現(xiàn)了重力梯度數(shù)據(jù)文件的選擇、預(yù)覽、輸入、修改與儲(chǔ)存等功能。

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

圖3 軟件主頁(yè)面Fig.3 Main interface of software

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

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

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

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

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

圖5 文件預(yù)覽窗口Fig.5 File preview window

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

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

圖8 測(cè)線切換窗口Fig.8 Line switching window

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

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

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

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

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

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

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

圖13 結(jié)果篩選標(biāo)準(zhǔn)輸入界面Fig.13 Results filtering standard input interface

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

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

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)

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

圖16 輔助線框坐標(biāo)輸入界面Fig.16 Auxiliary wireframe coordinate input interface

3 數(shù)據(jù)試驗(yàn)

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

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

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

圖17 邊界識(shí)別結(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個(gè)地質(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可以準(zhǔn)確地識(shí)別出同一區(qū)域內(nèi)不同深度的地質(zhì)體的邊界。同時(shí)ERCC也良好地壓制了噪聲的影響，結(jié)果中并未產(chǎn)生與實(shí)際情況不符的邊界。

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

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

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

3.2 實(shí)測(cè)數(shù)據(jù)試驗(yàn)

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

a—三維顯示；b—z方向視圖a—3D display; b—z-direction view圖19 文頓鹽丘數(shù)據(jù)聯(lián)合歐拉反褶積計(jì)算結(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è)逐漸變淺。綜上，軟件計(jì)算結(jié)果與其他學(xué)者研究結(jié)果[30-31]一致，驗(yàn)證了軟件計(jì)算的準(zhǔn)確性。

4 結(jié)論

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

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