基于Laplace方程的球面正交曲線格網生成方法

盧 付 強,宋 志 堯,雷 夢 玲,丁 凱 孟

(1.南京師范大學虛擬地理環境教育部重點實驗室，江蘇 南京 210023;2.南京師范大學大規模復雜系統數值模擬江蘇省重點實驗室，江蘇 南京 210023;3.江蘇省地理信息資源開發與利用協同創新中心，江蘇 南京 210023)

基于Laplace方程的球面正交曲線格網生成方法

盧 付 強1,2,3,宋 志 堯1,2,3,雷 夢 玲1,3,丁 凱 孟1,3

(1.南京師范大學虛擬地理環境教育部重點實驗室，江蘇 南京 210023;2.南京師范大學大規模復雜系統數值模擬江蘇省重點實驗室，江蘇 南京 210023;3.江蘇省地理信息資源開發與利用協同創新中心，江蘇 南京 210023)

為克服基于傳統經緯度格網的全球海洋及大氣數值模式產生的極點問題，同時顧及全球格網編碼便捷及信息查詢,該文將立方體球心投影的格網經球極投影轉換到二維平面區域，再應用Laplace方程數值求解生成正交曲線格網；然后采用球極逆投影,得到球面上一種新的正交曲線格網。通過與等角球心投影、保角變換兩種球面格網生成方法的結果相比較,從格網兩邊的長度、單元的面積、正交性等指標進行統計分析，對所生成格網的質量進行評價。結果表明,等角球心投影法生成的格網單元大小較均勻，但正交性不理想；Laplace方程法生成的正交曲線格網稍優于保角變換法，可作為應用有限差分方法和有限體積方法的全球海洋及大氣大尺度數值模擬的基礎格網。

Laplace方程；球極投影；正交曲線格網；Neumann正交化

0 引言

隨著全球氣候逐漸變暖，極端天氣現象發生的頻率呈升高趨勢,需對全球氣候變化進行及時預測，以便采取可能的措施減小惡劣天氣產生的危害。而這些預測在很大程度上依賴于全球大氣模式、海洋模式的發展，因此需要構建一個理想的全球大氣、海洋數值模型都適用的網格，以便進一步實現相互耦合的數值模擬。建立全球大氣、海洋數值模型的第一步就需要選擇一個好的全球離散網格。目前在大氣、海洋數值模式中所采用的網格分為兩大類[1]：一類是有結構網格，網格的排列結構縱橫有序，相鄰節點或網格單元的數目是固定的，易于結點(單元)的編號和布置變量，便于方程的離散和計算編程，如傳統或退化的經緯度格網[2,3]、立方體球心投影網格[4,5]、陰陽交疊網格[6,7]等。另一類是無結構網格，即由任意三角形或六邊形組成的網格，如基于正二十面體產生的球面三角網格[8,9]、球面正六邊形網格[10,11]。此類網格在全球大氣或海洋運動數值模型中得到了廣泛應用，其主要優點是：整體上網格大小較均勻。但在這類網格上對數據進行操作，需要預先對每個網格建立索引、拓撲關系表，并且在這類網格上建立數值模型的求解方法通常較復雜。

正交曲線網格屬于有結構網格，具有很多優點：首先，網格的編碼與領域查詢比較方便,模型方程組在正交曲線坐標系中的表達形式比較簡潔；其次，便于數值有限差分方法或者有限體積方法的設計和實現等。許多大氣數值模式建立在經緯度網格上，但是經緯度網格在趨于兩個極點時經線收斂，除了譜變換方法，在采用基于格點的方法進行數值計算時會引發極點問題[12]，同時也產生了解決這些問題的方法：高緯度濾波、半隱式半拉格朗日方法等。但是這些方法在全球海洋模式中卻難以應用[13]，所以必須考慮采用一種新的正交網格，采用數值求解Laplace方程組，結合球極投影方法，可生成球面正交曲線網格，既利于全球地理信息的表達，也適用于地球流體大尺度運動的數值模擬，同時也是生成三維球體準正交的曲面六面體網格的基礎。

