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

馮進凱，王慶賓，黃佳喜，張 超

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

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

DOI:10.19349/j.cnki.issn1006-7949.2017.12.010

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

馮進凱，王慶賓，黃佳喜，張 超

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

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

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

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

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

1 地面點重力梯度局部地形改正原理

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

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

z.

(1)

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

(2)

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

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

1.2 改進1：棱柱法求解梯度分量

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

(3)

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

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

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

(4)

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

1.3 改進2：劃分近遠區(qū)求解梯度分量

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

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

(5)

(6)

2 數(shù)值計算分析

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

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

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

2.1 方法正確性檢驗

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

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

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

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

通過數(shù)值初步估算，距離計算點20 km以外積分點對計算點的梯度影響幾乎為0，因此固定積分半徑為20 km，密度為2.67 g/cm3，中央?yún)^(qū)域采用棱柱法積分，遠區(qū)采用簡化積分，通過改變中央積分區(qū)域大小，得到在3種地形數(shù)據(jù)下，重力梯度3個方向的值隨中央?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ù)據(jù)表明：當(dāng)中央?yún)^(qū)積分半徑過小時，不能將高差看做距離的小值，所以計算結(jié)果存在一定的誤差。為突出最后收斂結(jié)果特征，將中央?yún)^(qū)半徑從1 km開始計算。

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

2.3 基于以上算法和準(zhǔn)則的大范圍計算點實驗

為保證積分法能在有限時間內(nèi)計算完成，所以本次試驗區(qū)地形數(shù)據(jù)分辨率為1′×1′，選取3塊不同地形區(qū)域(2°×3°)，現(xiàn)通過積分核累加法和本文提出算法進行計算，二者總積分區(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)，計算所用時間對比如表3所示。

表3 不同計算方法所用時間對比

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

3 結(jié) 論

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

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

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

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é)通報,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] 邊少鋒,張赤軍.地面擾動重力垂直梯度的確定[J].地球物理學(xué)進展,2005,20(4):969-973.

[6] 武凜,馬杰,周瑤,等.重力匹配導(dǎo)航的全張量重力梯度基準(zhǔn)圖模擬[J].系統(tǒng)仿真學(xué)報,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]. 物探化探計算技術(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

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