邱隆君，孫誠業(yè)*，楊亞斌，吳新剛

（1. 中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所，河北廊坊 065000；2. 國家現(xiàn)代地質(zhì)勘查工程技術(shù)研究中心，河北廊坊 065000）

0 引言

重力勘探是一種傳統(tǒng)的地球物理勘探方法，是對實測重力數(shù)據(jù)進(jìn)行自由空氣校正、中間層校正、地形校正、正常場校正等預(yù)處理獲得布格重力異常［1-2］，再用場分離技術(shù)和重力反演方法獲得地下密度異常體的位置信息或者密度信息［3-4］的一種間接勘探方法。

在數(shù)據(jù)預(yù)處理階段，為了消除起伏地形密度體的干擾，需要進(jìn)行地形改正。它的計算過程主要考慮兩個方面：第一，對起伏地形選擇恰當(dāng)?shù)拿芏葦?shù)值。既可選擇能代表研究區(qū)地層特征的密度值，也可以按照有關(guān)規(guī)范設(shè)定巖石密度、海水密度、淡水密度［5-6］。第二，結(jié)合實際的地形改正范圍選擇合適的地形數(shù)據(jù)庫。一般近區(qū)和中區(qū)地形改正范圍為0～2000 m，采用空間分辨率和精度較高的地形數(shù)據(jù)庫，如應(yīng)用高精度差分衛(wèi)星定位測地技術(shù)（GPS-RTK）或航攝獲得的高精度數(shù)字高程模型（DEM）數(shù)據(jù)［7］；遠(yuǎn)區(qū)地形改正范圍為2000～166700 m，需要覆蓋率較高的地形數(shù)據(jù)庫，如航天飛機(jī)雷達(dá)地形測繪使命（SRTM）或高級星載熱發(fā)射和反射輻射儀（ASTER）得到的DEM 數(shù)據(jù)。

地形數(shù)據(jù)庫的空間分辨率、精度、覆蓋率都會影響計算結(jié)果，分析某單一變量難以清晰地論述這些關(guān)系。因此，地形數(shù)據(jù)庫是全面、深入地研究地形改正的各種影響因素及其影響程度的關(guān)鍵要素。全站儀或GPS-RTK 方式獲取的高程數(shù)據(jù)精度最高，但成本高，工作效率低；航天遙感技術(shù)與空間測量分析技術(shù)的快速發(fā)展將地面測量獲得DEM 轉(zhuǎn)變?yōu)橛芍鲃舆b感方法自動生成DEM［8］，其中激光雷達(dá)（LiDAR）在獲取大范圍、高精度、多種DEM 方面的效率較高，明顯提升了大比例尺地形測繪和勘察效率［9］。

一般假設(shè)地形體的密度為常數(shù)，基于DEM 地形改正的主要研究內(nèi)容包括地形數(shù)據(jù)水平分辨率與計算精度之間的關(guān)系［10-12］、不同數(shù)據(jù)源的地形改正精度要求［7］。由于單一源的高分辨率數(shù)據(jù)的覆蓋范圍有限，地形改正方法的精度評價都是基于某一個數(shù)學(xué)模型的假設(shè)條件，只能討論高程測量誤差對地形改正的影響，不能分析由扇形錐、扇形柱等數(shù)學(xué)模型近似實際地形產(chǎn)生的誤差。目前，利用LiDAR 數(shù)據(jù)分析陸地重力近區(qū)地形改正數(shù)學(xué)模型產(chǎn)生的近似誤差以及中區(qū)地形改正的影響因素的文獻(xiàn)較少。

本文采用美國地調(diào)局三維高程計劃（3DEP）的LiDAR 高程數(shù)據(jù)，選擇12個數(shù)據(jù)區(qū)作為研究對象，利用LiDAR 高程數(shù)據(jù)為參考數(shù)據(jù)，以直棱柱模型重力解析式計算結(jié)果為地形改正參考值，從統(tǒng)計角度對比近區(qū)和中區(qū)地形改正的數(shù)學(xué)模型和計算方法的精度，這是本文的主要創(chuàng)新點；在此基礎(chǔ)上，分析不同計算方法、不同地形網(wǎng)格間距的差異，并且對比不同方法的計算效率。

1 數(shù)據(jù)說明

