顧及地形效應的重力向下延拓模型分析與檢驗

劉　敏，黃謨濤，歐陽永忠，鄧凱亮，翟國君，吳太旗

1. 信息工程大學地理空間信息學院，河南 鄭州 450001； 2. 海軍海洋測繪研究所，天津 300061； 3. 海軍工程大學導航工程系，湖北 武漢 430033

?

顧及地形效應的重力向下延拓模型分析與檢驗

劉敏1，黃謨濤2,3，歐陽永忠2，鄧凱亮2，翟國君2,3，吳太旗2

1. 信息工程大學地理空間信息學院，河南 鄭州 450001； 2. 海軍海洋測繪研究所，天津 300061； 3. 海軍工程大學導航工程系，湖北 武漢 430033

Foundationsupport：TheNationalBasicResearchProgramofChina(973Program) (No. 613219)；TheNationalNaturalScienceFoundationofChina(Nos. 41474012； 41174062； 41374018)；TheGreatScientificInstrumentDevelopmentProjectofChina(No.2011YQ12004503)

向下延拓是航空重力測量數(shù)據(jù)實際應用中必不可少的技術(shù)環(huán)節(jié)。向下延拓屬于不適定反問題，其解算過程具有較大的不確定性，故該問題一直是大地測量領域國內(nèi)外學者的研究熱點。本文深入分析研究了當前國內(nèi)外最具代表性的3種向下延拓計算模型的技術(shù)特點和適用條件，提出了應用超高階位模型、局部地形改正和移去—恢復技術(shù)顧及地形效應，以及位場延拓結(jié)果球面化曲面的工程化方法，重點探討了計算模型的穩(wěn)定性及數(shù)據(jù)觀測誤差對延拓計算結(jié)果的影響。通過理論分析、數(shù)值仿真和實測數(shù)據(jù)計算等手段，定量評估了不同向下延拓模型的解算精度及其可靠性。其主要結(jié)論是：傳統(tǒng)逆Poisson積分模型解嚴重受制于輸入數(shù)據(jù)觀測噪聲的干擾，在現(xiàn)有作業(yè)條件下，該模型至多只能用于1km以下高度的延拓解算；頻譜截斷積分和位模型加地改兩種延拓新模型具有良好的計算穩(wěn)定性，完全適用于2′分辨率和5km飛行高度條件下的航空重力測量數(shù)據(jù)向下延拓解算，其延拓計算精度可達2×10-5m/s2，可滿足各方面實際應用需求。

航空重力測量；向下延拓；地形效應；Poisson積分；有限帶寬頻譜；超高階位模型

精化大地水準面始終是現(xiàn)代物理大地測量學永恒的研究主題[1-2]。精密確定大地水準面需要聯(lián)合利用衛(wèi)星、航空、地面、海洋等多源重力和地形高數(shù)據(jù)，其計算過程涉及多源數(shù)據(jù)的融合處理問題。向下延拓是多源重力數(shù)據(jù)融合處理必不可少的技術(shù)環(huán)節(jié)，在數(shù)據(jù)準備階段，需要將衛(wèi)星和航空重力測量成果延拓到地面作聯(lián)合處理；在邊值問題解算階段，需要將地面重力數(shù)據(jù)延拓到大地水準面。由于物理場向下延拓在數(shù)學上屬于不適定反問題，求解此類問題存在很大的不確定性[3]，一直以來國內(nèi)外諸多學者為解決這一棘手問題付出了不懈的努力[4-28]。目前解決向下延拓問題主要有3種途徑，第1種是直接求逆泊松(Poisson)積分方程，是應用比較廣泛的主流途徑[29]，國內(nèi)外學者主要圍繞求逆過程引起的不穩(wěn)定性問題開展不同形式的正則化方法研究[10,13,16,18-23,25-26]；第2種是間接求逆途徑，包括最小二乘配置方法[30-32]、虛擬點質(zhì)量方法[1-2]和矩諧分析法[24]等，此類方法雖然避開了求逆Poisson方程，但仍涉及矩陣求逆過程，因此也不可避免存在不穩(wěn)定性問題[33-34]；第3種可統(tǒng)稱為非求逆途徑，文獻[17]提出的利用航空重力測量和DEM確定地面重力場的直接代表法，文獻[27—28]提出的聯(lián)合使用超高階位模型和地形信息確定向下延拓改正數(shù)的方法(以下簡稱差分延拓法)等都屬于非求逆途徑的范疇，此類方法不受傳統(tǒng)求逆過程不穩(wěn)定性的影響。文獻[29]基于帶限(band-limited)航空重力測量數(shù)據(jù)特有的頻譜特性，提出了向下延拓的直接積分公式(以下簡稱頻譜截斷法)，并將其推廣應用于大地水準面的直接解算[35]，其計算過程也避開了方程求逆問題。文獻[13]將頻譜截斷方法和基于求逆過程的各類正則化方法作了全面的數(shù)值比較和分析，得出的結(jié)論是，前者的計算精度和效率都明顯優(yōu)于其他方法。

如前所述，困擾向下延拓問題的關鍵是其本身固有的不適定性。文獻[9]曾對向下延拓的穩(wěn)定性問題做過深入分析，同時提出了改善計算穩(wěn)定性的地形效應補償方法；文獻[11，13，14，29，35]專門針對航空重力向下延拓問題，開展了大量有效的數(shù)值計算和比對工作，得出了一些極具參考價值的結(jié)論。由已有的研究成果得知，雖然向下延拓問題本身是不適定的，但其解算結(jié)果的不穩(wěn)定程度取決于延拓計算高度和數(shù)據(jù)分辨率(即網(wǎng)格間距大小)兩個方面。在一定條件下，向下延拓解算方程可以是良態(tài)的，而在超越一定界限以后，即使采用正則化處理技術(shù)，也無法取得有效的解算結(jié)果[11,13]。因此，對于實際應用來說，最要緊的是預先掌控前面所指的一定“條件”和“界限”，盡可能在規(guī)定的“條件”下開展航空重力測量作業(yè)，避免超越已知的“界限”，只有這樣才能獲得預期的測量成果。另一方面，地形效應對重力向下延拓解算結(jié)果具有重要影響，因為求逆過程的穩(wěn)定性不僅取決于積分方程系數(shù)矩陣的結(jié)構(gòu)，還取決于觀測向量的頻譜特性[9]，故可通過地形效應的“移去-恢復”運算來改變重力觀測量的頻譜特征，從而改善向下延拓解算的穩(wěn)定性。本文的目的是，通過理論分析、模擬仿真和實際數(shù)值計算等手段，對當前國內(nèi)外最具代表性的3種向下延拓模型進行全面的分析比較和適用性檢驗，力爭為工程應用提供具有可操作性的延拓計算方案。

