李忠平，戴廣凱，張茂輝

(1.中國(guó)地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢430074；2.中國(guó)冶金地質(zhì)總局 山東正元地質(zhì)勘查院，山東 濟(jì)南 250014； 3.山東省地質(zhì)調(diào)查院，山東 濟(jì)南 250013； 4.中國(guó)海洋大學(xué) 工程學(xué)院，山東 青島 266100)

0 引言

1∶5萬(wàn)重力遠(yuǎn)區(qū)地改原則上使用1∶5萬(wàn)衛(wèi)星DEM數(shù)據(jù)(節(jié)點(diǎn)距25 m)，也可采用RGIS自帶高程庫(kù)進(jìn)行遠(yuǎn)Ⅰ、Ⅱ區(qū)地改。由于RGIS自帶高程庫(kù)是基于“區(qū)域重力信息系統(tǒng)2006”自帶1 km×1 km節(jié)點(diǎn)高程，網(wǎng)格間距偏大[1-3]。本次對(duì)山東省棲霞市臧家莊幅1∶5萬(wàn)重力觀測(cè)數(shù)據(jù)采用以上兩種網(wǎng)格節(jié)點(diǎn)的高程數(shù)據(jù)進(jìn)行了遠(yuǎn)Ⅰ、Ⅱ區(qū)地改試驗(yàn)，以評(píng)價(jià)兩種方法的合理性和可行性。

1∶5萬(wàn)重力勘查跨越了不同的DEM數(shù)據(jù)圖幅，要獲取有重力勘查邊界的DEM數(shù)據(jù)區(qū)域時(shí)，必須對(duì)多個(gè)圖幅的DEM數(shù)據(jù)進(jìn)行拼接。由于數(shù)據(jù)拼接時(shí)插值方法不同，其數(shù)據(jù)拼接精度不同。為滿足重力遠(yuǎn)Ⅰ區(qū)地改的需要，采用Suffer8雙線性插值法快速進(jìn)行DEM數(shù)據(jù)拼接，解決了該區(qū)地形高差較大時(shí)插值效果差等問題[4-6]，從而使高精度GPS因信號(hào)較差而采集的個(gè)別低精度高程點(diǎn)對(duì)應(yīng)地改誤差減小[7-9]，同時(shí)使用Suffer8可更加方便地自定義不同的網(wǎng)格間距進(jìn)行DEM數(shù)據(jù)拼接，以滿足不同的探測(cè)需求，解決了地形對(duì)重力觀測(cè)異常的影響[10-12]。

試驗(yàn)使用Suffer8進(jìn)行DEM數(shù)據(jù)網(wǎng)格拼接，然后用GeoIPAS3.0軟件進(jìn)行遠(yuǎn)Ⅰ區(qū)地形改正時(shí)，由于DEM數(shù)據(jù)網(wǎng)格間距不同，1∶5萬(wàn)重力遠(yuǎn)Ⅰ區(qū)(2～20 km)地改的速度不同，網(wǎng)格間距大，地形改正效率高[13-14]，對(duì)應(yīng)地改均方誤差大。實(shí)際工作時(shí)，采用RGIS自帶高程庫(kù)進(jìn)行1∶5萬(wàn)重力遠(yuǎn)Ⅰ、Ⅱ區(qū)地改速度快，760個(gè)測(cè)點(diǎn)1～2 min即可完成。而采用GeoIPAS3.0軟件進(jìn)行遠(yuǎn)Ⅰ、Ⅱ區(qū)地改時(shí)760個(gè)測(cè)點(diǎn)一次完成要40 h(網(wǎng)格間距25 m)，數(shù)據(jù)處理效率較低。網(wǎng)格間距成為制約地改效率和精度、地質(zhì)解釋效果的主要因素[15-16]。為探索最佳網(wǎng)格間距，本文通過對(duì)山東省棲霞市臧家莊幅1∶5萬(wàn)重力遠(yuǎn)Ⅰ區(qū)地形改正數(shù)據(jù)分析，采用25、50、100、200 m四種網(wǎng)格節(jié)點(diǎn)距進(jìn)行了1∶5萬(wàn)DEM數(shù)據(jù)拼接和對(duì)應(yīng)重力遠(yuǎn)Ⅰ區(qū)地改、均方誤差計(jì)算。將1∶5萬(wàn)衛(wèi)星DEM數(shù)據(jù)節(jié)點(diǎn)距擴(kuò)大一倍，即用50 m的網(wǎng)格距拼接，能較好滿足地形改正效率與精度的要求，縮短了1∶5萬(wàn)重力遠(yuǎn)區(qū)地形改正時(shí)間，提高了資料整理效率，更有利于推斷解釋工作[17-19]。