美國地質(zhì)調(diào)查局從2016年至2023年實施3DEP，目前可公開訪問的高程數(shù)據(jù)的最高精度為1 m×1 m，即柵格尺寸為1 m×1 m，數(shù)據(jù)維度為10012×10012，高程值為由機(jī)載LiDAR 掃描獲得的初始數(shù)據(jù)經(jīng)濾波得到的裸露地表高程［13］。本文選用3DEP 高程數(shù)據(jù)集中地形起伏特征明顯的12 個DEM 為試驗數(shù)據(jù)區(qū)，其地理位置分布見圖1。在生產(chǎn)過程中會融合其他數(shù)據(jù)，因此最終數(shù)據(jù)產(chǎn)品的精度并不統(tǒng)一。基于地貌劃分原則，試驗區(qū)的地貌包括丘陵（海拔<500 m，相對高度<200 m）、中山（海拔為800～3500 m，相對高度>500 m）、高山（海拔>3500 m，相對高度>500 m），詳細(xì)的統(tǒng)計數(shù)據(jù)見表1，數(shù)據(jù)下載網(wǎng)址為：https：∥www.usgs.gov/3d-elevation-program。

表1 機(jī)載激光雷達(dá)DEM 統(tǒng)計表

3DEP提供的原始數(shù)據(jù)的水平分辨率為1 m×1 m，為平面投影后數(shù)據(jù)，與ASTER GDEM（全球DEM）等以經(jīng)緯度表示的數(shù)據(jù)不同，更適合作為重力近區(qū)和中區(qū)地形改正的試驗參考數(shù)據(jù)。為分析不同網(wǎng)格間距對地形改正精度的影響，將1 m×1 m 柵格數(shù)據(jù)進(jìn)行“粗網(wǎng)格”處理，處理過程的本質(zhì)是“平滑”原始數(shù)據(jù)，使地形更“平緩”。以柵格數(shù)據(jù)存儲DEM 的數(shù)據(jù)空間排列都是有規(guī)律的，有些情況需要平面坐標(biāo)變換，如ASTER GDEM 高程數(shù)據(jù)，具體方法見下文。每個柵格像元的數(shù)值表示該像元覆蓋面積內(nèi)的平均高程，這與節(jié)點高程數(shù)據(jù)只表示單一位置高程信息的本質(zhì)不同。

2 近區(qū)與中區(qū)地形改正方法

2.1 近區(qū)地形改正方法

根據(jù)區(qū)域重力調(diào)查規(guī)范，1∶5 萬比例尺近區(qū)地形改正半徑為20 m，其他比例尺的近區(qū)地形改正半徑為50 m［14-15］。采用八方位圓域法計算時，內(nèi)環(huán)一般采用扇形錐模型（圖2a），外環(huán)采用扇形柱模型（圖2b）。扇形錐模型的地形改正為［16］

圖2 近區(qū)和中區(qū)地形改正采用的數(shù)學(xué)模型

式中：G、RD、nc、αc分別為萬有引力常數(shù)、近區(qū)地形改正半徑、方位數(shù)、扇形錐頂面中心線與其平面投影線間的夾角；σ為模型密度。

扇形柱模型的地形改正為

式中：Ri為第i環(huán)的半徑；Δh為扇形柱相對于測點的平均高差；np為方位數(shù)。

對于方域近區(qū)地形改正，重力規(guī)范中一般采用由測點P、角點B以及邊長中點C構(gòu)成的斜頂面三角棱柱模型（圖2e和圖2f）計算地形改正值，模型整體與扇形錐模型類似，不同的是該模型高度主要受角點和邊長中點的控制。基于該模型的方域近區(qū)地形改正為

其中

式中：h1，j、h2，j（j=1，2，…，8）分別為第j個斜頂面三角棱柱單元的中點和角點相對于測點的高差；Rs為方域近區(qū)地形改正半徑。

2.2 圓域數(shù)據(jù)與方域數(shù)據(jù)接口處理

分環(huán)、分方位進(jìn)行重力近區(qū)地形改正要將實際地形數(shù)據(jù)按照圓域數(shù)據(jù)處理，而高精度重力中區(qū)地形改正通常采用規(guī)則網(wǎng)格的高程數(shù)據(jù)。這些數(shù)據(jù)通常以附帶平面投影坐標(biāo)或者附帶大地坐標(biāo)信息的柵格數(shù)據(jù)存儲，在平面上坐標(biāo)數(shù)據(jù)通常以方域顯示。在調(diào)用這些方域數(shù)據(jù)對某個測點進(jìn)行中區(qū)地形改正時，還需要計算近區(qū)圓域數(shù)據(jù)與中區(qū)方域數(shù)據(jù)的間隙產(chǎn)生的地形改正值（圖2c的灰色充填部分），這一部分質(zhì)量產(chǎn)生的地形改正值要累加到近區(qū)圓域地形改正值中——數(shù)據(jù)接口處理，在重力調(diào)查規(guī)范中也稱為內(nèi)接口補(bǔ)角計算。

區(qū)域重力調(diào)查規(guī)范給出的內(nèi)接口補(bǔ)角地形改正的簡化公式為

