基于調(diào)和分析法的全球地形球諧系數(shù)模型構(gòu)建

單建晨，李姍姍，范 雕，李新星，黃 炎，邢志斌

（1.信息工程大學(xué)，河南 鄭州 450051；2.航天工程大學(xué)，北京 102200）

地形球諧系數(shù)模型是指描述地球地形地貌的一組基本參數(shù)的集合，是對全球地形數(shù)據(jù)的解析表達(dá)。相較于數(shù)值模型，地形球諧系數(shù)模型具有可以表達(dá)全球地形頻譜特性的優(yōu)點，并且依托地球物理公式關(guān)系便于實現(xiàn)與重力位、擾動位[1]和地球密度的轉(zhuǎn)換運用，同時避免了在獲取特定區(qū)域、點位物理量時計算積分、內(nèi)插等復(fù)雜過程。目前地形球諧系數(shù)模型在地形正向建模、布格異?；蛑亓_動計算、球諧地形引力位模型構(gòu)建和未測區(qū)域地形數(shù)據(jù)預(yù)測等方面擁有實際運用及廣闊前景[2,3]。

通常球諧系數(shù)模型的求解方法可分為調(diào)和分析法（Spherical Harmonic Analysis，SHA）和最小二乘法（Least-squares，LSQ）。最小二乘法能夠提供求算球諧系數(shù)的參數(shù)方差、協(xié)方差信息，但解算甚（超）高階球諧系數(shù)時，需要大量的計算資源對超大型非對角法方程矩陣進行求逆，對計算機性能提出了非常高的要求。調(diào)和分析法的出發(fā)點是對積分計算的數(shù)值評估，計算易于實現(xiàn)，但因為離散勒讓德函數(shù)的非正交性使得精度欠佳，為了提高調(diào)和分析法求解模型精度國內(nèi)外學(xué)者做了大量研究[4,5]。

本文研究推導(dǎo)了現(xiàn)有調(diào)和分析方法，以國際發(fā)布地形球諧系數(shù)模型Earth2012 球諧綜合計算的地形值作為輸入數(shù)據(jù)，分別利用Gauss-Legendre 積分法和矩形離散積分法完成了2160 階地形球諧系數(shù)模型的構(gòu)建，探討分析了各個方法的優(yōu)劣，著重分析了不同積分法恢復(fù)的模型在高緯度區(qū)域的表達(dá)精度，并提出了聯(lián)合快速傅里葉變換（FFT）技術(shù)和基于OpenMP 多核并行技術(shù)的模型系數(shù)計算策略，大大提高了計算效率。

1 原理與方法

1.1 矩形離散積分法

根據(jù)展開理論，一個任意的球面上的函數(shù)可以展開為面球諧函數(shù)，因此全球地形可由球諧系數(shù)表達(dá)成關(guān)于球心經(jīng)緯度的完全正規(guī)諧函數(shù)形式：

其中f(θ,λ)為全球地形函數(shù)，θ為以北極方向為起始方向的極距，范圍為(0,π)，λ為球心經(jīng)度，為歸一化締和Legendre 函數(shù)，、為全球地形完全正規(guī)諧函數(shù)系數(shù)。

利用正交關(guān)系，可以依據(jù)式(1)得到、的表達(dá)式：

通?，F(xiàn)實情況下，f(θ,)λ表達(dá)為等間隔的離散數(shù)值模型，與式(2)中數(shù)值積分計算的連續(xù)函數(shù)要求相悖。調(diào)和分析法利用f(θ,)λ離散值代入求和逼近數(shù)值積分，從而避免了獲取連續(xù)全球地形函數(shù)的難題。通過調(diào)和分析法所得面球諧系數(shù)存在最大（截斷）階次nmax，nmax與輸入的全球地形數(shù)據(jù)分辨率R存在關(guān)系nmax=180°/R。