1　計算模型及穩(wěn)定性分析

1.1逆Poisson積分迭代解延拓模型

由文獻[6]知，Poisson積分向上延拓公式為

(1)

(2)

向下延拓計算是求式(1)的逆問題，即已知Δgp，要求Δgd。對式(1)作離散化處理，可得一線性方程組，用矩陣形式表示為

y=Ax

(3)

式中，y為由航空重力測量獲得的已知重力異常向量；x為待求的球面重力異常向量；A為由式(1)積分核函數(shù)確定的系數(shù)矩陣?？刹捎萌缦碌腏acobi迭代法解算方程式(3)[9]，將A改寫為

A=E-B

(4)

代入式(3)有

x=y+Bx

(5)

取

x0=y

(6)

求

xi=Bxi-1(i>0)

(7)

(8)

式中，E為單位矩陣；ε為給定的某個大于零的限差值；k為迭代次數(shù)。

在實際應用中，由于觀測數(shù)據(jù)有限，一般將球面積分區(qū)域劃分為近區(qū)和遠區(qū)，近區(qū)是以計算點為中心、ψ0為半徑的球冠區(qū)域σ0，近區(qū)影響直接由觀測數(shù)據(jù)計算；遠區(qū)是球面上的剩余部分(σ-σ0)，遠區(qū)影響由位系數(shù)模型按下式計算[11,13]

(9)

(10)

(11)

式中，GM為地球引力常數(shù)；L為位系數(shù)模型的最高階數(shù)；l為參考場位模型的最高階數(shù)；a為參考橢球長半軸；Δgn代表重力異常n階面球諧函數(shù)；Qn為Poisson核截斷系數(shù)；Rn,m(ψ0)可由已知的遞推公式計算，具體參見文獻[36—37]。此時，對應于式(3)的離散形式為

(12)

式中，N為待求點個數(shù)，一般取與已知點個數(shù)相同。為了減小中心數(shù)據(jù)塊(即與計算點重合的數(shù)據(jù)塊)離散化誤差的影響，可按下式計算系數(shù)矩陣A的對角線元素[13]

(13)

矩陣A的非對角線元素為

(14)

式中，N0為位于積分球冠區(qū)σ0內(nèi)的數(shù)據(jù)點個數(shù)；Δσj為數(shù)據(jù)塊面積。

如前所述，式(1)只是近似的球面解模型，式(1)只能反解得到等高度面(即球面或近似為大地水準面)上的重力異常，而非地形面上的重力異常。如果是向下延拓到大地水準面上，那么由于大地水準面外地形質(zhì)量的存在，使得其延拓解與物理意義上的現(xiàn)實性不相對應，即這里的延拓解只是一組虛擬的重力異常[6]。這類重力異?？蓡为殤糜谖锢泶蟮販y量參數(shù)計算，但不宜與其他類型重力異常數(shù)據(jù)聯(lián)合使用。為了求得地形面上的重力異常，提出如下位場延拓球面化曲面方法：

(1) 將計算區(qū)域的地形高度變化范圍(Hmax-Hmin)，按整百米間隔(如100 m或200 m)劃分為若干個等高度面；

(2) 利用前面的迭代法依次計算出各個等高度面上的重力異常；

(3) 通過內(nèi)插方法計算得到位于兩個等高面之間的地形面上的重力異常。

此外，為了改善向下延拓解算的穩(wěn)定性，可通過地形效應的“移去-恢復”運算改變重力觀測量的頻譜特征[9]，即首先在航空重力測量成果中移去地形質(zhì)量對空間點的作用，然后在延拓計算結(jié)果中恢復地形質(zhì)量對地面點的影響。地形效應改正的計算模型可參見文獻[1，28]。

1.2頻譜截斷積分延拓模型

由于受到飛行環(huán)境動態(tài)效應的影響，航空重力測量原始觀測數(shù)據(jù)中一般包含上千甚至上萬毫伽(mGal=10-5m/s2)的干擾加速度信息[38-39]。為了消除高頻干擾噪聲的影響，通常需要對航空重力測量數(shù)據(jù)作低通濾波處理，以獲取所需的重力異常信息。但經(jīng)濾波處理后的重力測量成果必定會損失掉一部分有用的高頻重力信息，損失量值大小取決于濾波截斷頻率的選擇，實際應用中需要平衡好測量精度和分辨率兩方面的需求。由此可見，航空重力測量成果是經(jīng)過高頻截斷后的重力異常場信息，不包含某個頻率(L)以上的高頻分量。而在對航空重力測量成果作進一步處理(如向下延拓)時，通常還需要引入全球重力場模型(GGM)進行移去-恢復運算，即事先需要從航空重力測量成果中移去GGM的影響，然后在解算結(jié)果中作反向的恢復運算。這說明，實際參與向下延拓解算的航空重力測量數(shù)據(jù)是一類有限帶寬的重力場信息(band-limited airborne gravity data)[29]。設GGM的最高階次為(l-1)，則有L=l+b，b稱為航空重力測量數(shù)據(jù)的帶寬。如何合理利用航空重力測量數(shù)據(jù)有限帶寬的頻譜特性，是有效提高向下延拓解算過程穩(wěn)定性的關鍵。文獻[9、11、13—14、29、35]全面研究了等高飛行條件下的航空重力測量數(shù)據(jù)延拓問題，提出了如下相對穩(wěn)定的頻譜截斷積分向下延拓公式。