式中Rm、Δhc分別為近區(qū)地形改正半徑、補(bǔ)角重心高程與測點高程差，而

其中

為權(quán)重系數(shù)。式中：xm、ym和Hm分別為補(bǔ)角范圍內(nèi)代表高程特征的點的平面坐標(biāo)和高程；xP、yP和HP分別為測點的平面坐標(biāo)和高程值，通常計算Δhc需要四個高程數(shù)值和對應(yīng)的權(quán)重系數(shù)，如圖2c 的A、B、C及D四點的高程數(shù)據(jù)及其權(quán)重系數(shù)。由式（7）可知，wA=

夏江海［17］給出補(bǔ)角地形改正的修正公式

2.3 中區(qū)地形改正方法

中區(qū)地形改正方法的實現(xiàn)取決于采用的數(shù)學(xué)模型。通常采用直棱柱模型解析式求解，或采用二重梯形積分快速計算，還可以采用傾斜頂面棱柱模型，利用表面積分法近似計算。理想模型數(shù)值模擬表明，高斯表面積分法的計算精度高于二重梯形積分公式?？蓪⒗庵w質(zhì)量壓縮至中心垂線上——質(zhì)量線模型，這種簡化模型的計算表達(dá)式簡潔、計算速度較快。以下討論直棱柱模型解析式的離散形式、質(zhì)量線模型方法、直棱柱模型的二重梯形積分方法及基于三角面元的表面積分方法。

2.3.1 直棱柱模型解析式

根據(jù)區(qū)域重力調(diào)查規(guī)范［14-15］，中區(qū)地形改正采用平頂面棱柱模型計算公式，某點的水平坐標(biāo)和高程為(xk，yq，hkq)，即表示x方向第k個、y方向第q個直棱柱頂面中心坐標(biāo)及其高程；測點的水平坐標(biāo)和高程值為(xP，yP，hP)。假設(shè)該測點進(jìn)行中區(qū)地形改正的范圍共包含J×I個矩形棱柱，那么全部直棱柱產(chǎn)生的總地形改正為

式中：(x，y，z)為相對坐標(biāo)，坐標(biāo)原點位于直棱柱頂面中心；r為直棱柱頂面中心到測點的相對距離；?h=hkq-hP，表示直棱柱頂面中心與測點的高程差；?x、?y為直棱柱的水平尺寸。

應(yīng)用上述公式計算測點的中區(qū)地形改正時，需要排除近區(qū)范圍的直棱柱模型的影響，避免重復(fù)計算。

2.3.2 質(zhì)量線模型方法

基于質(zhì)量線模型的計算是將三重積分進(jìn)一步簡化為離散的二重積分形式，即

式中：rkq為 質(zhì)量線到測點的水平距離；(?x)kq、(?y)kq為某一個直棱柱的水平尺寸，由于(?x)kq、(?y)kq在累加求和符號內(nèi)部，因此能夠?qū)Σ煌叽缒Ｐ蛦卧M(jìn)行計算。

2.3.3 直棱柱模型的二重梯形積分方法

利用高斯公式將地形質(zhì)量引力場的三重積分轉(zhuǎn)化為二重表面積分，在測點數(shù)千米范圍內(nèi)進(jìn)行地形改正可忽略地球表面彎曲影響。由于地形體的外邊界的側(cè)面與高程基準(zhǔn)面互相垂直，其法線與z軸垂直，因此在側(cè)表面的積分結(jié)果為零，其余為地形高程基準(zhǔn)面（大地水準(zhǔn)面）以及起伏地形表面的積分結(jié)果，即

式中：rkq為直棱柱頂面中心到測點的水平距離；Sd表示高程基準(zhǔn)面；St表示地形起伏表面。當(dāng)用直棱柱單元劃分地形時，式（11）近似為梯形積分表達(dá)式

式中Ckq為選用的梯形積分系數(shù)，其中內(nèi)節(jié)點系數(shù)為1，邊緣點系數(shù)為0.5，外角點系數(shù)為0.25，內(nèi)角點系數(shù)為0.75，系數(shù)具體分布見圖2d。當(dāng)近區(qū)和中區(qū)間的邊界通過原始柵格時，需要測點附近的四個高程節(jié)點（高程數(shù)值仍為測點的高程值）的地形改正值，再進(jìn)行雙線性插值，從而獲得測點的地形改正近似計算結(jié)果，此方法稱為共用點法。

2.3.4 基于三角面元的表面積分方法

基于散度定理，將式（11）的體積分轉(zhuǎn)換為起伏地形和地形高程基準(zhǔn)面包圍區(qū)域的表面積分，這兩部分通過三角面元離散，再利用高斯積分等數(shù)值積分方法估算每個三角面元的積分近似結(jié)果，最后累加求和。地形改正近似表示為［18］

