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

（1. 中國地質科學院地球物理地球化學勘查研究所，河北廊坊 065000；2. 國家現代地質勘查工程技術研究中心，河北廊坊 065000）

0 引言

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

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

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

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

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

1 數據說明

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

表1 機載激光雷達DEM 統計表

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

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

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

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

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

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

扇形柱模型的地形改正為

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

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

其中

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

2.2 圓域數據與方域數據接口處理

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

區(qū)域重力調查規(guī)范給出的內接口補角地形改正的簡化公式為

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

其中

為權重系數。式中：xm、ym和Hm分別為補角范圍內代表高程特征的點的平面坐標和高程；xP、yP和HP分別為測點的平面坐標和高程值，通常計算Δhc需要四個高程數值和對應的權重系數，如圖2c 的A、B、C及D四點的高程數據及其權重系數。由式（7）可知，wA=

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

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

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

2.3.1 直棱柱模型解析式

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

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

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

2.3.2 質量線模型方法

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

式中：rkq為 質量線到測點的水平距離；(?x)kq、(?y)kq為某一個直棱柱的水平尺寸，由于(?x)kq、(?y)kq在累加求和符號內部，因此能夠對不同尺寸模型單元進行計算。

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

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

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

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

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

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

其中

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

2.4 大地坐標高程數據的處理

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

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

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

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

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

3 數值試驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 內接口補角計算方法對比

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

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

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

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

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

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

對試驗區(qū)統一采用質量線模型法進行中區(qū)地形改正，得到平均差異值統計數據（圖10）?？梢姡孩佼斁W格間距小于20 m×20 m時，無論中區(qū)地形改正范圍內邊界是20 m（圖10a）還是50 m（圖10b），平均差異值均小于20×10-8m/s2。②當內邊界從20 m（圖10a）增到50 m（圖10b），平均差異值整體減小。③內邊界為20 m、網格間距大于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；內邊界為50 m、網格間距大于40 m×40 m 時（圖10b），C 區(qū)的平均差異值減小（小于60×10-8m/s2)。④內邊界為50 m時（圖10b），網格間距為40 m×40 m 與50 m×50 m的計算結果高度相近，因此內邊界增大到一定程度可減弱網格間距對計算結果的影響。

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

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

4 結論

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

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

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

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

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

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