測高重力異常中央區奇異積分數值求積公式

2021-11-25 09:59:38宗敬文李厚樸歐陽永忠

測繪學報 2021年10期

宗敬文，李厚樸，紀兵，歐陽永忠

1. 海軍工程大學導航工程教研室，湖北武漢 430033； 2. 自然資源部海洋環境探測技術與應用重點實驗室，廣東廣州 510300

海洋重力場數據的獲取手段有衛星測高反演、船載及航空重力測量，隨著人造地球軌道衛星技術的發展，衛星測高技術促進了大地測量學、地球物理學和海洋學等學科的交叉發展[1]，同時與傳統的海洋觀測手段相比，能夠周期性獲取除極地以外的高分辨率、全天候、長時間序列的全球海洋觀測數據[2-5]。截至目前，世界上各個國家共計發射測高衛星18顆，其中包括美國海軍的Geosat系列、美國宇航局和法國空間研究中心聯合發射的T/P系列[6]、歐空局的ERS系列[7]及中國的HY-2A[8]。衛星測高數據遍布全球85°S—85°N的海洋區域，分辨率達1′×1′，精度達1～2 cm，利用高精度、高分辨率衛星測高數據反演重力異常精度已達到10 mGal(1 Gal=10-2m/s2)以內。文獻[9]使用逆Stokes法反演中國近海30′×30′海洋重力異常，精度達到3.5 mGal。文獻[10]使用Laplace方程的垂線偏差法反演全球海域的重力異常，反演的全球海域1′×1′重力場模型與NGDC船測數據的檢核精度達到4～8 mGal。文獻[11—14]在衛星測高數據融合和精細處理方面也作出了大量的貢獻。

文獻[15]利用Stokes公式將衛星測高數據轉換為重力異常。文獻[16]基于垂線偏差對地形信息敏感的特點，根據邊值理論由重力與地形數據確定格網垂線偏差模型。文獻[17]提出了利用垂線偏差來反演重力異常的方法，并系統地給出了垂線偏差、擾動重力和重力異常之間的關系表達式，但在利用大地水準面高和垂線偏差反演重力異常實際計算中，計算點及其附近區域到計算點的理論距離近似于零，會導致逆Stokes公式和逆Vening-Meinesz公式中的積分函數產生奇異。文獻[18—19]提出一組“非奇異變換”，系統地解決了物理大地測量中高程異常、垂線偏差、地形改正及重力異常梯度等泛函的中央區計算問題。對于中央區積分問題的研究，文獻[20—21]將該積分區域視為圓形域，推導的逆Vening-Meinesz公式被廣泛應用于反演海洋重力異常。文獻[22]在利用逆Vening-Meinesz公式的FFT算法時，考慮了中央區效應的影響，研究表明中央區效應對重力異常的影響可以達到百微伽量級，在高精度重力場計算中是不能忽視的。文獻[23—24]分別推導了逆Stokes法和逆Vening-Meinesz法的衛星測高反演重力異常解析計算公式，解決了奇異積分問題，具有較高的計算精度，但解析計算公式系數繁多，計算量大，實際應用中有所不便。因此，為進一步完善衛星測高反演重力異常的計算方法，簡化計算過程，提高計算效率，本文采用數值求積公式，分別利用Simpson公式和Cotes公式對逆Stokes法和逆Vening-Meinesz法中的奇異積分問題進行了新的研究，系統地推導出了中央區重力異常的計算公式。此公式可直接利用格網節點處的大地水準面高和垂線偏差計算重力異常值，形式簡單，計算效率高，計算精度與解析法計算結果精度相當，可以滿足實際應用。

1 逆Stokes公式法奇異積分數值求積公式

由文獻[23—24]可知，大地水準面計算重力異常的逆Stokes公式平面近似形式為

(1)

(2)

假設中央區NQ可以展開為如下Taylor級數

NQ=N(x,y)=NP+Nxx+Nyy+

(3)

由式(3)可得

(4)

設中央區積分面積為σ∈[-a≤x≤a,-b≤y≤b]，如圖1所示。為方便推導計算，引入新的積分變量u=x/a和v=y/b，因此中央區積分區域面積為σ′∈[-1≤u≤1,-1≤v≤1]。