其中

式中：(xα，yα，zα)、(xβ，yβ，zβ)以及(xη，yη，zη) 分別為第n個三角面元的三個頂點的空間坐標(biāo)；為第u個節(jié)點的空間坐標(biāo)；1，2，3；u=1，2，…，NQ)為節(jié)點坐標(biāo)，NQ 為使用的高斯節(jié)點數(shù)目；Wu為求積權(quán)重系數(shù)；L為規(guī)則數(shù)字化網(wǎng)格的網(wǎng)格間距，L2/2為三角面元在水平面上的投影面積；NE 為離散的三角面元總數(shù)。本文采用三點高斯積分公式［19］，節(jié)點的面積坐標(biāo)分別為(0，0.5，0.5)、(0.5，0，0.5)、(0.5，0.5，0)，求積權(quán)重系數(shù)分別為1/3、 1/3、 1/3。

2.4 大地坐標(biāo)高程數(shù)據(jù)的處理

目前，公開訪問的ASTER GDEM、SRTM 等高程數(shù)據(jù)都采用WGS-84大地基準(zhǔn)，由理論地球位場模型EGM96 推導(dǎo)的大地水準(zhǔn)面為高程基準(zhǔn)，計算柵格像元的沿經(jīng)線方向尺寸為?Mx，與沿緯線方向尺寸為?My，即

式中：φe、?φe以及?L分別為柵格像元的地心緯度、地心緯度跨度以及經(jīng)度跨度；re為柵格像元到地心的距離。實際計算中，全部柵格像元水平尺寸的轉(zhuǎn)換過程是一項繁瑣的計算過程，由于中區(qū)范圍一般在測點周圍2 km 以內(nèi)，因此可以利用測點位置地形模型單元的水平尺寸代替測點周圍全部柵格像元的水平尺寸，可極大地減小計算量。某一測點位置附近的柵格像元尺寸為

式中：?B為柵格像元緯度跨度；NEGM96為利用WGS-84 橢球參數(shù)和EGM96 地球重力位模型計算的大地水準(zhǔn)面高度；φe可由測點的大地緯度B換算；rE為WGS-84 橢球面的半徑，它是φe與第二偏心率e′的函數(shù)，即

式中e′2=(a2-b2)/b2，其中參數(shù)a、b分別為參考橢球的長半軸和短半軸。φe與B滿足

綜上所述，近區(qū)地形改正方法既可采用扇形錐與扇形柱的混合模型，也可僅采用扇形錐模型，而斜頂面三角棱柱模型可避免圓域與方域數(shù)據(jù)接口問題。中區(qū)地形改正方法中直棱柱模型解析式的精度最高；梯形積分與三角面元積分都是三維體積分的二維化簡形式，其中三角面元積分通過調(diào)整積分階數(shù)提高計算精度，是較靈活的計算方法；質(zhì)量線模型方法是將模型化簡的方法。

3 數(shù)值試驗

3.1 誤差分析與數(shù)值試驗環(huán)境

重力勘探近區(qū)和中區(qū)地形改正的計算精度及其改善方法一直以來都是關(guān)注重點。圖3 為重力中區(qū)地形改正的工作流程。從實際地形通過某種測量手段（水準(zhǔn)測量、機(jī)載激光雷達(dá)等）獲得數(shù)字高程數(shù)據(jù)，或者直接從測繪局申請數(shù)字高程數(shù)據(jù)，該過程存在測量誤差（包括點位誤差和高程誤差）和DEM 處理誤差。在獲得研究區(qū)完整DEM 后，利用某種數(shù)學(xué)模型進(jìn)行重力近區(qū)和中區(qū)地形改正，該過程會引入模型近似誤差及其計算誤差，這兩種誤差通常難以分離。本文默認(rèn)模型近似誤差包括計算誤差，作為主要研究對象。

圖3 重力中區(qū)地形改正的工作流程

區(qū)域重力規(guī)范［14-15］一般都假定密度為常數(shù)，基于某種數(shù)學(xué)模型及其公式討論高程精度以及點位誤差對重力地形改正精度的影響，如楊亞斌等［20］利用扇形、錐形和扇形柱公式計算理論誤差。

本文以機(jī)載LiDAR 高程數(shù)據(jù)為參考，衡量不同數(shù)學(xué)模型產(chǎn)生的模型近似誤差以及網(wǎng)格間距對地形改正精度的影響。假定地形密度為常數(shù)，忽略實際地層密度不均勻產(chǎn)生的影響，全部數(shù)值試驗均在Intel Core i7-11800H CPU（8 核、基準(zhǔn)頻率為 2.30 GHz）、64 GB RAM 的移動工作站、Windows 10 操作系統(tǒng)的計算環(huán)境下完成，采用基于MEX 文件的Matlab 2020b 并行計算方法。

