馮進(jìn)凱，王慶賓，黃佳喜，張 超

(信息工程大學(xué)，河南 鄭州 450000)

引用著錄：馮進(jìn)凱，王慶賓，黃佳喜,等.重力梯度之地形改正方法與研究[J].測繪工程，2017，26(12)：50-54.

DOI:10.19349/j.cnki.issn1006-7949.2017.12.010

重力梯度之地形改正方法與研究

馮進(jìn)凱，王慶賓，黃佳喜，張 超

(信息工程大學(xué)，河南 鄭州 450000)

簡要介紹地形起伏對地面點(diǎn)重力梯度的影響，推導(dǎo)棱柱法計(jì)算重力梯度改正的嚴(yán)格積分公式并和積分核累加法對比驗(yàn)證了其正確性；將積分區(qū)域劃分為中央?yún)^(qū)和遠(yuǎn)區(qū)，中央?yún)^(qū)采用嚴(yán)格的棱柱法進(jìn)行計(jì)算，遠(yuǎn)區(qū)采用線性積分的方式逼近真實(shí)值，簡化了積分過程；根據(jù)重力梯度應(yīng)用的要求針對不同類型地形設(shè)定了中央?yún)^(qū)劃分準(zhǔn)則。實(shí)驗(yàn)表明：根據(jù)文中提出的方法和分區(qū)準(zhǔn)則，相比于傳統(tǒng)的積分法和棱柱法，計(jì)算速度得到提高，且精度要求相比于積分核累加法優(yōu)于0.01 E。

重力梯度；地形改正；分區(qū)準(zhǔn)則

重力梯度是重力沿某一方向的變化率，它的空間分布和變化與地下構(gòu)造、分布密度有著很強(qiáng)的聯(lián)系[1]，廣泛用于礦產(chǎn)探測、地下軍事目標(biāo)識別和重力梯度匹配導(dǎo)航等領(lǐng)域[2]。特別在航空重力梯度測量中，需要高精度的地面控制點(diǎn)進(jìn)行控制和聯(lián)合平差。高精度高分辨率的地形數(shù)據(jù)為充分利用地形信息進(jìn)行重力梯度改正提供了可能。在完成大范圍的梯度地形改正過程中，過去通常采用直接積分法[3]，即將地形格網(wǎng)進(jìn)行分割，然后對積分核累加求和，這種方法耗時(shí)嚴(yán)重，難以在短時(shí)間內(nèi)完成。例如：使用30″×30″的地形數(shù)據(jù)，積分范圍取2°×2°，在CD-4360小型計(jì)算機(jī)上完成百萬圖幅格網(wǎng)點(diǎn)梯度地形改正所需時(shí)間超過400 h[1]。近些年來，國內(nèi)外學(xué)者將快速Fourier變換(FFT)應(yīng)用其中。Tziavos利用FFT算法進(jìn)行了梯度地形改正的計(jì)算[3]；Kass也推進(jìn)了這一領(lǐng)域的研究[4]，使得計(jì)算速度成百倍的提高，但是精度略有損失，特別是在計(jì)算機(jī)性能高速發(fā)展的今天，其優(yōu)勢不太明顯；邊少鋒等將圓柱體模型和棱柱模型應(yīng)用于梯度改正，確定了梯度垂直分量Tzz的值[5]； 武凜等也進(jìn)行了相應(yīng)的研究[6]；LiZhi Zhu等將重力梯度正演模型計(jì)算方法進(jìn)行了歸類總結(jié)和比較[7]。

本文通過分析不同DEM數(shù)據(jù)特點(diǎn)的基礎(chǔ)上，將積分區(qū)域劃分為中央?yún)^(qū)和遠(yuǎn)區(qū)，并分別利用不同的方法進(jìn)行計(jì)算。中央?yún)^(qū)域采用棱柱積分的方法；遠(yuǎn)區(qū)采用線性卷積分的方法來逼近真實(shí)值。通過分析實(shí)驗(yàn)結(jié)果并針對不同地形數(shù)據(jù)提出了中央?yún)^(qū)劃分準(zhǔn)則，并在此基礎(chǔ)上進(jìn)行梯度地形改正計(jì)算。實(shí)驗(yàn)結(jié)果表明：此種改進(jìn)方法可以使得計(jì)算速度大幅提高，且精度與直接積分法相當(dāng)。

1 地面點(diǎn)重力梯度局部地形改正原理

1.1 傳統(tǒng)重力梯度地形改正的算法