(5)

圖1 中央區示意圖Fig.1 Sketch map of the innermost area

對σ1引入新的積分變量

(6)

對σ2引入新的積分變量

(7)

將式(6)和式(7)代入式(5)，可得中央區IN的表達式為

(8)

1.1 基于Simpson公式推導出的逆Stokes公式中央區數值求積公式

利用Simpson公式

(9)

考慮到式(4)，利用式(9)對式(8)進行u和v方向上數值積分后可得

(10)

式中，INS表示基于Simpson公式的逆Stokes法奇異積分數值求積公式。將N(u,v)在u,v=0,±1處的值記為Nij(i,j=0,±1)，并對式(10)在k和λ方向上繼續使用Simpson公式，其中

(11)

根據數值微分公式可由節點處大地水準面高求出

(12)

整理后可得

(13)

將式(13)代入式(2)，最后可得中央區重力異常的數值求積公式為

Nab-4NP)

]

(14)

1.2 基于Cotes公式推導出的逆Stokes公式中央區數值求積公式

當只考慮中央區范圍為σ∈[-a≤x≤a,-b≤y≤b]時，中央區的影響是不完整的，特別是在起伏變化較大的區域時，考慮將中央區的范圍擴大為σ∈[-2a≤x≤2a,-2b≤y≤2b]，因此，本文引入了Cotes公式來計算。推導過程類似于上一部分，首先利用Cotes公式

(15)

考慮到式(4)，利用式(15)對式(8)進行u和v方向上數值積分后可得

(16)

式中，INC表示基于Cotes公式的逆Stokes法奇異積分求積公式。

如圖2所示，將N(u,v)在u,v=0,±1,±2處的值記為Nij(i,j=0,±1,±2)，對式(16)在k和λ方向上繼續使用Cotes公式，并考慮到式(11)和式(12)，整理后可得

Na-2b+N-a2b+N-a-2b-8NP)

}

(17)

將式(17)代入式(2)中，最后可得中央區重力異常的數值求積公式為

N-2a-2b+N-2a2b+N2a-2b-4NP)+

N-2a-b+Na2b+Na-2b+N-a2b+

N-a-2b-8NP)

}

(18)

圖2 Cotes公式中央區示意圖Fig.2 Sketch map of the innermost area used the Cotes formula

利用式(14)和式(18)計算中央區重力異常時，可直接使用已知網格節點處的大地水準面高作為輸入值，相較于利用解析法的求解過程減少了求解多項式系數的步驟，大大簡化了逆Stokes法反演重力異常的過程，從而提高了計算效率。

2 逆Vening-Meinesz公式法奇異積分數值求積公式

由文獻[25—26]可知，垂線偏差計算重力異常的逆Vening-Meinesz公式平面近似形式為

(19)

式中，ξQ、ηQ為移動點處的垂線偏差；ξP、ηP為計算點處的垂線偏差。為便于計算，令

(20)

將垂線偏差ξQ、ηQ在中央區展成如下Taylor級數

(21)

由式(21)得出

(22)

由于積分區域的對稱性，并考慮到式(6)和式(7)，式(20)可改寫為

(23)

2.1 基于Simpson公式推導出的逆Vening-Meinesz公式中央區數值求積公式

利用式(9)，并考慮到式(22)，對式(23)中積分部分進行u和v方向上數值積分后可得

(24)

將ξ(u,v)、η(u,v)在u、v=0、±1處的值記為ξij、ηij(i,j=0、±1)，并考慮到根據數值微分公式可由節點處垂線偏差求出

(25)

最后對式(24)在k和λ方向上繼續使用Simpson公式，整理后得

(26)

式中，ΔgξS、ΔgηS表示基于Simpson公式的逆Vening-Meinesz法奇異積分數值求積公式。因此，式(26)即為基于Simpson公式的逆Vening-Meinesz法反演中央區重力異常的數值求積公式。

2.2 基于Cotes公式推導出的逆Vening-Meinesz公式中央區數值求積公式