1 球面四邊形網格生成方法

1.1 立方體球心投影方法

考慮球的內接立方體，分別在立方體的每個表面上建立局部二維笛卡爾坐標，然后做一個伸縮變換，將表面上每個點變換到外接球面上。假設球半徑為R，可知立方體的邊長2a與半徑有關系式：

(1)

選取立方體中心為三維笛卡爾坐標系的原點，坐標軸X、Y、Z分別與立方體的前面、右面、頂面垂直(圖1)，其中左子圖為立方體各面的展開圖，右子圖為立方體與外接球面在三維笛卡爾坐標系中的位置關系。記(x,y)為立方體每個表面上的局部二維笛卡爾坐標，有-a≤x,y≤a。

圖1 立方體6個表面的展開圖及其外接球面在三維笛卡爾坐標系中的位置

Fig.1 Layout of an open cube and the cube and its circum-sphere under 3D absolute Cartesian coordinates

從立方體表面變換到球面上有兩種方法：

(1)立方體球心等距投影網格。直接對局部坐標進行等距網格剖分(圖2)，其中顯示了在立方體表面X=a上的局部二維平面笛卡爾坐標(x,y)與半徑為R的球面上點(X,Y,Z)之間的關系：

(2)

(2)立方體球心等角投影網格。首先對局部笛卡爾坐標做一個變換：

(3)

圖2 立方體表面X=a上的點經過球心投影到球面上

Fig.2 Gnomonic projection of point on the cube surfaceX=ato the sphere surface

1.2 保角變換方法

應用保角變換將網格點在各復平面與球面之間變換時，網格的角度保持不變，并且除去8個角點，映射在每個網格點上是光滑的，保角變換方法的主要思想是先采用廣義球極投影將球面上的點映射到復平面上，然后在復平面之間采用Taylor 級數展開式實現一系列復變換及逆變換[14]。保角變換生成的網格總體上正交性很好并且網格的長寬比接近1，即形狀與正方形相近。但在8個奇點(角點)處，網格偏離正交性固定在30°，在角點附近的網格正交性也受到影響，正交性有所降低。盡管不如球面經緯網格在極點處收斂那么嚴重,嚴格的正交性要求將導致網格在8個奇點處收斂,可以對網格在局部進行拉伸以改善這種狀況。

2 球面正交曲線網格生成

采用某種投影方法將球面上的網格變換到二維平面區域，然后借鑒平面上數值網格生成方法生成正交曲線網格，最后采用逆變換，將生成的網格變換到球面上。這個過程的主要思路可由圖3描述。在一般情形下，采用投影變換后，網格的幾何形狀、面積、角度等屬性都會發生變化。在此選取常用的3種地圖投影變換方法：圓錐、圓柱及赤平投影[18,19]，這些變換是保持正交性質的，即網格在球面與平面之間變換時，格網的角度不發生變化。

圖3 算法的總體路線

Fig.3 The process of the presented algorithm

2.1 地圖投影變換

由于立方體6個表面上的網格是對稱的，所以只需在一個面上生成網格，然后采取合適的旋轉變換就可得到覆蓋整個球面的網格。在此選取赤道區域所在的一個球面區域作為基準面，圖4顯示了立方體球面前表面上初始網格分別經過Lambert正形圓錐投影、Mercator投影、赤平投影變換后的網格。立方體球面前表面上的格網關于格網中心具有四重對稱性，從圖4可看出：經過圓錐或圓柱投影，格網的四重對稱性有損失；而經過赤平投影變換，格網的對稱性得到保持，所以采用赤平投影或球極投影。

圖4 經過圓柱、圓錐以及赤平投影變換之后,立方體球面前表面上格網點的分布

Fig.4 Distribution of gird points on front face of cubed sphere under Mercator projection,Lambert conformal conic projection and stereographic projection