設在恒等高度H的飛行面上獲取了一組有限帶寬航空重力異常數(shù)據(jù)Δgb(R+H)，飛行面與大地水準面之間已經(jīng)不存在地形質(zhì)量，要求確定大地水準面(近似為半徑R的球面)上相對應的有限帶寬擾動位Tb(R)。理論上可將上述問題表述為如下的偽邊值問題[9,29]

從表2可以看出，參試品種的株高從72.5-88.0cm之間。對照(4個品種平均值)為83.2cm，青海13號最低，為72.5cm。4個品種均為直立生長型。4個品種的單株有效分枝數(shù)在1.3個-2.2個之間，對照為1.7個，以青海13號最高，為2.2個。4個品種的單株莢數(shù)在5.2個-7.5個之間，對照為6.0個，青海13號最多，為7.5個。4個品種的莢粒數(shù)在1.4個-2.2個之間，對照為1.8粒，青海13號最多，為2.2個。最少為青海12號，為1.4個。4個品種百粒重在79.8g-170.0g之間，對照為140.8g，青蠶14號最高，為170.0g。青海13號最低，為79.8g。

(15)

(16)

(17)

對式(16)中的雙邊值面作近似處理并轉(zhuǎn)換為單邊值面后[9]，可求得上述問題的擾動位解為[29]

Δgb(R+H,Ω′)dΩ′

(19)

Δgb(R+H,Ω′)dΩ′

(20)

(21)

(22)

(23)

式中，L0為位模型最高階次，且L0≤L，其他符號意義同前。對于近區(qū)計算，同樣需要對中心數(shù)據(jù)塊作精細化處理，以減小離散化誤差的影響。類似于式(13)和式(14)，此時對應于積分式(20)的系數(shù)矩陣Ab的對角線元素為

(24)

(25)

[Pn+1(cosψ0)-Pn-1(cosψ0)]/(2n+1)

(26)

矩陣Ab的非對角線元素為

(27)

式中其他符號意義同前。由式(20)可直接計算得到地形面上的重力異常，故無須作位場延拓平面化曲面處理，此時數(shù)據(jù)面高度保持不變(r=R+H)，計算面高度由R改變?yōu)?R+hi)，hi代表地面計算點的大地高。類似于前一小節(jié)的做法，這里同樣需要通過地形效應的“移去-恢復”運算，來消除計算面與數(shù)據(jù)面之間的地形質(zhì)量影響，以改善向下延拓解算的穩(wěn)定性。

1.3基于位模型和地形改正差分的延拓模型

為了規(guī)避傳統(tǒng)逆Poisson積分向下延拓解算過程的不適定性問題，文獻[27—28]借鑒導航定位中的“差分”概念，曾先后提出了利用超高階位模型和地形高信息，直接實施海域和陸部航空重力測量數(shù)據(jù)向下延拓計算的新方案。其中，陸部顧及地形效應延拓方案的核心思想是：以飛行高度面與地面對應點的位模型重力異常差分信息表征向下延拓總改正數(shù)的中長波分量，以相對應的局部地形改正差分修正量表征總改正數(shù)的中高頻成分，從而實現(xiàn)航空重力數(shù)據(jù)向地面點對點的全頻段延拓。具體延拓計算模型為[28]

(28)

(29)

(30)

1.4計算模型穩(wěn)定性分析

如前言所述，重力場向下延拓計算在數(shù)學上屬于不適定反問題，其解算過程存在不穩(wěn)定性是該問題本身固有的一種屬性[3]。理論上，由逆Poisson積分公式(1)組成的向下延拓解算方程屬于第一類弗雷德霍姆(Fredholm)積分方程，當且僅當滿足如下的Picard條件時，該方程存在唯一的收斂解[9]

(31)

式中，fi代表已知航空重力觀測量Δgp的傅里葉(Fourier)展開系數(shù)；λi代表Poisson核函數(shù)的特征值。式(31)說明，要想確保逆Poisson積分方程有解，從某個階次開始，F(xiàn)ourier系數(shù)fi的衰減速度必須快于特征值λi的衰減速度。

對于離散化形式的線性方程組即式(3)，通常采用如下的系數(shù)矩陣條件數(shù)指標來度量方程解的病態(tài)性[9,20]

(32)

式中，λmax和λmin分別代表系數(shù)矩陣A之特征值的最大值和最小值。一般認為[20]，當0<κ<100時，方程組是良態(tài)的；當100≤κ≤1000時，定義為中等程度病態(tài)；當κ>1000時，定義為嚴重病態(tài)。對于在恒等飛行高度H條件下的航空重力向下延拓問題，逆Poisson積分方程系數(shù)矩陣的條件數(shù)可近似表示為[9,11]

(33)

式中，π/ΔΩ稱為奈奎斯特(Nyquist)頻率，ΔΩ為數(shù)據(jù)網(wǎng)格間距。取ΔΩ=2′,5′，可分別計算得到對應于不同延拓高度H的條件數(shù)κ，具體結(jié)果如表1所示。必須指出的是，式(33)給出的條件數(shù)只是一種理論上的病態(tài)性度量指標，受各種擾動因素影響后，實際問題解的變化特性則要復雜得多，其穩(wěn)定性除了取決于系數(shù)矩陣的結(jié)構(gòu)外，還取決于已知觀測量的頻譜特性[9]。文獻[11]在完成大量的數(shù)值計算和分析比較后發(fā)現(xiàn)，航空重力數(shù)據(jù)向下延拓解算的穩(wěn)定性與數(shù)據(jù)網(wǎng)格間距和延拓高度的比值(α=ΔΩ/H)密切相關，在受到觀測噪聲干擾的情況下，只有當比值α>1.1時，直接利用航空重力數(shù)據(jù)確定的局部大地水準面解才是有效的。而航空重力異常直接向下延拓到地面結(jié)果的精度要比前者悲觀得多，即使采用一些正則化方法進行優(yōu)化處理，也無法徹底解決觀測噪聲的放大問題，要想獲得可接受的地面重力異常延拓結(jié)果，必須盡可能降低航空重力測量的飛行高度[11]。本文將在第2節(jié)繼續(xù)對此問題進行數(shù)值計算和檢驗。表1同時列出了數(shù)據(jù)網(wǎng)格間距分別取ΔΩ=2′,5′時，對應于不同延拓高度的比值α。