一種計算處理式(2)的方法是將關(guān)于球心經(jīng)緯度的二重積分轉(zhuǎn)換成全球區(qū)域等角格網(wǎng)數(shù)據(jù)的離散求和計算：

其中，n=0,1…nmax，m=0,1…n；、為格網(wǎng)中心點的極距和球心經(jīng)度；Δθ=π/nmax，Δλ=2π/2nmax。

1.2 Gauss-Legendre 積分法

由于離散勒讓德函數(shù)的非（弱）正交性，按照矩形離散積分方法處理積分式將會導(dǎo)致較大的精度損失，相比之下，Gauss-Legendre 積分法通過設(shè)置Gauss-Legendre 權(quán)重保證了從連續(xù)函數(shù)過渡到離散函數(shù)中勒讓德函數(shù)正交性。該方法引入數(shù)值分析中Gauss-Legendre 積分的概念對數(shù)據(jù)進行重新采樣、賦權(quán)并求和。區(qū)別于上文方法中的等間距網(wǎng)格（nmax×2nmax）的構(gòu)建，Gauss-Legendre 積分法依據(jù)最大（截斷）階數(shù)nmax建立高斯格網(wǎng)（（nmax+1）×（2nmax+1）），該格網(wǎng)的緯度圈為滿足方程Pnmax+1(cosθi)=0（θi為極矩）的Gauss-Legendre 節(jié)點，格網(wǎng)子午圈則等間距分布；同時為了便于編程實現(xiàn)對離散求和的公式進行了改進，將式(2)中的面球諧系數(shù)展開至經(jīng)緯度方向分步拆解計算，對于經(jīng)度方向：

其中m=0,1,2…,∞。應(yīng)用離散求和法計算式(4)中定積分近似值：

其中m=0,1,2…,nmax。對于緯度方向：

其中n=m,m+1,m+2…,∞。運用數(shù)值分析中Gauss-Legendre 方法對式(6)離散求積：

其中，n=m,m+1,m+2…∞；式(5)(7)中δm0為克羅內(nèi)克爾符號；式(7)中wi為Gauss-Legendre 權(quán)，定義為：

其中i=0,1,2…nmax。需要特別指出的是式(8)分母中(cosθ)為勒讓德多項式對cosθ的一階導(dǎo)數(shù)，并非對θ的一階導(dǎo)數(shù)。

2 試驗結(jié)果與分析

2.1 數(shù)據(jù)準(zhǔn)備和處理

本文選取科廷大學(xué)發(fā)布的Earth2012.topo_bathy（Earth2012 系列模型之一）2160 階地形球諧系數(shù)模型作為試驗數(shù)據(jù)，該球諧模型是基于SRTM V4.1、SRTM30_PLUS、Etopo1 等數(shù)據(jù)源所構(gòu)建的5 ′×5 ′數(shù)值模型并做球諧分析得到。為保證用于Gauss-Legendre 和矩形離散積分法兩種球諧分析的輸入數(shù)據(jù)具有一致性，試驗利用Earth2012.topo_bathy模型球諧綜合完成所構(gòu)建格網(wǎng)的數(shù)據(jù)采樣，并將Earth2012.topo_bathy 模型作為真值，以檢驗兩種方法所得地形球諧系數(shù)的還原效果。試驗采用改進的Belikov 方法計算Legendre 函數(shù)，可以較好地保證函數(shù)計算精度與穩(wěn)定性[6]。在編程實現(xiàn)球諧綜合計算時，使用OpenMP 并行算法和FFT 技術(shù)，實現(xiàn)了快速運算。

2.2 Gauss-Legendre 格網(wǎng)緯度圈計算

根據(jù)Gauss-Legendre 格網(wǎng)構(gòu)建要求，格網(wǎng)緯度圈取值為函數(shù)(x)在區(qū)間x∈[-1,1]的零點（nmax+1個）。試驗使用搜索迭代法（scan-iteration method,SIM）求算方程解[7]。因為函數(shù)(x)關(guān)于原點對稱，為減少計算量只需要考慮求算函數(shù)在區(qū)間x∈(0,1]的零點。SIM 計算步驟如下：