2.2 Laplace微分方程方法

Laplace微分方程是一種橢圓形偏微分方程，橢圓形網格生成是通過求解擬線性橢圓形偏微分方程組的一種網格生成方法，可以對網格的光滑性及正交性進行控制。在生成球面網格時，初始網格取立方體等距球心投影網格，將在球面上物理網格的每個表面通過球極投影變換到二維的平面空間，然后對網格施加所需的限制條件，在投影平面上求解橢圓型偏微分方程組。在大部分基于有限差分方法的地球流體動力學數值模擬中，網格的正交性及網格尺寸的均勻性是最重要的，對網格施加正交性有Neumann條件限制和Dirichlet條件限制[20]兩種方法。Neumann條件限制：網格點可以沿著每個面的邊界滑動，而Dirichlet條件限制則是固定在邊界上的點，通過調整網格內點并且保持網格尺寸的準一致性。Neumann條件限制通過要求未知量對兩個參數的偏導數的乘積等于零來計算每個網格點沿網格邊界的偏移量，但這種做法將導致網格單元在角點鄰近區域趨于收斂，網格單元的面積偏小，可在局部對網格進行拉伸來改善這種狀況。

假定一個二維初始變換W:(ξ,η)→(x,y),將計算空間[1,M]×[1,N]變換到物理空間Ω?R2,其中正整數M、N表示沿兩個方向的網格點數目，網格線索引為ξ=1,2,…,M;η=1,2,…,N。初始變換常由一般的方法如代數方法中的超限插值技術得到。對于給定的初始變換，二維平面正交曲線網格的生成方法多采用數值求解下列擬線性偏微分方程組[21]:

g22(xξξ+Pxξ)-2g12xξη+g11(xηη+Qxη)=0
g22(yξξ+Pyξ)-2g12yξη+g11(yηη+Qyη)=0

(4)

其中：

(5)

(6)

(7)

其中:(x,y)為原二維物理區域中的網格點坐標，(ξ,η)為變換區域中網格點坐標。網格的度量項g11、g22分別表示網格沿兩個方向的長度，g12可用于表示網格單元的角度大小。P、Q表示網格點疏密的控制函數，可通過改變P、Q對網格點的分布進行調整，若P=Q=0，則產生步長均勻的網格。

為了求解擬線性橢圓形方程組(4)，需給出必要的邊界條件，此處設置控制函數P=Q=0，取Neumann正交邊界條件，允許邊界點通過自由滑動達到邊界上的正交。采用有限差分方法求解方程組(4)，其中網格的度量項g11、g12、g22可以根據網格點坐標表達式預先計算出來，或者采用差分近似方法，系數中含有待求解函數的非線性項，采用前一次的迭代值計算這些系數，只對二階偏導數以及與控制函數有關的導數項進行離散，采用中心差分近似，選取映射到計算空間中的離散步長Δξ=Δη=1，則:

(8)

(9)

采用超松弛迭代方法(SOR)求解上述方程組：

(10)

其中wi,j為松弛因子，應滿足條件：0

2.3 Neumann正交邊界條件

在求解區域邊界上應該滿足正交性條件：

XξXη=0,在ξ=1,M

(11)

XξXη=0在η=1,N

(12)

其中X表示網格點的二維坐標，采用有限差分方法對式(11)、式(12)進行近似求解，在迭代求解方程組時，每次循環中邊界條件式都要進行更新。為了防止邊界點滑動距離過大，超過其鄰近點，需要設置臨界值，特別是在鈍角區域附近容易出現這種情形。合適的臨界值對邊界每個點的滑動距離起限制作用，一般要求在每次迭代過程中距離邊界上最近兩個鄰近點的距離至多減少一半。

3 結果分析與比較

3.1 網格質量評價指標

網格的4個指標有：網格分別沿參數ξ方向與η方向的長度hξ,hη網格單元的面積大小area以及角度的正交性偏差DO：