表1　不同延拓高度對應的條件數(shù)和網(wǎng)格間距與高度比值

2　數(shù)值計算與分析

為了檢驗上述各個計算模型的延拓效果，分別采用模擬仿真和實測數(shù)據(jù)兩種方式開展數(shù)值計算及分析比較研究。

2.1模擬仿真計算與分析

2.1.1試驗數(shù)據(jù)與計算方案

以美國本土一個3°×3°區(qū)塊作為試驗區(qū)，選用EGM2008位模型作為標準場模擬產(chǎn)生不同高度的2′×2′網(wǎng)格航空和大地水準面(即高度為零)重力異常。之所以選擇美國本土作為試驗區(qū)，是考慮到EGM2008位模型在美國地區(qū)具有更高的逼近度[40-41]，并與3.2節(jié)選用的實測數(shù)據(jù)試驗區(qū)取得一致。該區(qū)塊屬于地形變化比較劇烈的大山區(qū)，試驗效果具有一定的說服力。分別以1km、3km和5km3個高度的位模型重力異常作為觀測量，依次采用逆Poisson積分(簡稱模型1)和頻譜截斷積分(簡稱模型2)兩種計算模型，將3個不同高度面上的網(wǎng)格重力異常延拓到零高度面，并將其同由位模型計算得到的零高度面觀測量(作為基準值)進行比較，從而獲取相應的計算模型精度評估參數(shù)。由于基于位模型和地形改正差分的延拓模型(簡稱模型3)本身就是建立在超高階位模型之上的，因此該模型暫不參加本階段的數(shù)值計算檢驗。

為了進一步考察觀測噪聲對延拓計算結(jié)果的影響，特別設計在模擬觀測量中分別加入±1mGal和±3mGal的隨機噪聲，生成對應高度面上的帶噪聲的兩組觀測量，并重復前面的計算過程和對比分析。上述誤差量值與當前國內(nèi)外航空重力測量的精度水平相當[11,14,29,35,42]，因此其檢驗結(jié)論具有實用意義。本試驗統(tǒng)一采用EGM2008位模型的前360階次作為參考場(GGM)，對應于零高度面和3個不同高度面的EGM2008位模型(361～2160階次)殘差重力異常統(tǒng)計結(jié)果見表2。

表2　試驗區(qū)位模型殘差重力異常統(tǒng)計結(jié)果

2.1.2計算結(jié)果與分析

本試驗在一個高度面上共有網(wǎng)格觀測數(shù)據(jù)30×3×30×3=8100個，對應于模型1，需要求解8100階線性方程組，采用Jacobi迭代法解算時，迭代終止參數(shù)取ε=0.1mGal；對于模型2，模型參數(shù)分別取為l=361、L=L0=2160、b=1799；積分半徑統(tǒng)一取為ψ0=1°。按照前面設計的計算方案，對3個高度和兩個模型分別解算無誤差和有誤差干擾條件下的零高度面延拓重力異常，同時將其與零高度面基準值作比較，具體計算結(jié)果如表3所示。為了減小積分邊緣效應對統(tǒng)計結(jié)果的影響，計算區(qū)域邊緣1°范圍內(nèi)的數(shù)據(jù)不參加對比分析。

從表3結(jié)果可以看出，對于無誤差干擾條件下的輸入數(shù)據(jù)，模型1和模型2都能給出比較理想的輸出結(jié)果，基本不受計算高度大小的影響。但對于有誤差干擾條件下的輸入數(shù)據(jù)，模型1和模型2計算效果則表現(xiàn)出非常顯著的差異。在3km以下計算高度，模型2對數(shù)據(jù)誤差具有一定的抑制作用，只有當延拓高度增大到5km時，模型2才對數(shù)據(jù)誤差產(chǎn)生一定的放大效應，這一結(jié)果充分體現(xiàn)了模型2作為正解模型及其截斷核函數(shù)所具有的超強的抗干擾能力。模型1的解算結(jié)果則明顯受到數(shù)據(jù)誤差的干擾，在所有計算高度上，模型1對數(shù)據(jù)誤差都有放大作用，計算結(jié)果嚴重失真，當數(shù)據(jù)誤差大于1mGal，計算高度超過1km時，模型1的解算結(jié)果都是不可靠的。這一結(jié)果說明，雖然在理想情況下，通過迭代計算求解逆Poisson積分方程，也能收斂到問題的理論解[6,9]。但當觀測數(shù)據(jù)存在噪聲干擾時，由于反問題自身固有的不適定性，觀測噪聲將在迭代過程中得到累積和放大，使得較小的噪聲干擾也會引起問題解出現(xiàn)較大的變化，最終造成解算結(jié)果嚴重偏離真解。

表3　不同模型計算結(jié)果與基準值比對

2.2地面實測數(shù)據(jù)計算與分析

2.2.1試驗數(shù)據(jù)與計算方案

檢驗向下延拓計算模型適用性最有效的方法是，直接使用實測航空重力數(shù)據(jù)完成向下延拓計算，將計算結(jié)果與地面實測重力基準值進行比較，由此可得到不同計算模型的精度評價指標。遺憾的是，陸部同時擁有高精度高分辨率航空和地面實際重力觀測數(shù)據(jù)的區(qū)域并不多見，目前筆者還缺乏這方面的可靠資料。為此，本文改用向上與向下延拓比對方法對計算模型精度進行外部檢核，即首先利用地面網(wǎng)格重力和地形高數(shù)據(jù)，通過向上延拓方法計算得到一定高度面上的空中重力異常，將其作為航空重力測量的觀測量，進而使用不同的計算模型將其向下延拓到地面，求取延拓計算值與地面已知網(wǎng)格重力值的差異，便可獲得相應延拓計算模型的精度評價。