3.2 近區(qū)地形改正的模型近似誤差

以往評價近區(qū)地形改正精度時以扇形錐、扇形柱以及三角棱柱等理論模型獲得的地形改正值為參考值［20］，忽略數(shù)學(xué)模型近似產(chǎn)生的差異。本文以1 m×1 m 間距的高程數(shù)據(jù)為參考數(shù)據(jù)、直棱柱模型重力解析式計算結(jié)果為參考值，評價常用數(shù)學(xué)模型的近區(qū)地形改正精度。采用八方位圓域或者方域三角棱柱計算方法，近區(qū)地形改正的范圍分為0～20 m與0～50 m兩種，形成7種計算方案（表2）。可見，既可以單獨使用扇形錐或者三角棱柱單一模型進(jìn)行近區(qū)地形改正（Ⅰ、Ⅱ、Ⅲ與Ⅳ），也可以內(nèi)環(huán)采用扇形錐、中環(huán)和外環(huán)采用扇形柱的混合模型方案（Ⅴ、Ⅵ與Ⅶ）。為模擬野外情況做如下約定：基于1 m×1 m網(wǎng)格數(shù)據(jù)，以間隔45°讀取1 個高程數(shù)據(jù)。如方案Ⅶ中采用3環(huán)，環(huán)的外邊緣到測點距離分別為10、25 和50 m，利用計算機(jī)程序讀取原始地形每個方位的半徑為10、17.5、37.5 m 圓上的高程節(jié)點值［21］，內(nèi)環(huán)利用扇形錐、中環(huán)和外環(huán)利用扇形柱模型計算后累加求和。

表2 近區(qū)地形改正的七種計算方案

近區(qū)地形改正值試驗設(shè)定如下：每個試驗區(qū)中央6000 m×6000 m 范圍內(nèi)，平均分布間距為100 m的3721(61×61)個測點，測點高程為1 m×1 m 網(wǎng)格的高程值，12 個試驗區(qū)一共有44652 個測點。由于近區(qū)范圍為圓域，在近區(qū)分界線的直棱柱模型只有一部分處于近區(qū)范圍，為保證計算嚴(yán)謹(jǐn)性，將這些模型單元細(xì)分為更小的棱柱，利用質(zhì)量線模型方法計算小棱柱體在測點的地形改正值。

計算每種計算方案的計算值與參考值的差異，利用該差異的最大值（10-8m/s2）和平均相對誤差進(jìn)行評價。差異最大值反映了每一個試驗區(qū)3721 個測點在近區(qū)地形改正時由模型近似產(chǎn)生的最大誤差；平均相對誤差定義為試驗區(qū)內(nèi)各測點差異值的絕對值的總和與各測點參考值的總和的比值（用百分比表示），是整體評價不同方案計算精度的指標(biāo)。

圖4 為近區(qū)地形改正不同計算方案的差異最大值對數(shù)曲線。由圖可見：①不同計算方案得到的差異最大值與近區(qū)地形改正半徑有直接關(guān)系，差異最大值通常集中在Ⅱ和Ⅳ，即使八方位Ⅶ的最大差異值在某些試驗區(qū)仍高于Ⅰ、Ⅲ和Ⅴ。②單一模型方案（Ⅰ、Ⅱ、Ⅲ和Ⅳ）的最大差異值較大，混合模型方案有效降低了最大差異值。③差異最大值與地貌密切相關(guān)，如Ⅱ、Ⅳ和Ⅵ中，大于550×10-8m/s2的差異最大值出現(xiàn)在高山（C 區(qū)和G 區(qū)）。值得注意的是，雖然K 區(qū)（丘陵）的高程數(shù)據(jù)標(biāo)準(zhǔn)差只有36.2 m，但是其差異最大值與J 區(qū)（中山）相當(dāng)。由于人工開采，L 區(qū)形成環(huán)狀階梯式臺階，高程范圍為［1472 m，3081 m］，標(biāo)準(zhǔn)差為358.4 m，是相對復(fù)雜的地貌區(qū)，然而不同計算方案的差異最大值卻處于一個量級（最小為135×10-8m/s2，最大為187×10-8m/s2)，即不同方案的差異最大值變化不大。以上分析說明，不能單獨以高程數(shù)據(jù)的范圍和標(biāo)準(zhǔn)差判定地貌對近區(qū)地形改正的影響，如地形起伏較平緩的丘陵的差異最大值與中山相當(dāng)，而地形復(fù)雜的露天銅礦區(qū)的差異最大值大致處于同一量級。