將重力勘查數(shù)據(jù)與地質(zhì)、遙感資料綜合是發(fā)現(xiàn)深部隱伏斷裂的重要手段[20]。而遠(yuǎn)區(qū)地改精度在重力地改中占有重要地位，它影響對(duì)地質(zhì)和遙感資料的準(zhǔn)確解釋，對(duì)發(fā)現(xiàn)深大斷裂和深部礦體意義重大[21]。由不同網(wǎng)格節(jié)點(diǎn)的DEM數(shù)據(jù)和RGIS自帶高程庫(kù)做地改后的精度對(duì)比，研究其對(duì)布格重力異常的影響，很容易被技術(shù)人員忽略。它會(huì)影響對(duì)局部重力異常的解釋，有時(shí)會(huì)漏掉一些有意義的異常，從而錯(cuò)失找礦機(jī)會(huì)。

因此，在1∶5萬(wàn)重力勘查區(qū)地形較為復(fù)雜、布格重力異?？偩燃暗匦胃恼秸`差要求較高時(shí)，首選1∶5萬(wàn)衛(wèi)星DEM數(shù)據(jù)(節(jié)點(diǎn)距50 m)進(jìn)行重力遠(yuǎn)區(qū)地改。在滿足規(guī)范要求的情況下，也可采用RGIS自帶高程庫(kù)進(jìn)行重力遠(yuǎn)區(qū)地改。經(jīng)研究對(duì)比兩種遠(yuǎn)區(qū)地改方法各有優(yōu)缺點(diǎn)。對(duì)于大測(cè)區(qū)大數(shù)據(jù)量的重力遠(yuǎn)區(qū)地改，建議分區(qū)分幅進(jìn)行，同時(shí)重視DEM數(shù)據(jù)拼接時(shí)網(wǎng)格間距的選擇，避免重力遠(yuǎn)區(qū)地改耗時(shí)太長(zhǎng)，降低工作效率。

1 重力遠(yuǎn)區(qū)地改方法

將以測(cè)點(diǎn)為中心的四周地形分割成許多小塊，盡可能模擬實(shí)際地形。計(jì)算出每一小塊地形質(zhì)量對(duì)測(cè)點(diǎn)的重力值，再累加求和得出該點(diǎn)地形改正值。共用點(diǎn)法為其常用方法之一，共用點(diǎn)法是將實(shí)際地形分割成四棱柱體。首先計(jì)算重力測(cè)點(diǎn)(自由網(wǎng)點(diǎn))附近4個(gè)節(jié)點(diǎn)的地形改正值(計(jì)算時(shí)用測(cè)點(diǎn)高程值代替四個(gè)節(jié)點(diǎn)的高程值)，然后再將4個(gè)節(jié)點(diǎn)的地形改正值內(nèi)插到測(cè)點(diǎn)位置上作為測(cè)點(diǎn)地形改正值，計(jì)算公式為(1)[22]

(1)

式中：δg(x,y)為測(cè)點(diǎn)地形改正值；f為萬(wàn)有引力常數(shù)6.67×10-11m3/(kg·s2)；σ為地形校正密度；hij為節(jié)點(diǎn)與計(jì)算點(diǎn)的高程差；rij為節(jié)點(diǎn)與計(jì)算點(diǎn)的距離；Δx、Δy為網(wǎng)格距或節(jié)點(diǎn)距；Cij為積分常數(shù)，此處選用梯形系數(shù)Cij：

由上式可知，測(cè)點(diǎn)重力地形改正值δg(x,y)與網(wǎng)格距或節(jié)點(diǎn)距Δx、Δy有關(guān)，網(wǎng)格距或節(jié)點(diǎn)距選擇合適與否，決定了重力遠(yuǎn)區(qū)地改的精度高低[22-23]。

計(jì)算測(cè)點(diǎn)附近的4個(gè)節(jié)點(diǎn)的測(cè)點(diǎn)重力值后，采用雙線性插值公式計(jì)算出測(cè)點(diǎn)重力數(shù)據(jù)的遠(yuǎn)區(qū)改正值(圖1)。

圖1 雙線性插值法圖示Fig.1 Bilinear Interpolation Graphics

雙線性插值計(jì)算公式如下：

(2)

(3)

(4)

式中：g(x,y)為地改值雙線性插值結(jié)果；g(xi,yi)、g(xi,yi+1)、g(xi+1,yi)、g(xi+1,yi+1)為4個(gè)節(jié)點(diǎn)的地改值；xi、xi+1、yi、yi+1為節(jié)點(diǎn)坐標(biāo)值；Δx、Δy為節(jié)點(diǎn)距。

由上式可知，測(cè)點(diǎn)重力地形改正值δg(x,y)精度及地改速度受網(wǎng)格距和邊緣方向插值算法的制約，為解決這一問題，本文提出對(duì)雙線性插值法進(jìn)行改進(jìn)。

2 改進(jìn)的雙線性插值重力遠(yuǎn)區(qū)地改計(jì)算方法

針對(duì)傳統(tǒng)的雙線性插值在實(shí)現(xiàn)沿任意邊緣方向插值時(shí)圖像邊緣模糊及算法速度、精度問題，本次研究對(duì)雙線性插值法進(jìn)行改進(jìn)，基于邊緣方向的插值是一種有效圖像插值方法。 雙線性插值沿水平和垂直兩個(gè)方向進(jìn)行插值，不利于圖像的實(shí)時(shí)放大處理，或難以達(dá)到最佳效果。本文提出首先確定插值點(diǎn)的重力值和邊緣方向，再引入雙線性插值算法，沿邊緣方向進(jìn)行插值，有效地提高了插值圖像邊緣的清晰度及算法速度、地改精度[24-25]。

如圖2所示，g(xi,yi)、g(xi,yi+1)、g(xi+1,yi)、g(xi+1,yi+1)為4個(gè)節(jié)點(diǎn)的地改值，O點(diǎn)為位于非正交直線AB上的待求內(nèi)插點(diǎn)，非正交直線4個(gè)節(jié)點(diǎn)所構(gòu)成的四邊形非垂直相交。

圖2 改進(jìn)的雙線性插值法圖示Fig.2 Improved bilinear interpolation method

由A點(diǎn)在g(xi,yi)、g(xi+1,yi)之間進(jìn)行線性插值得到A點(diǎn)的地改值公式為：

(5)

由B點(diǎn)在g(xi,yi+1)、g(xi+1,yi+1)之間進(jìn)行線性插值得到B點(diǎn)的地改值公式為：

(6)

由O點(diǎn)在A、B之間進(jìn)行線性插值得到O點(diǎn)的地改值公式為：

(7)

使用matlab編程實(shí)現(xiàn)改進(jìn)的雙線性插值法，并進(jìn)行了重力遠(yuǎn)區(qū)地改值計(jì)算。

在測(cè)區(qū)西南角選擇了4個(gè)重力觀測(cè)點(diǎn)(節(jié)點(diǎn))(圖3)，改進(jìn)的雙線性插值法4節(jié)點(diǎn)編號(hào)分別為g(xi,yi)、g(xi,yi+1)、g(xi+1,yi)、g(xi+1,yi+1)，其坐標(biāo)見表1。依據(jù)圖3，取si=2 000 m，si+1=1 000 m，α=45°。

由圖3可知，在測(cè)區(qū)西南角，Δg一般在(-15～-10)×10-5m/s2之間變化，因此改進(jìn)的雙線性插值法較好反映了重力異常的變化趨勢(shì)，可以提高推斷解釋的精度。

圖3 山東煙臺(tái)地區(qū)區(qū)域布格重力異常平面Fig.3 Plane of Bouguer gravity anomaly in the Yantai area of Shandong Province

表1 改進(jìn)的雙線性插值法試驗(yàn)計(jì)算對(duì)比

3 重力遠(yuǎn)區(qū)地改精度對(duì)比分析

1∶5萬(wàn)衛(wèi)星DEM數(shù)據(jù)原始節(jié)點(diǎn)距為25 m，基于該節(jié)點(diǎn)距的山東省棲霞市臧家莊幅1∶5萬(wàn)重力勘查DEM數(shù)據(jù)網(wǎng)格拼接用于重力遠(yuǎn)區(qū)地改存在效率較低的事實(shí)。本次選取該區(qū)27個(gè)測(cè)點(diǎn)原始與檢查數(shù)據(jù)，在使用Suffer8進(jìn)行DEM數(shù)據(jù)拼接時(shí)選取了25、50、100、200 m四種網(wǎng)格節(jié)點(diǎn)距，使用金維GeoIPAS3.0軟件將原始與檢查數(shù)據(jù)分別進(jìn)行了重力遠(yuǎn)Ⅰ區(qū)(2～20 km)地改試驗(yàn)。地改精度與時(shí)效指標(biāo)情況見表2。由表2可知，網(wǎng)格節(jié)點(diǎn)距越大，重力遠(yuǎn)Ⅰ區(qū)(2～20 km)地改時(shí)間越短，地改均方誤差越大。經(jīng)試算，該區(qū)合適的網(wǎng)格節(jié)點(diǎn)距為50 m，其重力遠(yuǎn)Ⅰ區(qū)地改均方誤差為0.000 36，優(yōu)于使用RGIS自帶高程庫(kù)重力遠(yuǎn)Ⅰ區(qū)地形改正均方誤差(0.03),滿足規(guī)范要求。

對(duì)山東省棲霞市臧家莊幅1∶5萬(wàn)重力勘查760測(cè)點(diǎn)采用25 m×25 m網(wǎng)格節(jié)點(diǎn)距DEM拼接數(shù)據(jù)。金維GeoIPAS3.0軟件進(jìn)行了重力遠(yuǎn)Ⅰ、Ⅱ區(qū)地形改正，760個(gè)測(cè)點(diǎn)一次地改需要42 h?？梢?5 m×25 m網(wǎng)格節(jié)點(diǎn)距DEM數(shù)據(jù)用于地改效率較低。

表2 不同網(wǎng)格間距重力遠(yuǎn)Ⅰ區(qū)地改精度統(tǒng)計(jì)

4 應(yīng)用效果對(duì)比

4.1 重力異常特征對(duì)比

研究區(qū)為山東棲霞臧家莊盆地，區(qū)內(nèi)巖漿巖分布廣泛，以中生代郭家?guī)X序列羅家單元斑狀中細(xì)粒含黑云二長(zhǎng)花崗巖為主。玲瓏—招風(fēng)頂脈巖帶通過該區(qū)，脈體產(chǎn)狀為：走向20°～40°，傾向NW向?yàn)橹?，SE向次之，單脈寬2～50 m，長(zhǎng)度多在1 100～8 000 m。巖石類型有煌斑巖、閃長(zhǎng)玢巖、石英閃長(zhǎng)玢巖、石英正長(zhǎng)斑巖、細(xì)晶閃長(zhǎng)巖等。臧家莊盆地西緣紫現(xiàn)頭斷裂和北緣西林?jǐn)嗔呀粎R于本區(qū)，構(gòu)造復(fù)雜，次級(jí)斷裂眾多。區(qū)內(nèi)NE向斷裂構(gòu)造較發(fā)育，且多以張性為主，與金礦關(guān)系密切。

研究區(qū)區(qū)域重力異常特征：區(qū)內(nèi)重力場(chǎng)值西北低，東南高(圖3)，變化范圍為(-14～14)×10-5m/s2，異常等值線主要呈現(xiàn)NNE向展布，并疊加一些局部重力異常。重力場(chǎng)特征為密集的重力異常梯級(jí)帶。根據(jù)區(qū)域地質(zhì)資料推斷為典型的深大斷裂的反映。在其西北部區(qū)域，重力梯級(jí)帶變得寬緩，局部區(qū)段扭動(dòng)頻率增強(qiáng)，說(shuō)明局部構(gòu)造發(fā)育，但重力場(chǎng)構(gòu)造線總體為NNE，與地質(zhì)構(gòu)造格架吻合。

研究區(qū)地形特點(diǎn)：在研究區(qū)西部和中東部均有2處相對(duì)較高地形(圖4)，高程一般在111～316 m。全區(qū)最高點(diǎn)在東部，高程360 m；最低點(diǎn)在東北角，高程70 m，全區(qū)地形高差為290 m。這些相對(duì)較高地形為山區(qū)，中部和東部共有3條帶狀相對(duì)低值地形區(qū)，為山溝和第四系。

研究區(qū)已知礦床地質(zhì)特征：該區(qū)內(nèi)已知的礦產(chǎn)主要是金礦，其中小型礦床1處，礦點(diǎn)2處，礦化點(diǎn)1處。郭家店小型金礦，賦存于NE向張性破碎帶中，產(chǎn)于一組NE向的石英脈中。成礦時(shí)代為燕山期，圍巖為斑狀中?；◢忛W長(zhǎng)巖，主要蝕變?yōu)辄S鐵礦化、絹云母化、褐鐵礦化等，化學(xué)分析金含量在(3～10)×10-6之間(圖5)。

研究區(qū)重力異常特征：重力異常總體走向?yàn)镹E向。重力場(chǎng)表現(xiàn)為東高西低特征，局部稍有變化。低值區(qū)位于工作區(qū)西北角，布格重力異常值一般在(-14～0)×10-5m/s2之間變化，與二長(zhǎng)花崗巖和古元古界粉子山群祝家夼組透閃大理巖和變粒巖對(duì)應(yīng)；高值區(qū)位于工作區(qū)東北部，布格重力異常值一般在(0～15)×10-5m/s2之間變化，與二長(zhǎng)花崗巖和石英閃長(zhǎng)玢巖脈對(duì)應(yīng)(圖5、圖6)。

研究區(qū)已知金礦脈的重力異常特征：西北部3條和西南、中南部9條金礦脈均處在布格異常等值線梯級(jí)帶和扭曲變形之處，布格重力異常值一般在(-13～5)×10-5m/s2之間變化，屬相對(duì)低緩布格重力異常(圖1、圖6)。

圖6與圖7是分別基于RGIS自帶高程庫(kù)(網(wǎng)格距1 km)和基于DEM(網(wǎng)格距50 m)改進(jìn)的雙線性插值中遠(yuǎn)區(qū)地改后布格重力異常分布特征。二者布格重力異常展布形態(tài)基本相同，只在基于RGIS自帶高程庫(kù)中遠(yuǎn)區(qū)地改后布格重力異常圖的北東部有兩處重力高值異常與基于DEM改進(jìn)的雙線性插值中遠(yuǎn)區(qū)地改后布格重力異常稍有不同，在圖7中兩處重力高值異常消失，顯示為正重力異常梯級(jí)帶，提高了地質(zhì)解釋的精度。產(chǎn)生這一差別的原因與這兩處高差變化之大有關(guān)， 是RGIS自帶高程庫(kù)的網(wǎng)格距遠(yuǎn)大于DEM網(wǎng)格距所致。

圖4 研究區(qū)地形Fig.4 Topographic map of the study area

1—推測(cè)斷裂; 2—壓扭性斷裂; 3—張性斷裂; 4—韌性斷層; 5—構(gòu)造破碎帶; 6—金礦脈; 7—鉆孔1—speculated fracture; 2—compression torsion fracture; 3—tensile fracture; 4—ductile fault; 5—structural fracture zone; 6—gold vein vein; 7—drilling圖5 研究區(qū)構(gòu)造及礦點(diǎn)分布Fig.5 Structure of the study area and Titanium Vein

圖6 基于RGIS自帶高程庫(kù)中遠(yuǎn)區(qū)地改后布格重力異常Fig.6 Based on the Bouguer gravity anomaly map of the remote area of RGIS with its own elevation

圖7 基于DEM改進(jìn)的雙線性插值中遠(yuǎn)區(qū)地改后布格重力異常Fig.7 Modified Bouguer gravity anomaly map of far area by bilinear interpolation based on DEM

研究區(qū)構(gòu)造特征：布格重力高異常等值線密集、梯度變化較大及錯(cuò)斷部位，是推斷深大斷裂的依據(jù)。區(qū)內(nèi)具多個(gè)圈閉中心的異常，表明局部構(gòu)造發(fā)育。位于研究區(qū)西部摩天嶺—南交毛寨—艾山湯一帶斷裂，異常呈梯級(jí)帶和扭曲收縮特征，沿NE向展布，布格重力異常值在-11×10-5m/s2左右；位于研究區(qū)西部高家—盛家—疃頂一帶斷裂，異常呈梯級(jí)帶和扭曲收縮特征，沿NE向展布，布格重力異常值在-10×10-5m/s2左右(圖5、圖7)。

4.2 重力異常相關(guān)分析

研究區(qū)的地形變化對(duì)重力異常的影響不容忽視，無(wú)論地形的隆起或是凹陷，都將減小異常的實(shí)測(cè)值導(dǎo)致重力異常圖在形態(tài)上的直接變化，影響對(duì)異常的客觀解釋。將布格重力異常和地形兩個(gè)物理量進(jìn)行相關(guān)分析，以反映研究區(qū)內(nèi)地形起伏對(duì)重力異常的影響情況，以及地形改正是否完全或改算方法是否完善。重力異常和地形之間的相關(guān)分析依據(jù)式(8)[27]：

(8)

圖8、圖9分別為基于RGIS自帶高程庫(kù)(網(wǎng)格距1 km)和DEM(網(wǎng)格距50 m)改進(jìn)的雙線性插值中遠(yuǎn)區(qū)地改后布格重力異常和高程相關(guān)系數(shù)分布特征。圖8的北東部有兩處相關(guān)系數(shù)與圖9的不同，在圖9中這兩處相關(guān)系數(shù)明顯增大，并且連續(xù)。產(chǎn)生這一差別的原因亦與該處高差變化之大有關(guān)，是RGIS自帶高程庫(kù)的網(wǎng)格距遠(yuǎn)大于DEM網(wǎng)格距所致，說(shuō)明基于DEM改進(jìn)的雙線性插值遠(yuǎn)區(qū)地改完全度高于基于RGIS自帶高程庫(kù)遠(yuǎn)區(qū)地改，并且改算方法精度較高。

5 結(jié)論與討論

本次試驗(yàn)重力遠(yuǎn)區(qū)地改精度以均方誤差體現(xiàn)。經(jīng)對(duì)比分析：基于1∶5萬(wàn)衛(wèi)星DEM數(shù)據(jù)50 m×50 m網(wǎng)格節(jié)點(diǎn)距的1∶5萬(wàn)重力遠(yuǎn)區(qū)地改精度優(yōu)于使用RGIS自帶高程庫(kù)(網(wǎng)格距1 km)重力遠(yuǎn)Ⅰ區(qū)地形改正精度。在不同地形和不同精度要求的測(cè)區(qū)可選擇使用以上兩種方法。該試驗(yàn)證明：對(duì)重力遠(yuǎn)區(qū)地改不能完全以DEM數(shù)據(jù)中默認(rèn)網(wǎng)格節(jié)點(diǎn)距為準(zhǔn)，要通過DEM數(shù)據(jù)不同網(wǎng)格節(jié)點(diǎn)距試算重力遠(yuǎn)區(qū)地改值。改進(jìn)的雙線性插值法計(jì)算遠(yuǎn)區(qū)地改值證明，該方法提高了算法速度和精度。

研究結(jié)果表明，基于DEM改進(jìn)的雙線性插值法遠(yuǎn)區(qū)地改布格重力異常和高程相關(guān)度更高，對(duì)線性構(gòu)造和地層、巖體的邊界的識(shí)別精度高。從本區(qū)的布格重力異常圖對(duì)比中發(fā)現(xiàn)：隨著基于DEM改進(jìn)的雙線性插值法地改不同網(wǎng)格距的選擇和高差變大，異常局部會(huì)產(chǎn)生不同的形態(tài)變化，從而提高地質(zhì)解釋的精度，這是值得注意的一點(diǎn)。

圖8 基于RGIS自帶高程庫(kù)布格重力異常和高程相關(guān)系數(shù)分布Fig.8 Based on the distribution of elevation kupger gravity anomaly and elevation correlation coefficient in RGIS

圖9 基于DEM改進(jìn)的雙線性插值布格重力異常和高程相關(guān)系數(shù)分布Fig.9 Improved bilinear interpolation kupger gravity anomaly and elevation correlation coefficient distribution based on DEM