這里繼續(xù)選用與2.1節(jié)完全相同的美國本土3°×3°區(qū)塊作為試驗區(qū)，開展航空重力向下延拓的實際數(shù)值計算和對比分析。該區(qū)塊同時擁有2′×2′網(wǎng)格地面觀測重力異常和30″×30″網(wǎng)格地形高數(shù)據(jù)，兩組數(shù)據(jù)的變化特征如表4所示。

表4　試驗區(qū)塊實測重力和地形數(shù)據(jù)變化特征

使用“先向下后向上”方法[6,43]，首先將地面重力向下延拓到大地水準面，然后將其向上延拓到5km高度面。為了考察地形效應對向下延拓計算結(jié)果的影響，這里采用兩種途徑完成地面重力向下延拓解算，得到兩組大地水準面上的重力值，一種途徑是直接采用原始觀測重力異常Δg0(1)和式(3)進行純粹數(shù)學意義上的向下延拓計算，得到一組虛擬的重力異常Δg*(1)，此方法不擾亂外部重力場[6]；另一種途徑是首先從地面重力中扣除掉大地水準面以外地形質(zhì)量引力的影響，得到調(diào)整后的重力異常Δg0(2)，然后將其向下延拓到大地水準面，得到一組消除地形效應影響后的延拓重力異常Δg*(2)。最后利用Δg*(1)和Δg*(2)可由式(1)向上延拓得到對應5km高度面的兩組航空重力觀測量Δgp(1)和Δgp(2)，并根據(jù)模型1和模型2完成后續(xù)的向下延拓計算。但此時需要確定的是地形面上的重力異常，故向下延拓計算高度面是對應于地面網(wǎng)格點高度的等高面。對于模型1，需要利用前面提出的位場延拓球面化曲面方法，將球面上的計算結(jié)果內(nèi)插到地面上的各個測點；模型2可直接計算得到測點上的重力異常；模型3通過局部地形改正差分方式顧及了地形效應的影響，故只需參加由Δg0(1)生成的觀測量檢驗?？紤]到本文關注的重點是計算模型的適用性和地形效應的影響，不必過分追求計算結(jié)果的絕對精度，故可對向上和向下延拓模型參數(shù)設置作統(tǒng)一的近似處理。這里將大地水準面近似為球面，統(tǒng)一取ri=R+hi，rp=R+H，hi為地面點高程，其他符號意義同前。此時模型2的計算參數(shù)取為l=361，L=5400，b=5039,L0=2160，即與2′×2′數(shù)據(jù)分辨率相一致，其他模型參數(shù)保持不變。

2.2.2計算結(jié)果與分析

按照前面設計的計算方案，分別采用上述3個計算模型完成5km高度面各兩組觀測數(shù)據(jù)的向下延拓解算，將其分別與相對應的地面基準值Δg0(1)和Δg0(2)作比較，具體計算結(jié)果如表5所示。

表5　不同模型5 km高度向下延拓結(jié)果與基準值比對

表5結(jié)果顯示，模型1的檢核效果反而明顯好于模型2和模型3，這是由于作為觀測量的空中重力異常是由地面重力異常通過Poisson積分向上延拓得到的，它與模型1使用的逆Poisson積分形成閉環(huán)關系，其計算模型誤差有相互抵消作用，因此該比對結(jié)果不能作為評價模型一計算性能的依據(jù)，但相應結(jié)果至少說明本文提出的位場延拓球面化曲面方法和顧及地形效應策略是可行有效的。由于模型2和模型3的計算原理與Poisson積分向上延拓過程無關，故表5中對應于模型2和模型3的檢核結(jié)果是獨立有效的。從表5看出，當不顧及地形效應影響時，模型2的檢核精度為±2.63mGal；而當使用移去-恢復方法顧及地形影響時，模型2的檢核精度提高到±1.73mGal?？紤]到使用Poisson積分向上延拓解作為觀測量可能存在1～2mGal的數(shù)據(jù)誤差[44]，與新型航空重力測量系統(tǒng)精度水平基本相當[42]，對照表3的仿真檢驗結(jié)果，不難看出，表5給出的實測數(shù)據(jù)檢核結(jié)果完全符合預期，同時說明顧及地形效應對提高模型2的計算精度具有顯著成效。模型3的計算輸入量與航空重力測量觀測量無關，計算過程穩(wěn)定可靠，且模型自身就具備顧及地形效應的功能，因此具有較高的計算精度，達到±1.83mGal，與顧及地形效應后的模型2計算精度水平相當。

3　結(jié)　論

根據(jù)前面的模型適用性分析和數(shù)值計算比對結(jié)果，可得出以下基本結(jié)論：

(1) 基于逆Poisson積分的傳統(tǒng)向下延拓模型明顯受到數(shù)據(jù)分辨率、計算高度和觀測噪聲等多種因素的制約，即使在比較理想的數(shù)據(jù)精度條件下，最大延拓高度也只能達到1km。使用正則化方法可在一定程度上提高傳統(tǒng)延拓模型的穩(wěn)定性，但其解算過程仍存在較多的不確定性因素，不利于該模型的推廣應用。

(2) 基于頻譜截斷積分的直接向下延拓模型具有良好的計算穩(wěn)定性，對高頻觀測噪聲干擾具有很好的抑制作用。在當前航空重力測量通常采用的2′分辨率和3km飛行高度條件下，該模型幾乎可以不損失觀測數(shù)據(jù)精度的能力恢復地球表面的重力異常，因此具有良好的應用前景。在實際應用中需要注意把握的問題是積分核函數(shù)截斷階數(shù)與航空重力測量數(shù)據(jù)濾波截止頻率的匹配關系。解算結(jié)果存在一定的邊緣效應，同時觀測數(shù)據(jù)網(wǎng)格化處理給測線布設帶來的三維空間約束等因素是影響該模型推廣應用的主要障礙。