(13)

(14)

3.2 計算結果分析

選取立方體的頂面Z=a作為基礎，在其上建立局部二維笛卡爾坐標系(x,y)，然后變換到半徑為R的球面上。取a=1，采用北極投影，將立方體球心等距投影網格映射到與北極點相切的二維平面上，數值求解橢圓形微分方程組(4)，采用Neumann正交邊界條件(式(11)、式(12))，控制函數取值P=Q=0；然后再采用北極逆投影得到球面上的正交網格。圖5給出了在球極投影下，各網格點處地圖放大系數的等值線圖，可以看出投影系數接近1 ，即網格經過球極投影后，幾何度量變化不是很大。

目前，也可采用等角或等距球心投影方法以及保角映射方法生成立方體球面網格，但等角或等距球心投影產生的網格不具有正交性，而保角變換方法或求解橢圓形方程產生的網格正交性較好。下面給出上述方法生成網格的度量性質的比較結果。

選取網格單元數分別為45×45、90×90、180×180的情形時進行計算，相對應的空間分辨率大約為2°、1°、30°。表1給出了不同網格數時，分別采用等角投影方法、保角變換、Laplace方程方法生成網格的統計結果，包括網格單元最大邊與最短邊的長度、單元角度的正交性最大偏差以及平均偏差、網格單元面積的最大值和最小值。

圖5 各網格點處的地圖放大因子等值線

Fig.5 Map factor of each grid point after elliptic-smoothing under north polar stereographic projection

如果頂面的網格已經生成，則采用合適的旋轉變換生成其他5個表面上的正交網格。圖6顯示了網格數在45×45×6時，采用等角投影方法、保角變換、 Laplace方程方法生成球面網格的正交性偏差及網格單元面積的變化。可看出保角變換、Laplace方程方法生成的球面網格非常相似，絕大多數網格的正交性偏差不超過10°，只是接近8個角點的網格單元的正交性偏差增大了一些，最大偏差發生在8個角點處。另外，從表1也可看出，保角變換與Laplace方程方法生成網格的統計量差別相當小，基本在同一個數量級。綜合起來，在相同情形下，等角投影方法生成球面網格單元面積的大小較均勻，但正交性偏差較大；而保角變換、Laplace方程方法生成球面網格單元的面積總體上不如等角投影方法的均勻，尤其是在8個角點附近區域網格單元的面積偏小。

圖 7 給出了網格逐次分半加密時，分別采用三種方法，Fortran編程的CPU消耗時間隨網格數目的變化，可以看出在同等格網數目的條件下，求解Laplace方程組耗時最多，等角球心投影方法耗時最少，保角變換方法耗時稍多于等角球心投影方法。即便如此，生成空間分辨率為幾分的格網最多時間耗費不過10 s，并且在發展全球大氣或海洋數值模式過程中，網格生成是一次性的，這個時間是可接受的。

4 結語

基于立方體球心投影網格，結合球極投影技術，采用Laplace方程組在球面上生成了一種廣義正交曲線網格，并且對網格的質量進行分析評價。通過與等角球心投影、保角變換方法的結果進行比較，表明等角球心投影生成的網格單元大小較均勻，但網格不具有正交性；而保角變換方法、數值求解微分方程組生成的網格的正交性質較好，但在8個角點鄰域的網格單元面積偏小。在應用于實際大氣或海洋大尺度運動的具體數值模型中，這個缺陷可以通過對網格進行拉伸，或采用隱式時間積分方法來彌補。因此,這種廣義正交曲線網格可用于地理空間數據的表達以及全球大氣或海洋運動的數值模擬。

表1 在不同網格數情形下三種方法生成網格的統計結果

Table 1 Quality statistics of grids generated by three different methods with different cell-numbers