圖4 近區(qū)地形改正不同計算方案的差異最大值對數(shù)曲線

圖5 為近區(qū)地形改正不同計算方案的平均相對誤差曲線。由圖可見：①Ⅰ、Ⅲ、Ⅴ的平均相對誤差曲線的變化趨勢一致，計算精度從高到低的排序為Ⅴ、Ⅲ、Ⅰ，特別是Ⅴ可以保證K 區(qū)以外的11 個試驗區(qū)的平均相對誤差小于20%。②Ⅱ、Ⅳ、Ⅵ、Ⅶ的平均相對誤差曲線的變化趨勢相近，計算精度從高到低的排序為Ⅶ、Ⅵ、Ⅳ、Ⅱ，其中八方位Ⅶ使全區(qū)的平均相對誤差小于13%，計算精度最高。③在相同計算方案情況下，K 區(qū)的平均相對誤差高于其他區(qū)；L 區(qū)的平均相對誤差總體上小于15%。

圖5 近區(qū)地形改正不同計算方案的平均相對誤差曲線

除了差異最大值和平均相對誤差外，也可以使用絕對誤差值小于50×10-8m/s2的測點數(shù)與總測點數(shù)的比值（頻率）評價近區(qū)地形改正的精度。圖6為A 區(qū)Ⅳ、Ⅵ和Ⅶ的絕對誤差直方圖及絕對誤差—地形標(biāo)準(zhǔn)差散點圖。由圖可見：①相對于Ⅳ，Ⅵ和Ⅶ的絕對誤差更趨向于零誤差區(qū)間分布；Ⅳ、Ⅵ和Ⅶ的頻率分別為0.554、0.772 和0.926，因此八方位Ⅶ是非常有效的計算方案，將使約90%測點的模型絕對誤差小于50×10-8m/s2。②地形標(biāo)準(zhǔn)差小于20 m 的測點數(shù)占總測點數(shù)的98.1%，地形標(biāo)準(zhǔn)差小于10 m 的測點數(shù)占總測點數(shù)的45.2%；由于Ⅵ和Ⅶ引入扇形柱體模型，因此誤差散點逐漸向零誤差線收斂、聚集。文中將Ⅳ、Ⅵ和Ⅶ用于12 個試驗區(qū)并統(tǒng)計了頻率（圖7）。可見，除了C 區(qū)與G 區(qū)，在其他試驗區(qū)Ⅶ的頻率都大于0.97；與Ⅳ相比，Ⅶ明顯提升了12 個試驗區(qū)的頻率，即頻率提升最小值為0.059、最大值為0.539、平均值為0.274。

圖6 A 區(qū) Ⅳ（上）、Ⅵ（中）和Ⅶ（下）的絕對誤差直方圖（左）及絕對誤差—地形標(biāo)準(zhǔn)差散點圖（右）

圖7 近區(qū)地形改正不同計算方案的頻率

3.3 內(nèi)接口補(bǔ)角計算方法對比

根據(jù)式（5）與式（8），對12 個試驗區(qū)的全部測點的東北補(bǔ)角進(jìn)行計算，采用質(zhì)量線公式計算補(bǔ)角地形改正的參考值。對比每個試驗區(qū)3721 個測點的計算結(jié)果與參考值，并統(tǒng)計每個試驗區(qū)全部測點的差異最大值（10-8m/s2）以及平均相對誤差（圖8）。可見：式（8）降低了補(bǔ)角地形改正的差異最大值，幅值降低了約(2 ～ 8)×10-8m/s2(圖8a)；若不計入J 區(qū)與K 區(qū)的計算結(jié)果，由式（8）得到的平均相對誤差都小于20%（圖8b）。表1 表明，J 區(qū)與K 區(qū)為地形最平緩的地區(qū)，因此補(bǔ)角地形改正值很小，一般都在10-11m/s2量級。值得注意的是，K 區(qū)作為典型的丘陵區(qū)，地形相對起伏只有36.2 m，但由式（8）得到的平均相對誤差大于90%，因此推斷式（8）不適用于丘陵地區(qū)。相比而言，式（5）簡潔、計算穩(wěn)定性更好，適用性更強(qiáng)。

圖8 內(nèi)接口補(bǔ)角地形改正的差異最大值（a）及平均相對誤差（b）

3.4 高程數(shù)據(jù)網(wǎng)格間距和計算方法對中區(qū)地形改正的影響