(3) 基于位模型和地形改正差分的直接向下延拓模型具有良好的應用前景，其計算過程完全獨立于航空重力觀測數(shù)據(jù)，因此不受觀測噪聲的干擾。其計算精度取決于超高階位模型和局部地形改正差分的相對精度，在現(xiàn)有技術(shù)條件下，該模型延拓計算精度可達2mGal，與當前航空重力測量精度水平相當。與模型2相比較，不需要對觀測數(shù)據(jù)作網(wǎng)格化預處理，也不存在積分邊緣效應，可實現(xiàn)空中和地面點對點延拓計算是該模型的主要優(yōu)勢。但關于該模型物理意義的嚴密解釋還有待進一步的研究，并通過更多的實際觀測數(shù)據(jù)計算來驗證其應用效果。

[1]李建成, 陳俊勇, 寧津生, 等. 地球重力場逼近理論與中國2000似大地水準面的確定[M]. 武漢: 武漢大學出版社, 2003.

LIJiancheng,CHENJunyong,NINGJinsheng,etal.TheoryoftheEarth’sGravityFieldApproximationandDeterminationofChinaQuasi-geoid2000[M].Wuhan:WuhanUniversityPress, 2003.

[2]黃謨濤, 翟國君, 管錚, 等. 海洋重力場測定及其應用[M]. 北京: 測繪出版社, 2005.

HUANGMotao,ZHAIGuojun,GUANZheng,etal.TheDeterminationandApplicationofMarineGravityField[M].Beijing:SurveyingandMappingPress, 2005.

[3]王彥飛. 反演問題的計算方法及其應用[M]. 北京: 高等教育出版社, 2007.WANGYanfei.ComputationalMethodsforInverseProblemsandTheirApplications[M].Beijing:HigherEducationPress, 2007.

[4]PELLINENLP.AccountingforTopographyintheCalculationofQuasigeoidalHeightsandPlumb-lineDeflectionsfromGravityAnomalies[J].BulletinGéodésique, 1962, 63(1): 57-65.

[5]BJERHAMMARA.GravityReductiontoaSphericalSurface[R].Stockholm:RoyalInstituteofTechnology, 1962.

[6]HEISKANENWA,MORITZH.PhysicalGeodesy[M].SanFrancisco:FreemanWHandCompany, 1967.

[7]JEKELIC.TheDownwardContinuationofAerialGravimetricDatawithoutDensityHypothesis[J].BulletinGéodésique, 1987, 61(4): 319-329.

[8]WANGYM.TheEffectofTopographyontheDeterminationoftheGeoidUsingAnalyticalDownwardContinuation[J].BulletinGéodésique, 1990, 64(3): 231-246.

[9]MARTINECZ.StabilityInvestigationsofaDiscreteDownwardContinuationProblemforGeoidDeterminationintheCanadianRockyMountains[J].JournalofGeodesy, 1996, 70(11): 805-828.

[10]XUPeiliang.TruncatedSVDMethodsforDiscreteLinearIll-posedProblems[J].GeophysicalJournalInternational, 1998, 135(2): 505-514.

[11]NOVKP,KERNM,SCHWARZKP.NumericalStudiesontheHarmonicDownwardContinuationofBand-limitedAirborneGravity[J].StudiaGeophysicaetGeodaetica, 2001, 45(4): 327-345.

[12]MUELLERF,MAYER-GUERRT.ComparisonofDownwardContinuationMethodsofAirborneGravimetryData[C]∥ProceedingsoftheInternationalAssociationofGeodesyIAGGeneralAssemblySapporo.Berlin:Springer, 2005, 128: 254-258.

[13]KERNM.AnAnalysisoftheCombinationandDownwardContinuationofSatellite,AirborneandTerrestrialGravityData[D].Calgary:UniversityofCalgary, 2003.

[14]ALBERTSB,KLEESR.AComparisonofMethodsfortheInversionofAirborneGravityData[J].JournalofGeodesy, 2004, 78(1-2): 55-65.

[15]TZIAVOSIN,ANDRITSANOSVD,FORSBERGR,etal.NumericalInvestigationofDownwardContinuationMethodsforAirborneGravityData[C]∥GGSM2004IAGInternationalSymposiumPorto.BerlinHeidelberg:Springer, 2005, 129: 119-124.

[16]ALBERTSB.RegionalGravityFieldModelingUsingAirborneGravimetryData[C]∥PublicationsonGeodesy70.Delft:NetherlandsGeodeticCommission, 2009.

[17]石磐, 王興濤. 利用航空重力測量和DEM確定地面重力場[J]. 測繪學報, 1997, 26(2): 117-121.

SHIPan,WANGXingtao.DeterminationoftheTerrainSurfaceGravityFieldUsingAirborneGravimetryandDEM[J].ActaGeodaeticaetCartographicaSinica, 1997, 26(2): 117-121.

[18]沈云中, 許厚澤. 不適定方程正則化算法的譜分解式[J]. 大地測量與地球動力學, 2002, 22(3): 10-14.

SHENYunzhong,XUHouze.SpectralDecompositionFormulaofRegularizationSolutionforIll-posedEquation[J].JournalofGeodesyandGeodynamics, 2002, 22(3): 10-14.

[19]王興濤, 石磐, 朱非洲. 航空重力測量數(shù)據(jù)向下延拓的正則化算法及其譜分解[J]. 測繪學報, 2004, 33(1): 33-38.

WANGXingtao,SHIPan,ZHUFeizhou.RegularizationMethodsandSpectralDecompositionfortheDownwardContinuationofAirborneGravityData[J].ActaGeodaeticaetCartographicaSinica, 2004, 33(1): 33-38.

[20]王振杰. 測量中不適定問題的正則化解法[M]. 北京: 科學出版社, 2006.

WANGZhenjie.RegularizationMethodsforIll-posedProbleminSurvey[M].Beijing:SciencePress, 2006.