45×4590×90180×180等角投影方法保角變換方法Laplace方程方法網格邊長單元面積正交性偏差(°)網格邊長單元面積正交性偏差(°)網格邊長單元面積正交性偏差(°)最大值3.4906E-021.7453E-028.7266E-03最小值2.4684E-021.2342E-026.1707E-03均值3.2367E-021.6212E-028.1132E-03方差6.9822E-061.6732E-064.0929E-07最大值1.2000E-033.0648E-047.7482E-05最小值8.7658E-042.1682E-045.3487E-05均值1.0343E-032.5857E-046.4642E-05方差8.5667E-095.364E-0103.3550E-011最大值28.814229.421529.7109均值8.06318.06568.0661最大值5.1088E-022.5973E-021.3097E-02最小值1.7561E-027.1262E-032.8604E-03均值3.2365E-021.6188E-028.1008E-03方差3.2343E-058.1043E-062.0272E-06最大值1.3519E-033.4438E-048.7732E-05最小值4.2396E-046.9325E-051.1316E-05均值1.0343E-032.5856E-046.4643E-05方差2.2336E-081.3998E-098.759E-011最大值12.922512.930312.9310均值0.45220.22610.1131最大值3.6401E-021.8622E-021.0004E-02最小值6.5321E-032.5311E-038.9817E-04均值3.1804E-021.5916E-027.9523E-03方差1.7828E-054.7045E-061.5216E-06最大值1.3000E-033.4978E-041.0225E-04最小值7.5305E-051.1096E-051.8763E-06均值1.0343E-032.5857E-046.4643E-05方差5.2956E-083.7047E-093.137E-010最大值12.688912.564712.5055均值0.73910.35000.2063

圖6 不同投影生成網格的角度正交性偏差與網格單元面積大小 圖7 三種方法的計算效率隨網格數目的變化

Fig.6 Angle deviation orthogonality and cell area on equi-angular sphere grids based on different method Fig.7 The efficiency of three methods under different grid-numbers

[1] STANIFORTH A,THUBURN J.Horizontal grids for global weather and climate prediction models:A review[J].Quarterly Journal of the Royal Meteorological Society,2012,138:1-26.

[2] HOLLOWAY L,SPELMAN M J,MANABE J.Latitude-longitude grid suitable for numerical time integration of a global atmospheric model[J].Monthly Weather Review,1973,101(1):69-78.

[3] FADEEV R Y.Algorithm for reduced grid generation on a sphere for a global finite difference atmospheric model[J].Computational Mathematics and Mathematical Physics,2013,53(2):237-252.

[4] RONCHI C,IACONO R,PAOLUCCI P S.The "Cubed Sphere":A new method for the solution of partial differential equations in spherical geometry[J].Journal of Computational Physics,1996,124(1):93-114.

[5] MCGREGOR J L.Semi-Lagrangian advection on a cubic gnomonic projection of the sphere[J].Atmosphere-Ocean,1997,35(1):153-169.

[6] KAGEYAMA A,SATO T."Yin-Yang grid":An overset grid in spherical geometry[J].Geochem.Geophys.Geosyst.,2004,5,Q09005,doi:10.1029/2004GC000734.

[7] 艾細根,劉宇迪,崔新東,等.基于球面陰陽網格數值計算研究進展[J].氣象科技,2014,42(1):13-22.

[8] TOMITA H,SATOH M,GOTO K.An optimization of the Icosahedral Grid modified by Spring Dynamics[J].Journal of Computational Physics,2002,183(1):307-331.

[9] IGA S,TOMITA H.Improved smoothness and homogeneity of icosahedral grids using the spring dynamics method[J].Journal of Computational Physics,2014,258(1):208-226.

[10] MIURA H,KIMOTO M.A comparison of grid quality of optimized spherical hexagonal-pentagonal geodesic grids[J].Monthly Weather Review,2005,133(10):2817-2833.

[11] 童曉沖,奔進,汪瀅.利用數值投影變換構建全球六邊形離散格網[J].測繪學報,2013,42(2):268-276.