同樣，當只考慮中央區范圍為σ∈[-a≤x≤a,-b≤y≤b]時，中央區的影響是不完整的，特別是在起伏變化較大的區域時，考慮將中央區的范圍擴大為σ∈[-2a≤x≤2a,-2b≤y≤2b]，因此，本文引入了Cotes公式來計算。利用式(15)并考慮到式(22)，對式(23)中積分部分進行u和v方向上數值積分后可得

(27)

如圖2所示，將ξ(u,v)、η(u,v)在u,v=0,±1,±2處的值記為ξij、ηij(i,j=0,±1,±2)。對式(27)在k和λ方向上繼續使用Cotes公式，并考慮到式(11)和式(25)，整理后可得

(28)

式中，ΔgξC、ΔgηC表示基于Cotes公式的逆Vening-Meinesz法奇異積分數值求積公式，因此，式(28)即為基于Cotes公式的逆Vening-Meinesz法反演中央區重力異常數值求積公式。利用式(26)和式(28)計算中央區重力異常時，可直接使用已知網格節點處的垂線偏差作為輸入值，相較于利用解析法的求解過程減少了求解插值多項式系數的步驟，從而簡化了逆Vening-Meinesz公式反演重力異常的計算過程，提高了計算效率。

3 算例分析

3.1 理論模型檢驗

為驗證本文推導的基于數值求積公式的測高反演重力異常普適計算公式的準確性和可靠性，假定大地水準面高和垂線偏差為理論模型值，計算了非奇異變換前、后求積公式計算結果，分析比較了非奇異變換后兩種不同普適數值求積公式的計算精度?，F設計如下兩個算例。

3.1.1 理論模型1

取中央區大地水準面高的數學模型為

(29)

中央區垂線偏差的數學模型為

(30)

根據本文推導出的普適數值求積公式取a=1,b=1，因此中央區范圍為σ∈[-2≤x≤2,-2≤y≤2]，δ為平滑因子，則NP=NQ(0,0)=δ，ξP=ξQ(0,0)=0，ηP=ηQ(0,0)=0。以δ取10時為例，大地水準面高如圖3所示。

圖3 大地水準面高Fig.3 Geoidal height

(31)

為驗證非奇異變換前、后積分表達式精度，分別取δ=1、5、10時，IN1和IN2在不同等分區間下的計算結果列于表1。

由表1可知，非奇異變換后IN250×50等分計算結果的準確度和非奇異變換前IN1500×500等分計算結果的準確度基本相同。由圖4(a)可以看出，非奇異變換前在低等分區間精度較低且收斂速度慢，非奇異變換后的收斂速度要大于變換前且計算結果更接近真值。因此，非奇異變換后的積分公式計算效率更高，收斂速度更快且精度更高。為進一步分析比較基于Simpson公式和Cotes公式的逆Stokes法反演中央區重力異常普適數值計算公式，利用大地水準面高理論模型1，取IN24000×4000等分的計算結果為準確值，式(13)和式(17)的計算結果分別記為INS和INC，將δ分別取δ=3,5,8,10時的計算結果列于表2。

表1 非奇異變換前后奇異積分計算結果

圖4 δ=10時非奇異變換前后不同等分區間積分數值Fig.4 Results of singular integral before and after the non-singular transformation

由表2可知，當δ=3，5時，Cotes公式法的相對誤差小于Simpson公式法；當δ=8，10時Simpson公式法的相對誤差要小于Cotes公式法。Simpson公式法數值計算公式的相對誤差波動較大，與平滑因子δ取值有關，大地水準面高變化緩和地區相對誤差較小，精度較高；Cotes公式法數值計算公式的相對誤差比較穩定，隨著平滑因子δ取值變化，精度受影響較小，但精度稍低。相比較于Cotes公式法，Simpson公式法在平緩地區具有更好的計算精度，但其計算精度波動較大；相比較于Simpson公式法，Cotes公式法具有更好的計算穩定性，計算精度受地形影響較小。

(32)

為驗證非奇異變換前、后積分表達式精度，分別取δ=1、5、10，IV1和IV2在不同等分區間下的計算結果列于表3。