（1）設(shè)定區(qū)間寬度w，根據(jù)給定區(qū)間寬度，通過區(qū)間端點函數(shù)值的異號性判斷是否含有零點，并標(biāo)記含零點區(qū)間；

（2）不斷縮小搜索區(qū)間寬度，直至 Int[nmax/2]個零點獨立分布在標(biāo)記區(qū)間；

（3）在每一個標(biāo)記區(qū)間使用割線迭代法計算零點值。

2.3 計算效率

綜合上述式(5)(7)，Gauss-Legendre 積分法計算球諧系數(shù)的實質(zhì)運算公式為：

其中，n=0,1…nmax、m=0,1…n。若直接按照式(9)累加的方式串行編譯實現(xiàn)球諧系數(shù)的求解，計算效率低下，計算耗時約為121.82h，難以適應(yīng)應(yīng)用需求，實際試驗中分別采用FFT 技術(shù)和OpenMP 并行算法對式(5)(7)中的Am(θ)、Bm(θ)與Cnm、Snm的計算求解。

式(5)為卷積表達(dá)形式，可以按照式(10)對Am(θ)、Bm(θ)進行快速傅里葉變換處理，其中Re、Im 分別代表變換所得實部與虛部[8,9]。

其中m=0,1,2…,nmax。考慮式(7)中各階次待求球諧系數(shù)僅僅與唯一變量極矩θ相關(guān)，因而考慮以極矩θ為單元并行計算，此種并行方式還可以將勒讓德函數(shù)計算次數(shù)降到最低，即一個采樣緯度圈上的所有采樣點只需計算一次勒讓德函數(shù)。為避免內(nèi)存位置讀寫沖突，利用pravite()子句對簡單變量實現(xiàn)線程私有化處理，對勒讓德函數(shù)[n][m]、s球諧系數(shù)Cnm[n][m]、Snm[n][m]等數(shù)組升維成[v][n][m]、Cnm[v][n][m]、Snm[v][n][m]，通過對v的控制，實現(xiàn)各線程中對這些數(shù)組的獨立訪問與運算[10]。

試驗所用計算器處理器為Intel?Core? i7-9750H CPU @ 2.60 GHz，內(nèi)存32.00 GB，采用VS2013 編程平臺10 個線程并行計算，總耗時466.083 s，相比于串行累加算法耗時，加速比約為983，較好地解決了串行編譯耗時多、效率低的問題。

2.4 精度評價

按照上文原理算法，建立矩形等角格網(wǎng)和Gauss-Legendre 格網(wǎng)，并利用格網(wǎng)數(shù)據(jù)恢復(fù)地形球諧模型。其中Gauss-Legendre 積分法恢復(fù)模型記為Sph.GL 模型，矩形離散積分法恢復(fù)模型記為Sph.Rect模型，試驗以原模型Earth2012 為標(biāo)準(zhǔn)，引入模型系數(shù)階方差、系數(shù)階誤差等標(biāo)準(zhǔn)檢驗?zāi)Ｐ途取?/p>

以Earth2012 作為標(biāo)準(zhǔn)模型，計算兩個恢復(fù)模型的階誤差RMSn，公式如下：

其中，n=0,1,2…nmax；、為Earth2012模型球諧系數(shù)；、為恢復(fù)模型球諧系數(shù)。階誤差反應(yīng)了恢復(fù)模型與標(biāo)準(zhǔn)模型之間的差異，差值越小，說明模型符合程度越高，恢復(fù)效果更好。

如圖1所示用階誤差RMSn對模型進行對比，為便于直觀比較，同時繪制了Earth2012、Sph.Rect、Sph.GL 模型的階方差曲線。

圖1 模型球諧系數(shù)階誤差Fig.1 Order error of spherical harmonic models