梯度地形改正就是地面點(diǎn)(計(jì)算點(diǎn))周圍地形起伏對該點(diǎn)梯度的影響[6]，以計(jì)算點(diǎn)為原點(diǎn)，建立局部笛卡爾坐標(biāo)系，x軸正方向指向北，y軸正方向指向東，為表示方便x1=x，x2=y，x3=z。圖1所示在計(jì)算點(diǎn)P(ip,jp)處，重力梯度的地形改正可以表示為

z.

(1)

在計(jì)算中，通常將其化為累加形式，積分核累加形式如下：

(2)

其中:N，M分別為積分區(qū)域中積分單元的行數(shù)、列數(shù)，Q為高度方向劃分的單元格數(shù)。

圖1 局部地形重力梯度改正

1.2 改進(jìn)1：棱柱法求解梯度分量

根據(jù)式(1)可得梯度三方向的分量如下[6]：

(3)

現(xiàn)將式(3)中的函數(shù)直接積分求得原函數(shù)，以格網(wǎng)的邊界為上下限，每個(gè)地形格網(wǎng)只計(jì)算一次，從而避免大量的累加運(yùn)算。

地形網(wǎng)格數(shù)據(jù)步長近似為Δx=RΔB，Δy=RcosBmΔL，Bm為積分區(qū)域的平均緯度。ΔL和ΔB分別為球面網(wǎng)格的分辨率。選擇積分區(qū)域?yàn)镹×M，N為行數(shù)，M為列數(shù),為避免積分中的強(qiáng)奇異性問題，在此將計(jì)算點(diǎn)選在格網(wǎng)中心。

對式(3)積分后去離散求和形式如下：

(4)

其中：(i,j)為計(jì)算點(diǎn)編號，(ip,jp)為積分點(diǎn)編號，至此將復(fù)雜的三重積分化為了簡單形式，即在一定分辨率的地形數(shù)據(jù)下，每一個(gè)地形格網(wǎng)數(shù)據(jù)只需計(jì)算一次，避免了式(2)中復(fù)雜的累加運(yùn)算，將原來的N×M×Q次累加運(yùn)算降低為N×M運(yùn)算。

1.3 改進(jìn)2：劃分近遠(yuǎn)區(qū)求解梯度分量

在上述計(jì)算中，所有的距離量都是相對于積分點(diǎn)而言的，在距離計(jì)算點(diǎn)較遠(yuǎn)區(qū)域，h(i,j)-h(ip,jp)相對于距離量來說為小值。為加快計(jì)算速度，將

積分區(qū)σ劃分為中央?yún)^(qū)σ0和遠(yuǎn)區(qū)σ1=σ-σ0， 最后結(jié)果為中央?yún)^(qū)與遠(yuǎn)區(qū)地形改正之和：

(5)

(6)

2 數(shù)值計(jì)算分析

為研究不同地形數(shù)據(jù)對中央?yún)^(qū)劃分準(zhǔn)則的影響以及劃分準(zhǔn)則對計(jì)算速度的影響，本文搜集了3塊區(qū)域(1°×1°)，3塊區(qū)域的高程統(tǒng)計(jì)特性見表1。

本次采集數(shù)據(jù)為DEM數(shù)據(jù)，分辨率為3″×3″，選取了3塊地形特征明顯的區(qū)域。平原地區(qū)地形較為平緩，高程較低；山地地區(qū)地形落差大，高程普遍較高，地形起伏明顯。

表1 3塊區(qū)域高程統(tǒng)計(jì)信息

2.1 方法正確性檢驗(yàn)

積分核累加法是嚴(yán)格的數(shù)值計(jì)算方法，其正確性已經(jīng)得到廣泛檢驗(yàn)?？紤]到積分核累加法對時(shí)間的要求，本實(shí)驗(yàn)采取固定積分區(qū)域。本文提出算法直接使用中央?yún)^(qū)積分法。在積分核累加法中，將地形網(wǎng)格劃分為底面積1 m×1 m，高度為1 m的正方體，計(jì)算點(diǎn)區(qū)域?yàn)?′×3′，使用兩種方法計(jì)算得到差值見表2。

表2 兩種方法計(jì)算結(jié)果對比

根據(jù)式(3)嚴(yán)格的數(shù)學(xué)推導(dǎo)和上述實(shí)驗(yàn)兩種方法的差值比較，本文方法的正確性得到檢驗(yàn)。

2.2 中央?yún)^(qū)劃分準(zhǔn)則