[12] NAUGHTON M,COURTIER P. A pole problem in the reduced Gaussian grid[J].Quarterly Journal of the Royal Meteorological Society,1994,120:1389-1407.

[13] ADCROFT A,CAMPIN J M,HILL C,et al.Implementation of an atmosphere-ocean general circulation model on the expanded spherical cube[J].Monthly Weather Review,2004,132:2845-2863.

[15] NAIR R D,THOMAS S J,LOFT R D.A discontinuous Galerkin transport scheme on the cubed sphere[J].Monthly Weather Review,2004,133(4):814-828.

[16] LU F,SONG Z,FANG X.Generation of Orthogonal Curvilinear Grids on the Sphere Surface Based on Laplace-Beltrami Equations[C].IOP Conference Series:Earth and Environmental Science,2014,(19) 012012.

[17] PUTMAN W M,LIN S J.Finite-volume transport on various cubed-sphere grids[J].Journal of Computational Physics,2007,227(1):55-78.

[18] 沈桐立,田永祥,葛孝貞.數值天氣預報[M].北京:氣象出版社,2003.60-71.

[19] 王美玲,付夢印.地圖投影與坐標變換[M].北京:電子工業出版社,2014.51-125.

[20] THOMPSON J F,SONI B K,WEATHERILL N P.Handbook of Grid Generation[M].CRC Press,1999,ISBN 0-8493-2687-7.

[21] SPEKREIJSE S P. Elliptic grid generation based on Laplace equations and algebraic transformations[J].Journal of Computational Physics,1995,118(1):38-61.

Generation Method of Orthogonal Curvilinear Grids on Sphere Surface Based on Laplace Equations

LU Fu-qiang1,2,3,SONG Zhi-yao1,2,3,LEI Meng-ling1,3,DING Kai-meng1,3

(1.Key Laboratory of Virtual Geographic Environment(Ministry of Education),Nanjing Normal University，Nanjing 210023;2.Jiangsu Key Laboratory for Numerical Simulation of Large Scale Complex Systems,Nanjing Normal University,Nanjing 210023;3.Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application,Nanjing 210023,China)

In order to overcome the pole-problem induced by grid-based methods on traditional longitude-latitude grid for global atmosphere and ocean circulation models,a new kind of orthogonal curvilinear grids based on the gnomonic cube grids on the sphere was presented.Firstly,the gnomonic cube grids was on a two-dimensional plane by stereographic projection,then Laplace equations are solved to generate quasi-orthogonal curvilinear grids on the projection plane.Finally,the grids are transformed back to the sphere surface by the inverse of stereographic projection.The produced grid configuration has quasi-uniform cell-size on the whole sphere surface except the neighbourhood of eight major corner vertices.Moreover,another two methods which include equi-angular gnomonic projection and conformal mapping were also introduced to generate spherical grids.The resulting comparisons between three methods would be presented at last.Some quantities such as the grid length along two directions,the angle deviation from orthogonality and the area of the cell were used to evaluate the quality of the grids generated by three methods.The results show that the equi-angular gnomonic projection produces grids of the most uniform size,but are non-orthogonal,while elliptic-smoothing combined with stereographic projection method produces quasi-orthogonal curvilinear grids that are much like that of conformal mapping method.The small grid cells in the neighbourhood of eight corner vertices could be improved by stretching method,or this small defect could be circumvented by implicit time integration method,which demonstrate the kind of grid newly produced is fit to be a model grid on which the finite difference method or finite volume method can be implemented for numerical simulation of global atmosphere or ocean dynamics on large scale.

Laplace equations;stereographic projection;orthogonal curvilinear grids;Neumann orthogonal boundary condition

2015-06-18；

2015-10-13

國家自然科學基金資助項目(41306078)

盧付強(1984-),男,博士,研究方向為地圖學與地理信息系統。E-mail:lufuqiang2000@163.com

10.3969/j.issn.1672-0504.2016.01.004

P208

A

1672-0504(2016)01-0018-06