圖1(a)-(d)顯示，各模型階方差與階誤差曲線均連續(xù)且光滑，所恢復(fù)的球諧系數(shù)模型Sph.Rect、Sph.GL的階方差曲線在各頻段與原模型Earth2012 近乎吻合，可以認(rèn)為球諧分析結(jié)果正確。

比較圖1(a)-(d)可以看出，模型Sph.GL 階誤差在各頻段都小于模型Sph.Rect，模型Sph.GL 比模型Sph.Rect 更接近于模型Earth2012。其中圖1(a)中，盡管各模型零階球諧系數(shù)量級達(dá)到 103的情況下，計算數(shù)據(jù)顯示還原模型Sph.GL 與標(biāo)準(zhǔn)模型零階系數(shù)的階誤差達(dá)到 10-6量級。圖1(a)-(c)顯示模型Sph.GL 在這些頻段的階誤差維持在 10-1以下，而模型Sph.Rect 的階誤差在 10-1以上。值得關(guān)注的是，在圖1(d)中，約2150 階以后，模型Sph.Rect 的階方差和階誤差曲線均出現(xiàn)了急劇下降的現(xiàn)象。經(jīng)過反復(fù)研究比較模型數(shù)據(jù)，發(fā)現(xiàn)直至2159 階次，模型Sph.Rect 球諧系數(shù)的量級與模型Earth2012 保持一致，然而在2160 階的高次項處模型Sph.Rect的系數(shù)量級出現(xiàn)了偏移：、量級為10-27，遠(yuǎn)小于認(rèn)為不再有效，因此按照試驗中矩形離散積分法所求算模型嚴(yán)格意義上只精確到了2159 階。相比之下Sph.GL 模型球諧系數(shù)的量級與Earth2012 模型在全階次保持一致，可以認(rèn)為精確到了2160 階，達(dá)到模型恢復(fù)階數(shù)要求。

為了進一步比較分析兩種方法恢復(fù)模型的精度，給出了球諧系數(shù)的誤差譜，如圖2所示。比較圖2(a)和圖2(b)，在中、低階次部分，模型Sph.GL 球諧系數(shù)誤差小于模型Sph.Rect，兩者在高階次部分差異較小。

圖2 模型Sph.Rect 和模型Sph.GL 球諧系數(shù)誤差譜Fig2 Coefficient error spectral of models

綜合以上分析，模型Sph.GL 的階誤差較小且中、低階次球諧系數(shù)誤差明顯小于模型Sph.Rect，因此相較于矩形離散積分法，Gauss-Legendre 積分法具有更好的球諧系數(shù)模型還原效果。

分別利用模型Earth2012、還原模型Sph.Rect 及Sph.GL 的球諧系數(shù)進行球諧綜合計算，獲得分辨率為5 ′×5 ′的全球數(shù)值模型SHS-Earth2012、SHS-Rect 和SHS-GL。以SHS-Earth2012 為標(biāo)準(zhǔn)，分別記SHS-Rect和 SHS-GL 與其差值為“ DIF-SHS-Rect”、“DIF-SHS-GL”，統(tǒng)計信息示于表1。

表1 球諧綜合結(jié)果差值統(tǒng)計（單位：米）Tab.1 Difference statistics of spherical harmonic synthesis results (Unit:m)

表1顯示，球諧模型GL 恢復(fù)的數(shù)值模型與Earth2012 更加接近，兩者差值絕對值的均值為1.857×10-6m，遠(yuǎn)小于DIF-SHS-Rect 的4.605 m；兩者的標(biāo)準(zhǔn)差差異較大，DIF-SHS-GL 差值大小分布更為集中；DIF-SHS-GL 極大、小值均較小。就所對比的各項指標(biāo)而言，可以認(rèn)為模型SHS-GL 均優(yōu)于模型SHS-Rect，更為接近標(biāo)準(zhǔn)模型。

統(tǒng)計DIF-SHS-Rect、DIF-SHS-GL 兩者在±500 m之間的點數(shù)量，繪制柱狀分布統(tǒng)計圖如圖3所示。