重力勘探中區(qū)地形改正主要分為圓域計算或方域計算，前者利用扇形柱體模型，多用于早期手工數(shù)據(jù)成圖階段或圓域電子量板地形改正，后者利用規(guī)則網(wǎng)格劃分的直棱柱模型，便于計算機(jī)編程計算。當(dāng)前，方域計算已成為主要方式。參考地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn)［15］，區(qū)域重力調(diào)查中區(qū)地形改正范圍一般為50～2000 m，而1∶5 萬比例尺的重力中區(qū)地形改正范圍則為20 ～2000 m。

對12 個試驗區(qū)的原始1 m×1 m 高程數(shù)據(jù)進(jìn)行粗網(wǎng)格化操作，分別得到5 m×5 m、10 m×10 m、20 m×20 m、40 m×40 m 以及50 m×50 m 共5 種網(wǎng)格間距的高程數(shù)據(jù)，對每個試驗區(qū)分別采用直棱柱模型解析式法、質(zhì)量線模型、梯形積分法以及三角面元高斯積分法進(jìn)行計算，進(jìn)而分析網(wǎng)格間距和計算方法對中區(qū)地形改正的影響。每個試驗區(qū)平均分布間距為400 m 的256(16×16) 個測點，測點高程采用試驗區(qū)相同水平位置的高程值。將粗網(wǎng)格處理得到的低分辨率DEM 的地形改正值與測點參考值的差異作為研究對象，即利用每個試驗區(qū)的256 個測點的平均差異值評價四種計算方法的誤差。為了簡單起見，給出C 區(qū)、I 區(qū)、K 區(qū)和L 區(qū)的計算結(jié)果（圖9）?？梢姡孩俨煌W(wǎng)格間距計算結(jié)果差異較大。當(dāng)網(wǎng)格間距為1 m×1 m 和5 m×5 m 時，四種計算方法的平均差異值均小于20×10-8m/s2，三角面元高斯積分法和質(zhì)量線模型法均優(yōu)于梯形積分法；當(dāng)網(wǎng)格間距大于或等于20 m×20 m時，四種計算方法的平均差異值明顯增大，不同計算方法的增幅不同。②K 區(qū)計算結(jié)果（圖9c）的統(tǒng)計特征明顯不同于其他區(qū)。C 區(qū)（圖9a）、I 區(qū)（圖9b）和L 區(qū)（圖9d）的平均差異曲線趨勢類似；當(dāng)內(nèi)邊界為20（圖9c 左）、50 m（圖9c 右）時，K 區(qū)三角面元高斯積分法的平均差異值分別大于280×10-8、360×10-8m/s2，其他方法的平均差異值均小于15×10-8m/s2，說明三角面元高斯積分法不適用于丘陵區(qū)。③網(wǎng)格間距對不同計算方法的影響不同。隨著網(wǎng)格間距增大，總體上平均差異值也增大，其中梯形積分法的平均差異值曲線斜率最小，說明該方法較穩(wěn)定；當(dāng)網(wǎng)格間距為1 m×1 m～20 m×20 m 時，三角面元高斯積分法的曲線斜率大于其他計算方法，說明該方法受網(wǎng)格間距的影響較大。④內(nèi)邊界從20 m（圖9a 左～圖9d 左）增大到50 m（圖9a 右～圖9d 右），平均差異值整體降低。當(dāng)網(wǎng)格間距小于20 m×20 m時，圖9a 右～圖9d 右的平均差異值總體上與圖9a左～圖9d 左相同；當(dāng)網(wǎng)格間距小于40 m×40 m 時，三角面元高斯積分法的計算精度與質(zhì)量線模型法接近，且明顯高于梯形積分法，說明三角面元高斯積分法受中區(qū)地形改正內(nèi)邊界的影響較大；網(wǎng)格間距由40 m×40 m 增大到50 m×50 m 時，平均差異值基本不變，說明內(nèi)邊界增至50 m 可降低網(wǎng)格間距對計算結(jié)果的影響。

直棱柱模型解析式法、梯形積分法、質(zhì)量線模型法以及三角面元表面積分法的計算時間分別為457.93、48.69、57.87 和157.01 s，因此梯形積分法的計算效率最高，其中三角面元高斯積分法的計算時間和計算精度受高斯節(jié)點數(shù)量的影響（本文采用三個節(jié)點）。綜合考慮計算效率和計算精度，認(rèn)為質(zhì)量線模型法可作為中區(qū)地形改正的首選方法，特別是當(dāng)網(wǎng)格間距為5 m×5 m～20 m×20 m 時，該方法的平均差異值與直棱柱模型解析式法相近。