表2 非奇異變換后數值解計算結果誤差

表3 非奇異變換前后奇異積分計算結果

表4 非奇異變換后數值解計算結果誤差

由表4可知，Cotes公式法數值計算公式的相對誤差明顯小于Simpson公式法數值計算公式，Simpson公式法數值計算公式的相對誤差波動較大，與平滑因子δ取值有關；Cotes公式法數值計算公式的相對誤差比較穩定，當δ分別取δ=3,6,8,10時，相對誤差均小于1%，相比較于Simpson公式法數值計算公式，Cotes公式法數值計算公式具有較好的計算精度和穩定性。

文獻[23]和文獻[24]分別推導了雙二次多項式插值表示下的大地水準面高和垂線偏差非奇異積分解析計算公式。記錄其結果分別為IN3和IV3

(α02+α20))

(33)

(34)

式中,IN3為雙二次多項式插值表示下的大地水準面高非奇異積分解析表達式；IV3為雙二次多項式插值表示下的垂線偏差非奇異積分解析表達式。以非奇異變換后4000×4000等分區間計算結果作為準確值，其中逆Stokes法準確值記為IN2，逆Vening-Meinesz法準確值記為IV2，式(13)、式(17)、式(26)和式(28)的計算結果分別為INS、INC、IVS和IVC，計算結果及其相應計算誤差見表5。

表5 非奇異變換后解析解與數值解的計算誤差

分析表5可知，逆Stokes法重力異常求積公式解析解與準確值的相對誤差為0.64%，Simpson公式數值解與準確值的相對誤差為0.15%，Simpson公式數值解精度高于解析解計算精度，二者相對誤差均小于1%，而Cotes公式數值解計算精度較低，相對誤差為1.62%；逆Vening-Meinesz法重力異常求積公式解析解與準確值的相對誤差為0.62%，Cotes公式數值解與準確值的相對誤差為0.46%，Simpson公式數值解與準確值的相對誤差為0.50%，解析解計算精度略低于數值解計算精度，三者相對誤差均小于1%。

3.1.2 理論模型2

在理論模型1的分析基礎上，理論模型2僅用于分析兩種不同數值積分方法下測高反演重力異常普適計算公式的準確性和可靠性。取中央區大地水準面高的數學模型為

(35)

中央區垂線偏差的數學模型為

(36)

現取a=1,b=0.7，因此中央區范圍為σ∈[-2≤x≤2,-1.4≤y≤1.4]，δ為平滑因子，則NP=NQ(0,0)=δ，ξP=ξQ(0,0)=0，ηP=ηQ(0,0)=0。以δ取10時為例，大地水準面高如圖5所示。

圖5 大地水準面高Fig.5 Geoid height

(37)

(38)

為分析比較基于Simpson公式法和Cotes公式法的逆Stokes法反演中央區重力異常普適數值計算公式，利用大地水準面高理論模型2，式(13)和式(17)的計算結果分別記為INS和INC，將δ分別取δ=3，5，8，10時的計算結果列于表6。

表6 非奇異變換后數值解計算結果誤差

由表6可知，當δ=3時，Cotes公式法的相對誤差小于Simpson公式法；當δ=5、8、10時，Simpson公式法的相對誤差小于Cotes公式法。Simpson公式法數值計算公式的相對誤差波動較大，當平滑因子變大時，計算精度較高，相對誤差小于0.5%；Cotes公式法數值計算公式的相對誤差比較穩定，隨著平滑因子δ取值變化，精度受影響較小，但精度稍低。相比較于Cotes公式法，Simpson公式法在平緩地區具有更好的計算精度，但其計算精度波動較大；相比較于Simpson公式法，Cotes公式法具有更好的計算穩定性，計算精度受地形影響較小。

利用垂線偏差理論模型2，式(26)和式(28)的計算結果分別為IVS和IVC，將δ分別取δ=3,5,8,10時的計算結果列于表7。

表7 非奇異變換后數值解計算結果誤差