圖3 全球球諧模型誤差分布圖Fig.3 Error distribution of global spherical harmonic models

圖3顯示，DIF-SHS-GL 數(shù)據(jù)更集中地分布在0附近。為了更好地比較兩個模型差值的分布情況，分別統(tǒng)計差值在±1 m、±10 m、±100 m、±500 m 區(qū)間的點數(shù)量和占比情況，如表2所示。

表2 全球球諧模型誤差數(shù)值分布統(tǒng)計Tab.2 Numerical distribution statistics of global spherical harmonic model error

表2顯示，兩組差值在|δh|＜500m范圍內(nèi)的點數(shù)量都遠(yuǎn)超過90%。定義差值在|δh|＜1m的點為高精度還原點，比較兩表，DIF-SHS-GL 的高精度點數(shù)超過DIF-SHS-Rect 5 倍，達(dá)到總點數(shù)的12.160%。同時DIF-SHS-GL 數(shù)值在|δh|＜100m范圍內(nèi)的點數(shù)占比超過了 90%，占到了試驗點的絕大多數(shù)。認(rèn)為|δh|＞ 100m 的點還原效果較差，定義為異常點，其在兩個模型占比分別約為23%、10%。繪制DIF-SHS-Rect、DIF-SHS-GL 兩模型中異常點位置分布圖，示于圖4、圖5。

圖4、圖5中大陸架輪廓清晰可見，兩模型中異常點多數(shù)沿著海陸交界處分布，筆者認(rèn)為可能的原因是陸地與海洋交界處地形變化劇烈，但緯度圈上的格網(wǎng)點是等角距分布，所包含信息不夠豐富，不能細(xì)致地描繪這些區(qū)域；DIF-SHS-GL 中異常點的數(shù)量、分布密度都明顯小于DIF-SHS-Rect；在南北方向高緯度區(qū)域SHS-GL 與SHS-Earth2012 吻合地很好，對比與矩形離散積分法，Gauss-Legendre 積分方法在這些區(qū)域具有獨特優(yōu)勢。

圖4 DIF-SHS-Rect 異常點位分布圖Fig.4 Distribution of anomalous points in DIF-SHS-Rect

圖5 DIF-SHS-GL 中異常點位分布圖Fig.5 Distribution of anomalous points in DIF-SHS-GL

為進一步比較兩種方法所恢復(fù)模型在高緯度區(qū)域的精度，分別研究了模型DIF-SHS-Rect、DIF-SHS-GL在區(qū)域60°N～90°N、60°S～90°S的表現(xiàn)效果，北半球區(qū)域記為DIF-SHS-Rect-N、DIF-SHS-GL-N，南半球區(qū)域記為DIF-SHS-Rect-S、DIF-SHS-GL-S，兩個模型高緯度區(qū)域差值分布如圖6所示。比較圖6(a)與圖6(b)、圖6(c)與圖6(d)，DIF-SHS-GL 在南、北高緯度區(qū)域的值更加集中于 0m 附近，明顯優(yōu)于DIF-SHS-Rect。

圖6 高緯度區(qū)域差值分布圖Fig.6 Difference distribution diagram at high latitudes

表3為兩組球諧綜合結(jié)果在高緯度的差值統(tǒng)計情況，結(jié)合兩種方法的原理，由于Gauss-Legendre 積分法賦予采樣點Gauss-Legendre 權(quán)重，從而很好地改善了積分法在高緯度區(qū)域表達(dá)精度問題。如表3所示，DIF-SHS-GL 在南、北高緯度區(qū)域的標(biāo)準(zhǔn)差約為DIF-SHS-Rect 的1/6。對比表1和表3，還可以發(fā)現(xiàn)研究區(qū)域DIF-SHS-GL-N、DIF-SHS-GL-S 的標(biāo)準(zhǔn)差遠(yuǎn)小于全球區(qū)域的標(biāo)準(zhǔn)差20.920 m，約為其1/2，表達(dá)精度優(yōu)于全球平均標(biāo)準(zhǔn)；而 DIF-SHS-Rect-N、DIF-SHS-Rect-S 的標(biāo)準(zhǔn)差比表1中全球區(qū)域標(biāo)準(zhǔn)差大，低于全球平均標(biāo)準(zhǔn)。