[21]顧勇為, 歸慶明. 航空重力測量數(shù)據(jù)向下延拓基于信噪比的正則化方法的研究[J]. 測繪學報, 2010, 39(5): 458-464.

GUYongwei,GuiQingming.StudyofRegularizationBasedonSignal-to-noiseIndexinAirborneGravityDownwardtotheEarthSurface[J].ActaGeodaeticaetCartographicaSinica, 2010, 39(5): 458-464.

[22]蔣濤, 李建成, 王正濤, 等. 航空重力向下延拓病態(tài)問題的求解[J]. 測繪學報, 2011, 40(6): 684-689.JIANGTao,LIJiancheng,WANGZhengtao,etal.SolutionofIll-posedProbleminDownwardContinuationofAirborneGravity[J].ActaGeodaeticaetCartographicaSinica, 2011, 40(6): 684-689.

[23]鄧凱亮, 黃謨濤, 暴景陽, 等. 向下延拓航空重力數(shù)據(jù)的Tikhonov雙參數(shù)正則化法[J]. 測繪學報, 2011, 40(6): 690-696.

DENGKailiang,HUANGMotao,BAOJingyang,etal.TikhonovTwo-parameterRegulationAlgorithminDownwardContinuationofAirborneGravityData[J].ActaGeodaeticaetCartographicaSinica, 2011, 40(6): 690-696.

[24]蔣濤, 黨亞民, 章傳銀, 等. 基于矩諧分析的航空重力向下延拓[J]. 測繪學報, 2013, 42(4): 475-480.

JIANGTao,DANGYamin,ZHANGChuanyin,etal.DownwardContinuationofAirborneGravityDataBasedonRectangularHarmonicAnalysis[J].ActaGeodaeticaetCartographicaSinica, 2013, 42(4): 475-480.

[25]孫文, 吳曉平, 王慶賓, 等. 航空重力數(shù)據(jù)向下延拓的波數(shù)域迭代Tikhonov正則化方法[J]. 測繪學報, 2014, 43(6): 566-574.DOI: 10.13485/j.cnki.11-2089.2014.0080.

SUNWen,WUXiaoping,WANGQingbin,etal.WaveNumberDomainIterativeTikhonovRegularizationMethodforDownwardContinuationofAirborneGravityData[J].ActaGeodaeticaetCartographicaSinica, 2014, 43(6): 566-574.DOI: 10.13485/j.cnki.11-2089.2014.0080.

[26]劉曉剛, 李迎春, 肖云, 等. 重力與磁力測量數(shù)據(jù)向下延拓中最優(yōu)正則化參數(shù)確定方法[J]. 測繪學報, 2014, 43(9): 881-887.DOI: 10.13485/j.cnki.11-2089.2014.0160.

LIUXiaogang,LIYingchun,XIAOYun,etal.OptimalRegularizationParameterDeterminationMethodinDownwardContinuationofGravimetricandGeomagneticData[J].ActaGeodaeticaetCartographicaSinica, 2014, 43(9): 881-887.DOI: 10.13485/j.cnki.11-2089.2014.0160.

[27]黃謨濤, 歐陽永忠, 劉敏, 等. 海域航空重力測量數(shù)據(jù)向下延拓的實用方法[J]. 武漢大學學報(信息科學版), 2014, 39(10): 1147-1152.HUANGMotao,OUYANGYongzhong,LIUMin,etal.PracticalMethodsforDownwardContinuationofAirborneGravityDataintheSeaArea[J].GeomaticsandInformationScienceofWuhanUniversity, 2014, 39(10): 1147-1152.

[28]黃謨濤, 寧津生, 歐陽永忠, 等. 聯(lián)合使用位模型和地形信息的陸區(qū)航空重力向下延拓方法[J]. 測繪學報, 2015, 44(4): 355-362.DOI: 10.11947/j.AGCS.2015.20130751.HUANGMotao,NINGJinsheng,OUYANGYongzhong,etal.DownwardContinuationofAirborneGravimetryonLandUsingGeopotentialModelandTerrainInformation[J].ActaGeodaeticaetCartographicaSinica, 2015, 44(4): 355-362.DOI: 10.11947/j.AGCS.2015.20130751.

[29]NOVKP,HECKB.DownwardContinuationandGeoidDeterminationBasedonBand-limitedAirborneGravityData[J].JournalofGeodesy, 2002, 76(5): 269-278.

[30]MORITZH.AdvancedPhysicalGeodesy[M].Karlsruhe:HerbertWichmannVerlag, 1980.

[31]HWANGC,HSIAOYS,SHIHHC,etal.GeodeticandGeophysicalResultsfromaTaiwanAirborneGravitySurvey:DataReductionandAccuracyAssessment[J].JournalGeophysicalResearch, 2007, 112(B4): 93-101.

[32]FORSBERGR,OLESENAV.AirborneGravityFieldDetermination[M]∥XUGuochang.SciencesofGeodesy-I.Berlin:Springer-Verlag, 2010: 83-103.

[33]黃謨濤, 歐陽永忠, 翟國君, 等. 融合多源重力數(shù)據(jù)的Tikhonov正則化配置法[J]. 海洋測繪, 2013, 33(3): 6-12.

HUANGMotao,OUYANGYongzhong,ZHAIGuojun,etal.TikhonovRegularizationCollocationforMulti-sourceGravityDataFusionProcessing[J].HydrographicSurveyingandCharting, 2013, 33(3): 6-12.

[34]黃謨濤, 歐陽永忠, 劉敏, 等. 融合海域多源重力數(shù)據(jù)的正則化點質(zhì)量方法[J]. 武漢大學學報(信息科學版), 2015, 40(2): 170-175.

HUANGMotao,OUYANGYongzhong,LIUMin,etal.RegularizationofPoint-massModelforMulti-sourceGravityDataFusionProcessing[J].GeomaticsandInformationScienceofWuhanUniversity, 2015, 40(2): 170-175.

[35]NOVKP,KERNM,SCHWARZKP,etal.OnGeoidDeterminationfromAirborneGravity[J].JournalofGeodesy, 2003, 76(9-10): 510-522.