對試驗區(qū)統(tǒng)一采用質(zhì)量線模型法進(jìn)行中區(qū)地形改正，得到平均差異值統(tǒng)計數(shù)據(jù)（圖10）?？梢姡孩佼?dāng)網(wǎng)格間距小于20 m×20 m時，無論中區(qū)地形改正范圍內(nèi)邊界是20 m（圖10a）還是50 m（圖10b），平均差異值均小于20×10-8m/s2。②當(dāng)內(nèi)邊界從20 m（圖10a）增到50 m（圖10b），平均差異值整體減小。③內(nèi)邊界為20 m、網(wǎng)格間距大于40 m×40 m 時（圖10a），僅有D 區(qū)、E 區(qū)、H 區(qū)、G 區(qū)和K 區(qū)的平均差異值小于40×10-8m/s2，其他試驗區(qū)的平均差異值較大，其中C 區(qū)的平均差異值大于100×10-8m/s2；內(nèi)邊界為50 m、網(wǎng)格間距大于40 m×40 m 時（圖10b），C 區(qū)的平均差異值減?。ㄐ∮?0×10-8m/s2)。④內(nèi)邊界為50 m時（圖10b），網(wǎng)格間距為40 m×40 m 與50 m×50 m的計算結(jié)果高度相近，因此內(nèi)邊界增大到一定程度可減弱網(wǎng)格間距對計算結(jié)果的影響。

圖10 由質(zhì)量線模型法得到的地形改正范圍分別為20～2000 m（a）、50～2000 m（b）的中區(qū)地形改正平均差異值

綜上所述：為了保證高山區(qū)中區(qū)地形改正精度，應(yīng)采用較高水平分辨率的DEM，并設(shè)定地形改正范圍的內(nèi)邊界為50 m；采用網(wǎng)格間距為20 m×20 m 以及質(zhì)量線模型方法可以獲得較高的地形改正精度，并確保平均差異值小于20×10-8m/s2。

4 結(jié)論

本文基于包含美國典型地貌（喀斯喀特山脈、落基山脈、阿巴拉契亞山脈、德克薩斯丘陵、賓漢峽谷銅礦區(qū)）的1 m×1 m 激光雷達(dá)DEM，以直棱柱模型解析式計算結(jié)果為地形改正參考值，對比、分析了重力勘探近區(qū)和中區(qū)地形改正的模型近似誤差和計算方法誤差；基于5 種低分辨率網(wǎng)格高程數(shù)據(jù)和4 種計算方法，分析了網(wǎng)格間距、中區(qū)地形改正內(nèi)邊界和計算方法對中區(qū)地形改正的影響。得到如下認(rèn)識。

（1）采用八方位圓域、方域等7種計算方案的近區(qū)地形改正統(tǒng)計結(jié)果表明，近區(qū)地形改正范圍與計算模型是影響計算結(jié)果的主要因素。使用混合模型方案可有效降低模型近似的計算誤差，當(dāng)近區(qū)地形改正半徑為20 m 時，采用混合模型的方案Ⅴ可以保證K 區(qū)以外的11 個區(qū)的平均相對誤差小于20%；當(dāng)近區(qū)地形改正半徑為50 m 時，采用八方位方案Ⅶ可以將全區(qū)的平均相對誤差控制在13%以內(nèi)，是計算精度最高的方案。

（2）近區(qū)地形改正計算結(jié)果表明，地形起伏較平緩的丘陵（K 區(qū)）的差異最大值與中山（J 區(qū)）相當(dāng)；在復(fù)雜地表的露天銅礦區(qū)（L區(qū)），7種計算方案的差異最大值大致處于相同量級。

（3）補(bǔ)角地形改正的修正公式（式（8））不適用于丘陵。相比而言，區(qū)域重力調(diào)查規(guī)范給出內(nèi)接口補(bǔ)角地形改正的簡化公式（式（5））的平均相對誤差約為30%，表達(dá)式簡潔、計算穩(wěn)定性更好，適用性更強(qiáng)。

（4）影響中區(qū)地形改正計算結(jié)果的因素包括高程數(shù)據(jù)網(wǎng)格間距、中區(qū)地形改正內(nèi)邊界和計算方法，其中網(wǎng)格間距是主要因素。當(dāng)網(wǎng)格間距小于20 m×20 m 時，質(zhì)量線模型法地形改正結(jié)果的平均差異值均小于20×10-8m/s2，中區(qū)地形改正內(nèi)邊界明顯影響三角面元高斯積分法的計算精度。

（5）當(dāng)中區(qū)地形改正范圍為50～2000 m 時，計算精度從高到低的排序為直棱柱模型解析式法、質(zhì)量線模型法、三角面元高斯積分法、梯形積分法，計算時間分別為457.93、57.87、157.01、48.69 s。綜合計算精度和計算時間兩個因素，質(zhì)量線模型法為首選計算方法。