表3 高緯度區(qū)域球諧綜合結(jié)果差值統(tǒng)計（單位：米）Tab.3 Difference statistics of spherical harmonic synthesis results at high latitudes (Unit:m)

表4為球諧模型誤差在南、北高緯度區(qū)域數(shù)值分布統(tǒng)計信息，比較發(fā)現(xiàn)DIF-SHS-GL 在各數(shù)值區(qū)間的數(shù)量占比均明顯大于DIF-SHS-Rect。對比表4中高緯度誤差與表2中全球誤差分布，各模型在高緯度區(qū)域數(shù)量占比均有提高，其中：DIF-SHS-Rect 和DIF-SHS-GL 在|δh|＜1m區(qū)間數(shù)量占比分別平均提高1.678%和16.444%；在|δh|＜10m 區(qū)間數(shù)量占比分別平均提高12.041%和30.474%；在|δh|＜100m 區(qū)間數(shù)量占比分別平均提高 13.037%和 8.596%。除|δh|＜100m 區(qū)間外，DIF-SHS-GL 提高更顯著，但DIF-SHS-GL 在南、北高緯度區(qū)域|δh|＜100m 區(qū)間的差值數(shù)量占比均超過98%，遠(yuǎn)大于DIF-SHS-Rect。分析球諧模型表達(dá)高緯度與全球精度存在差異原因為：球面采樣中高緯度采樣點更密集，富含更多地形特征信息，因而恢復(fù)模型描繪高緯度區(qū)域能力更強，在各區(qū)間數(shù)量占比也會高于全球平均水平。

表4 高緯度區(qū)域差值模型數(shù)值分布統(tǒng)計Tab.4 Numerical distribution statistics of difference models at high latitudes

3 結(jié) 論

本文首先建立矩形等角格網(wǎng)和Gauss-Legendre 格網(wǎng)；然后依據(jù)Earth2012 模型球諧綜合計算離散格網(wǎng)點地形數(shù)據(jù)，并分別用矩形離散積分法、Gauss-Legendre 積分法構(gòu)建全球2160 階地形球諧系數(shù)模型Sph.Rect、Sph.GL，并還原出全球5′×5 ′地形數(shù)值模型SHS-GL、SHS-Rect；最后討論了兩種方法構(gòu)建地形球諧模型的精度。可得到以下結(jié)論：

（1）比對 Sph.GL 模型和 Sph.Rect 模型與Earth2012 模型的階誤差表明，在輸入數(shù)據(jù)近似相等的情況下，Gauss-Legendre 積分法還原效果更好，與標(biāo)準(zhǔn)模型更接近；

（2）比對數(shù)值模型SHS-GL 和SHS-Rect 表明，在海陸交界等地形變化劇烈區(qū)域，矩形離散積分法和Gauss-Legendre 積分法描繪地形細(xì)節(jié)能力均顯不足；

（3）相比矩形離散積分法，Gauss-Legendre 積分方法表達(dá)高緯度區(qū)域地形更具優(yōu)勢，且該方法描繪高緯度地區(qū)地形精度高于描繪全球的平均水平，能較好解決兩極附近存在采樣點較密、數(shù)據(jù)不均勻的問題；

（4）通過數(shù)組拓展升維等方法可以有效解決并行中內(nèi)存讀寫沖突問題，合理的并行方案及快速傅里葉變換技術(shù)可以充分利用計算機計算資源，進而大大提高調(diào)和分析計算中耗時長、效率低的問題?？蔀楦唠A地形球諧模型構(gòu)建方法選擇提供參考。