[36]PAULMK.AMethodofEvaluatingtheTruncationErrorCoefficientsforGeoidalHeight[J].BulletinGéodésique, 1973, 110(1): 413-425.

[37]陸仲連. 球諧函數(shù)[R]. 鄭州: 測繪學院, 1984.LUZhonglian.SphericalHarmonicFunction[R].Zhengzhou:SurveyingandMappingInstitute, 1984.

[38]OLESENAV.ImprovedAirborneScalarGravimetryforRegionalGravityFieldMappingandGeoidDetermination[D].Copenhagen:UniversityofCopenhagen, 2002.

[39]孫中苗. 航空重力測量理論、方法及應用研究[D]. 鄭州:信息工程大學, 2004.

SUNZhongmiao.Theory,MethodsandApplicationsofAirborneGravimetry[D].Zhengzhou:InformationEngineeringUniversity, 2004.

[40]PAVLISNK,HOLMESSA,KENYONSC,etal.TheDevelopmentandEvaluationoftheEarthGravitationalModel2008 (EGM2008)[J].JournalGeophysicalResearch, 2012, 117(B4): 531-535.

[41]章傳銀, 郭春喜, 陳俊勇, 等.EGM2008地球重力場模型在中國大陸適用性分析[J]. 測繪學報, 2009, 38(4): 283-289.

ZHANGChuanyin,GUOChunxi,CHENJunyong,etal.EGM2008andItsApplicationAnalysisinChineseMainland[J].ActaGeodaeticaetCartographicaSinica, 2009, 38(4): 283-289.

[42]歐陽永忠, 鄧凱亮, 陸秀平, 等. 多型航空重力儀同機測試及其數(shù)據(jù)分析[J]. 海洋測繪, 2013, 33(4): 6-11.

OUYANGYongzhong,DENGKailiang,LUXiuping,etal.TestsofMulti-typeAirborneGravimetersandDataAnalysis[J].HydrographicSurveyingandCharting, 2013, 33(4): 6-11.

[43]MORITZH.TheComputationoftheExternalGravityFieldandtheGeodeticBoundary-ValueProblem[C]∥ProceedingsoftheSymposiumontheExtensionofGravityAnomaliestoUnsurveyedAreas.GeophysicalMonographSeriesNumber9. [S.l.]:AmericanGeophysicalUnionoftheNationalAcademyofSciences, 1966: 127-136.

[44]ARGESEANUVS.UpwardContinuationofSurfaceGravityAnomalyData[C]∥ProceedingsoftheIAGSymposiumonAirborneGravityFieldDetermination.Boulder: [s.n.], 1995: 95-102.

(責任編輯：陳品馨)

TestandAnalysisofDownwardContinuationModelsforAirborneGravityDatawithRegardtotheEffectofTopographicHeight

LIUMin1,HUANGMotao2,3，OUYANGYongzhong2，DENGKailiang2,ZHAIGuojun2,3，WUTaiqi2

1.InstituteofGeospacialInformation,InformationEngineeringUniversity,Zhengzhou450001,China； 2.NavalInstituteofHydrographicSurveyingandCharting，Tianjin300061，China； 3.DepartmentofNavigation,NavalUniversityofEngineering,Wuhan430033,China

Downwardcontinuationisanessentialtechnicalstepofdataprocessinginairbornegravimetryforfurtherapplications.Itisknownthatthesolutionofdownwardcontinuationisuncertainduetoitsill-posedness.Soithasbeenatopicofgeneralinterestformanyscholarsathomeandabroadingeodesy.Themainpurposeofthispaperistogive3representativemodelsfordownwardcontinuationincludingtraditionalinversePoissonintegrationandtwomodernmethods,andmakeacomprehensivecomparisonandanalysisontheirpropertyandapplicabilityamongthedifferentmodels.Ultra-high-degreegeopotentialmodel,localtopographiccorrectionandremove-restoretechniquearesuggestedtobeusedforregardtotheeffectoftopographicheight,andfortherealizationofdownwardcontinuationcombiningwithatransformationfromsphericaltoundulatingsurface.Wepayourattentiontotheinfluenceofsurveyeddataerrorsonthestabilityofdownwardcontinuationsolutions.Theoreticalanalysis,simulateddataandrealnumericalcomputationsarecarriedouttoevaluatethestabilityandaccuracyofdownwardcontinuationmodels.Andsomeusefulconclusionsareobtained.Underexistingworkingconditions,thetraditionalinversePoissonintegrationmethodcanonlybeusedtothecontinuationcomputationunder1kmduetotheseriousdisturbingofsurveyingnoise.Excellentcomputationstabilitycanbeachievedbyusingtheband-limitedspectrumandthegeopotentialmodelplustopographiccorrectionmethods.Thetwonewmodelscanbeusedtothedownwardcontinuationofairbornegravitydataon5kmheightand2′dataresolution.Andtheaccuracyofcorrespondingcontinuationsolutionscanbereach2×10-5m/s2.Itcanmeettherequirementsfromdifferentapplications.

airbornegravimetry;downwardcontinuation;terraineffect;Poissonintegration;band-limitedspectrum;ultra-high-degreegeopotentialmodel

LIUMin,HUANGMotao，OUYANGYongzhong，etal.TestandAnalysisofDownwardContinuationModelsforAirborneGravityDatawithRegardtotheEffectofTopographicHeight[J].ActaGeodaeticaetCartographicaSinica,2016,45(5):521-530.DOI:10.11947/j.AGCS.2016.20150453.

P223

A

1001-1595(2016)05-0521-10

國家973計劃(613219)；國家自然科學基金(41474012；41174062；41374018)；國家重大科學儀器設備開發(fā)專項(2011YQ12004503)

引文格式：劉敏，黃謨濤，歐陽永忠，等.顧及地形效應的重力向下延拓模型分析與檢驗[J].測繪學報，2016,45(5)：521-530.DOI:10.11947/j.AGCS.2016.20150453.