單建晨,李姍姍,范 雕,李新星,黃 炎,邢志斌
(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ù)計算策略,大大提高了計算效率。
根據(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。
由于離散勒讓德函數(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ù)。
本文選取科廷大學(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)了快速運算。
根據(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ū)間使用割線迭代法計算零點值。
綜合上述式(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,較好地解決了串行編譯耗時多、效率低的問題。
按照上文原理算法,建立矩形等角格網(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途取?/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
本文首先建立矩形等角格網(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)建方法選擇提供參考。