通過數(shù)值初步估算，距離計(jì)算點(diǎn)20 km以外積分點(diǎn)對計(jì)算點(diǎn)的梯度影響幾乎為0，因此固定積分半徑為20 km，密度為2.67 g/cm3，中央?yún)^(qū)域采用棱柱法積分，遠(yuǎn)區(qū)采用簡化積分，通過改變中央積分區(qū)域大小，得到在3種地形數(shù)據(jù)下，重力梯度3個(gè)方向的值隨中央?yún)^(qū)半徑的變化。

圖2 區(qū)域1X 方向分量隨中央?yún)^(qū)的變化

圖3 區(qū)域1Y 方向分量隨中央?yún)^(qū)的變化

圖4 區(qū)域1Z 方向分量隨中央?yún)^(qū)的變化

圖5 區(qū)域2X 方向分量隨中央?yún)^(qū)的變化

圖6 區(qū)域2Y 方向分量隨中央?yún)^(qū)的變化

圖7 區(qū)域2Z 方向分量隨中央?yún)^(qū)的變化

圖8 區(qū)域3X 方向分量隨中央?yún)^(qū)的變化

圖9 區(qū)域3Y 方向分量隨中央?yún)^(qū)的變化

圖10 區(qū)域3Z 方向分量隨中央?yún)^(qū)的變化

圖2—圖4、圖5—圖7、圖8—圖10，分別對應(yīng)區(qū)域1(丘陵)、區(qū)域2(平原)、區(qū)域3(山地)的重力梯度三方向的改正值。實(shí)驗(yàn)數(shù)據(jù)表明：當(dāng)中央?yún)^(qū)積分半徑過小時(shí)，不能將高差看做距離的小值，所以計(jì)算結(jié)果存在一定的誤差。為突出最后收斂結(jié)果特征，將中央?yún)^(qū)半徑從1 km開始計(jì)算。

從圖中分析可得：隨著中央?yún)^(qū)的增大，重力梯度3個(gè)方向分量值呈收斂狀態(tài)，趨近真實(shí)值；當(dāng)積分半徑較小時(shí)，可以清楚看出，地形對梯度值的影響較為劇烈，而且不同地形的收斂節(jié)點(diǎn)不同，因此設(shè)定不同地形的不同中央?yún)^(qū)準(zhǔn)則是合理的；每一組三方向分量的計(jì)算值滿足拉普拉斯方程，即Txx+Tyy+Tzz=0，從側(cè)面證明了計(jì)算結(jié)果的正確性；三區(qū)域梯度值對比可以得出平原地區(qū)梯度值的收斂速度最快，丘陵地區(qū)次之，山區(qū)收斂速度最慢；要滿足梯度控制點(diǎn)和礦產(chǎn)探測的精度要求(±0.01 E)，現(xiàn)提出平原地區(qū)中央?yún)^(qū)積分半徑不小于5 km，丘陵地區(qū)中央?yún)^(qū)積分半徑不小于7 km，山地區(qū)域中央?yún)^(qū)積分半徑不小于9 km的中央?yún)^(qū)劃分準(zhǔn)則。這樣的劃分準(zhǔn)則，既滿足了點(diǎn)位的精度要求，又在一定程度上加快了計(jì)算速度；通過對比以上數(shù)據(jù)，不難發(fā)現(xiàn)，與地形高低起伏最為相關(guān)的是重力的垂直方向的絕對值，垂直梯度反映了重力在徑向的變化速率，因此有正有負(fù)。

2.3 基于以上算法和準(zhǔn)則的大范圍計(jì)算點(diǎn)實(shí)驗(yàn)

為保證積分法能在有限時(shí)間內(nèi)計(jì)算完成，所以本次試驗(yàn)區(qū)地形數(shù)據(jù)分辨率為1′×1′，選取3塊不同地形區(qū)域(2°×3°)，現(xiàn)通過積分核累加法和本文提出算法進(jìn)行計(jì)算，二者總積分區(qū)域半徑選為20 km，積分核累加法中將累加單元設(shè)為底面積邊長為300 m，高為1 m的正方體，密度為2.67 g/cm3。本文提出算法中，中央?yún)^(qū)化半徑分別為三塊區(qū)域設(shè)定的劃分準(zhǔn)則(平原5 km，丘陵7 km，山地9 km)，計(jì)算所用時(shí)間對比如表3所示。

表3 不同計(jì)算方法所用時(shí)間對比