由表7可知，Cotes公式法數值計算公式的相對誤差比較穩定，當δ分別取δ=3,5,8,10時，計算結果的相對誤差均小于0.5%；Simpson公式法數值計算公式的相對誤差波動較大，且當δ=3,5時，相對誤差較大，當δ=8,10時，相對誤差小于0.5%。對比發現，Cotes公式法數值計算公式相比較于Simpson公式法數值計算公式具有較好的計算精度和計算穩定性。

將式(33)和式(34)中解析解計算結果與本文推導公式進行比較，取δ=10時，以式(37)、式(38)非奇異變換后4000×4000等分區間計算結果作為準確值，其中逆Stokes法準確值記為IN22，逆Vening-Meinesz法準確值記為IV22，IN3、INS、INC、IV3、IVS和IVC計算結果及其相應計算誤差見表8。

表8 非奇異變換后解析解與數值解的計算誤差

分析表8可知，其中逆Stokes法重力異常數值求積公式解析解與準確值的相對誤差為0.57%，Simpson公式法數值解與準確值的相對誤差為0.15%，Simpson公式法數值解精度高于解析解計算精度，二者相對誤差均小于0.6%，而Cotes公式數值解計算精度較低，相對誤差為2.20%；逆Vening-Meinesz法重力異常數值求積公式解析解與準確值的相對誤差為0.44%，Cotes公式法數值解與準確值的相對誤差為0.02%，Simpson公式法數值解與準確值的相對誤差為0.31%，數值解計算精度略高于解析解計算精度。

本文通過利用兩種不同理論模型，分析比較基于Simpson公式法和Cotes公式法兩種不同數值求積方法推導出的測高反演重力異常普適計算公式，發現在逆Stokes法重力異常求積公式計算中，Simpson公式法計算精度較好，但受地形影響其誤差波動較大，適用于中央區變化比較平緩地區，Cotes公式法計算誤差穩定性較高，但計算精度略低；在逆Vening-Meinesz法重力異常求積公式中，Cotes公式法計算精度和計算誤差穩定性都優于Simpson公式法，且計算精度優于解析法計算精度，可以滿足應用要求。

3.2 試驗數據檢驗

理論模型下的精度檢驗表明，本文推導出的中央區奇異積分數值求積公式有較高的計算精度且形式簡單。為進一步說明中央區奇異積分數值求積公式的計算精度及中央區對重力異常的影響，選定南海16°N-20°N，112°E-116°E海域作為計算區域，以EGM2008地球重力場模型計算得到的2′×2′分辨率的大地水準面高數據作為試驗數據，計算海域大地水準面高分布圖如圖6所示，計算海域共記119×119個格網，中央區為8′×8′(4×4個格網)，采用逆Stokes法利用式(33)、式(14)和式(18)實際計算了該區域的中央區重力異常，結果分別記為Δg1、Δg2、Δg3列于表9，計算結果之間的比較情況見表10。

由表9和表10可以看出，3種計算方法計算出的4×4個格網大小的中央區重力異常平均值分別為-1.088、-1.089、-1.087 mGal，三者差異較?。粯藴什罘謩e為2.60、2.63、2.83 mGal；Simpson公式法與解析法計算得到的中央區重力異常差值的標準差為0.030 mGal，最大值為0.195 mGal，最小值為-0.129 mGal,因此當考慮中央區范圍擴大為4×4個格網時，Simpson公式法和解析法計算結果相當；Cotes公式法計算出的重力異常最大值為14.39 mGal，最小值為-13.23 mGal，與Simpson公式法和解析法計算得到的中央區重力異常差值的標準差分別為0.533 mGal和0.545 mGal，最大值分別為3.806 mGal和4.254 mGal，最小值分別為-4.167 mGal和-3.879 mGal,Cotes公式法與另外兩種方法計算結果具有一定差距，這是由于Cotes公式法計算過程中計算點更多，反映出利用該方法計算中央區重力異常變化更為明顯,可以更好地反映出起伏較大區域中央區重力異常變化。同時，試驗數據說明，分別利用3種不同方法計算出的4×4個格網中央區重力異常范圍為-13.5 mGal～+14.5 mGal，因此在高精度重力異常反演時，中央區有不容忽視的貢獻。