積分核累加法雖然形式簡單且具有嚴(yán)格的數(shù)學(xué)意義，但是小單元累加中，長方體劃分越小，累加次數(shù)越多，計(jì)算速度越慢，山區(qū)尤為明顯。而本文介紹的棱柱法，每個(gè)網(wǎng)格數(shù)據(jù)中只需計(jì)算一次，采用本文介紹的方法，使得計(jì)算時(shí)間大大縮短，而且基于上述提出的分區(qū)準(zhǔn)則，計(jì)算時(shí)間還會進(jìn)一步縮減。

3 結(jié) 論

本文主要針對地形起伏對地表點(diǎn)的重力梯度的影響，采用不同方法進(jìn)行試驗(yàn)，得到如下結(jié)論：

1)在以計(jì)算點(diǎn)為原點(diǎn)的局部坐標(biāo)系中推導(dǎo)了梯度地形改正積分公式，大大簡化了計(jì)算流程，且試驗(yàn)驗(yàn)證了該方法的正確性，兩種方法差值均值絕對值不超過±0.01 E，且計(jì)算速度呈數(shù)十倍的提高。

2)推導(dǎo)了遠(yuǎn)區(qū)的積分公式，即在遠(yuǎn)區(qū)用線性卷積分的方法去逼近近似值，實(shí)驗(yàn)表明，采用該種分區(qū)方法計(jì)算時(shí)間會進(jìn)一步縮短。

3)提出基于不同地形條件下的中央?yún)^(qū)分區(qū)準(zhǔn)則，對以后采用更高分辨率地形數(shù)據(jù)細(xì)化分區(qū)準(zhǔn)則有一定的指導(dǎo)意義。

[1] YULE D E. Microgravity investigation of foundation conditions[J]. Geophysics, 1998, 63(1).

[2] 張赤軍.用地形數(shù)據(jù)確定重力異常垂直梯度[J].科學(xué)通報(bào),1999,44(6)：656-661.

[3] TZIAVO I N. The effect of the terrain on airborne gravity and gradiometry [J].Journal of Geophysical Research,1988,93 (B8):173-186.

[4] KASS M A, LI Y. Practical aspects of terrain correction in airborne gravity gradiometry surveys[J]. Exploration Geophysics，2008, 39: 198-203.

[5] 邊少鋒,張赤軍.地面擾動(dòng)重力垂直梯度的確定[J].地球物理學(xué)進(jìn)展,2005,20(4):969-973.

[6] 武凜,馬杰,周瑤,等.重力匹配導(dǎo)航的全張量重力梯度基準(zhǔn)圖模擬[J].系統(tǒng)仿真學(xué)報(bào),2009,21(22):7037-7041.

[7] CHRISTOPHER J, ZHU L Z. Comparison of methods to model the gravitational gradients from topographic data bases[J]. Geophys.J.Int,2006,166:999-1014.

[8] DAVID H．Evaluation of a full tensor gravity gradiometer for kimberlite exploration[A]. Airborne Gravity 2014一 Abstracts from the ASEG-PESA Airborne Gravity 2014 Workshop: Geoscience Australia Record[C]. 2014,18:73-79.

[9] 邊少鋒,張赤軍.地形起伏對重力垂直梯度影響的計(jì)算[J]. 物探化探計(jì)算技術(shù),1999, 21(2): 133-140.

[10] PAJOT G, DE V O,DIAMENT M, et al. LequentrecLalaneette4, and V. Mikhailov. Noise reduction through joint processing of gravity and gravity gradient data[J]. Geophysics, 2008, 73(3):123-134.

[責(zé)任編輯：劉文霞]

Methodsofterraincorrectionongravitygradiometry

FENG Jinkai, WANG Qingbin, HUANG Jiaxi, ZHANG Chao

(Information Engineering University,Zhengzhou 45000,China)

This paper presents the influence of gravity gradient caused by topographic relief on ground points briefly. Strict integral formula of gradient correction is deduced in Performance Prism, whose correctness is verified compared with Integral Nuclear Accumulation Method. The integration zone is divided in two parts, one of which is the center zone, and the other is the outside zone. The integration process is simplified using linear integral approximation to the real values in outside zone and time is shorten using the rules of partition for different types of terrain. Experiments show that computing time is shorter and accuracy is better than the 0.01 E compared with Integral Nuclear Accumulation Method based on the method and the partition criterion proposed in this paper.

gravity gradient；terrain correction；partition

P631

A

1006-7949(2017)12-0050-05

2016-10-22

馮進(jìn)凱(1992-)，男，碩士研究生.