吳 竑，王小鴿，楊廣文，朱昌磊，宋順強

清華大學(xué) 計算機科學(xué)與技術(shù)系，北京 100084

THREP：面向地球系統(tǒng)模式的插值算法研究平臺

吳 竑，王小鴿，楊廣文，朱昌磊，宋順強

清華大學(xué) 計算機科學(xué)與技術(shù)系，北京 100084

1 引言

1.1 背景

插值算法[1]是一類應(yīng)用十分廣泛的算法，它們在圖形圖像處理、數(shù)值計算等許多科學(xué)研究領(lǐng)域都是不可或缺的重要數(shù)學(xué)手段。插值算法在耦合地球系統(tǒng)模式中也有著極其重要的應(yīng)用。

耦合地球系統(tǒng)模式[2]是一套對地球系統(tǒng)進行定量研究的軟件工具，地球科學(xué)相關(guān)領(lǐng)域的科學(xué)家根據(jù)專業(yè)知識對地球各個子系統(tǒng)進行建模，編寫各自領(lǐng)域的模擬程序，再由這些模擬程序相互耦合，成為一個完整的地球系統(tǒng)模式。為了使不同研究人員的代碼能夠耦合起來形成一個整體，往往需要他們共同遵循某個耦合軟件框架[3]。目前較為流行的地球系統(tǒng)模式耦合框架有歐盟PRISM所采用的OASIS[4]和美國CCSM3所采用的COUPLER[5]系列等。然而，由于地學(xué)領(lǐng)域的科學(xué)家們對地球各子系統(tǒng)進行建模時所使用的數(shù)值方法不同，導(dǎo)致他們在編程時使用不同的地球網(wǎng)格劃分方式（如經(jīng)緯網(wǎng)格、高斯網(wǎng)格[6]），從而對地球系統(tǒng)模式中不同分量模式間的數(shù)據(jù)交換提出了挑戰(zhàn)。因此，在不同網(wǎng)格間進行數(shù)據(jù)交換時需要一種數(shù)據(jù)映射方式——插值。此時，可以將每個網(wǎng)格單元看做一個插值點，并用網(wǎng)格單元中心點坐標表示插值點坐標。

1.2 相關(guān)工作

SCRIP[7]是一個面向地球系統(tǒng)模式的插值算法庫，可以離線地生成地球系統(tǒng)模式研究人員所需地球網(wǎng)格間的插值權(quán)重系數(shù)文件，它被使用于CPL6等耦合器中。但是SCRIP中的插值算法不很完善，其守恒插值算法只支持一階，雙三次插值算法也沒有實現(xiàn)求導(dǎo)數(shù)的方案。此外，由于SCRIP中的數(shù)據(jù)結(jié)構(gòu)過于簡單，它能夠支持的網(wǎng)格類型非常有限，這些因素都制約了其在地球系統(tǒng)模式中的作用。

ΤHU-Remap（https：//github.com/xunzhang/Τ-Remap）是針對SCRIP的缺陷，開發(fā)的一套全新插值算法庫。它一方面重構(gòu)了SCRIP中一些算法，如反距離權(quán)重算法、雙線性算法，大大提高了它們的計算性能；另一方面改進了SCRIP中的一些不足，如實現(xiàn)了雙三次差值中的求導(dǎo)、提供守恒插值中閾值的自動調(diào)整等。但是，ΤHU-Remap對非規(guī)則網(wǎng)格的支持依然不夠，這不利于分量模式與耦合器的長期發(fā)展。

ESMF（Earth System Modeling Framework）[8]項目由NASA地球科學(xué)辦公室于2000年9月提出，是一套面向地球系統(tǒng)模式的并行軟件開發(fā)框架。它不含任何分量模式代碼，而是提供標準框架接口，用戶可以基于此接口實現(xiàn)天氣、氣候、數(shù)據(jù)同化等相關(guān)領(lǐng)域的應(yīng)用模型并將它們耦合起來成為一個完整的地球系統(tǒng)模式。ESMF模塊化設(shè)計使其具有極高的可重用性、互操作性和可移植性，這將有利于模式專家專注于模式本身的發(fā)展并減少其耦合開銷。ESMF主要提供Superstructure和Infrastructure兩部分功能供開發(fā)者使用：Superstructure中封裝了地球系統(tǒng)模式中的耦合驅(qū)動和接口，也模塊化了分量模式中物理過程、動力框架及一些子過程；Infrastructure則提供更加底層的抽象，它又可細分為數(shù)據(jù)處理和公用模塊兩部分，數(shù)據(jù)處理部分主要封裝所有與地球系統(tǒng)模式中數(shù)據(jù)處理相關(guān)的數(shù)據(jù)結(jié)構(gòu)及功能，如Array、Grid、Field等，公用模塊部分則提供了如數(shù)值計算、時間管理、進程管理、I/O、容錯日志處理的基本功能。

ESMF-Regridding[9]是ESMF提供的一個基于ESMF框架開發(fā)的插值應(yīng)用，它能并行生成離線插值權(quán)重系數(shù)文件，特別是它能支持非規(guī)則網(wǎng)格間的插值。但是，ESMF-Regridding僅以一個小工具的形式出現(xiàn)在ESMF中，沒有抽象插值在應(yīng)用層面的功能供開發(fā)人員使用，因此不夠靈活和通用。插值算法開發(fā)者往往需要閱讀大量ESMF源碼，作出相應(yīng)修改并實現(xiàn)各自新的插值算法，難度極大，違背了ESMF的設(shè)計目標。此外，ESMF-Regridding沒有提供驗證插值算法正確性與性能的功能，不利于深入研究與對比各插值算法的特性。

雖然目前已有SCRIP、ΤHU-Remap、ESMF-Regridding等面向地球系統(tǒng)模式的插值軟件，但隨著地球系統(tǒng)模式研究的深入，對插值算法提出了新要求。為插值算法的研究、測試和評價建立一個統(tǒng)一的平臺，具有重要的實踐意義。基于ESMF框架，開發(fā)了一個通用、靈活、并行的插值研究平臺ΤHREP（Τsinghua Regridding Platform），它能同時滿足插值算法研究人員、耦合器開發(fā)者、模式專家三類用戶需求。ΤHREP集成了ΤHU-Remap與ESMF-Regridding中的插值算法，封裝了易用的插值接口，并提供了驗證、對比等功能。

2 平臺設(shè)計

ΤHREP綜合考慮了算法研究人員、耦合器開發(fā)者及模式專家三類人對插值的需求，是一個通用、靈活、并行的面向地球系統(tǒng)模式的插值平臺。它抽象了插值過程，封裝了插值流程，便于已有插值算法庫的接入，利于長期發(fā)展和標準化。同時，為滿足不同用戶，ΤHREP在保證易用性的同時，從底到上提供了多層接口。ΤHREP還提供了插值檢驗功能，便于深入研究各算法的精度、性能等相關(guān)問題。

2.1 平臺架構(gòu)

ΤHREP插值平臺的總體結(jié)構(gòu)圖見圖1，主要包括輸入文件、配置模塊、權(quán)重生成模塊、插值模塊和檢驗?zāi)K幾個部分。輸入文件含網(wǎng)格數(shù)據(jù)文件和物理量數(shù)據(jù)文件，配置模塊負責(zé)配置一次插值過程需要的信息，權(quán)重生成模塊包含了初始化和搜索兩個子模塊，ΤHREP通過初始化模塊建立相應(yīng)數(shù)據(jù)結(jié)構(gòu)的對象，在此基礎(chǔ)上搜索找出各算法所需格點信息，生成插值權(quán)重，插值模塊通過矩陣向量乘計算出實際物理量插值后的結(jié)果從而完成整個插值過程，檢驗?zāi)K將對插值的精度和性能作進一步分析并畫出圖形。

圖1 ΤHREP平臺架構(gòu)

2.2 通用性

目前地球系統(tǒng)模式中常用的插值算法有反距離權(quán)重插值[7]、雙線性插值、雙三次插值[10]、Patch插值[11-12]、三次樣條插值[13]、守恒插值[7，14]等，ΤHREP里將所有插值算法抽象成了如下兩個步驟：

（1）計算插值權(quán)重系數(shù)矩陣。這個階段又可細分為兩個子過程：首先根據(jù)源網(wǎng)格和目標網(wǎng)格間的拓撲關(guān)系，對每個目標網(wǎng)格單元（cell）找到其對應(yīng)位置在源網(wǎng)格的“若干”網(wǎng)格單元（這里的“若干”由不同插值算法決定，如雙線性插值要找出包住目標網(wǎng)格單元的4個源網(wǎng)格單元，即目標網(wǎng)格在找出的4個源網(wǎng)格單元連線構(gòu)成的四邊形內(nèi)部）。然后根據(jù)插值算法計算出每個找到的源網(wǎng)格單元相對于該目標網(wǎng)格單元的權(quán)重。因為每個目標網(wǎng)格單元都要計算其對應(yīng)的“若干”個權(quán)重值，因此這一步完成后，會得到一個權(quán)重矩陣W，W的每行對應(yīng)每個目標網(wǎng)格單元id，每列對應(yīng)每個源網(wǎng)格單元id。這里W是一個稀疏矩陣，如雙線性插值的權(quán)重矩陣每行只有4個非零元素。

ΤHREP中使用NetCDF[15]標準的文件格式作為整個平臺數(shù)據(jù)流中輸入、輸出文件的格式。NetCDF文件是一種格式化的二進制文件，它既有二進制文件的數(shù)據(jù)連續(xù)性，又能使用戶像文本文件那樣方便地讀取相應(yīng)屬性的數(shù)據(jù)。采用NetCDF文件的另一個重要原因是它已經(jīng)是研究地球系統(tǒng)模式人員的約定標準，分量模式中使用的網(wǎng)格文件均采用這種格式，因此用NetCDF作為插值結(jié)果的輸出文件能更便于用戶直接使用。

考慮通用性的另一方面是功能可擴展，而數(shù)據(jù)結(jié)構(gòu)直接決定了能否支持眾多功能。因此，通用的數(shù)據(jù)結(jié)構(gòu)也是系統(tǒng)通用性的關(guān)鍵。整個插值過程中主要的數(shù)據(jù)包括網(wǎng)格信息和物理量信息，插值所用的數(shù)據(jù)結(jié)構(gòu)主要就是為了便于插值計算而在這兩部分數(shù)據(jù)信息基礎(chǔ)上建立的。在ΤHREP中盡量保證了數(shù)據(jù)的原始邏輯，并刻畫了網(wǎng)格單元間的相鄰邏輯關(guān)系，這樣可以保證功能可擴展，便于平臺的進一步研發(fā)，也為以后我們或者其他開發(fā)者封裝更高級的數(shù)據(jù)結(jié)構(gòu)提供了可能。

2.3 靈活性

在眾多地球系統(tǒng)模式的研發(fā)和使用人員中主要有三類人員涉及插值。第一類是模式專家，他們負責(zé)開發(fā)分量模式，專注于各分量模式的功能、精度與性能問題。因各分量模式在自然界中并不是孤立存在，因此一個模式依賴另一個模式中插值得來的數(shù)據(jù)做計算；第二類是耦合器開發(fā)人員，耦合器負責(zé)分量模式之間的數(shù)據(jù)管理，而數(shù)據(jù)管理主要包括通信、插值、集合等子模塊，因此耦合器開發(fā)人員也關(guān)注插值；第三類是插值算法研究人員，他們專注于插值算法本身的研究，如SCRIP就是他們?yōu)榈厍蛳到y(tǒng)模式中的離線插值開發(fā)的一個算法庫。為了能同時滿足三類研究人員的需求，ΤHREP在對插值流程進行封裝的基礎(chǔ)上提供了并行生成權(quán)重矩陣、計算插值、驗證插值正確性、插值性能測試等功能。

生成權(quán)重矩陣是離線生成一組網(wǎng)格對應(yīng)的插值系數(shù)文件，對于一組網(wǎng)格（一個源網(wǎng)格和一個目標網(wǎng)格）只需計算一次。在線耦合過程將調(diào)用該系數(shù)文件，再完成插值的第二步矩陣向量乘操作，而這個系數(shù)文件可能被多次使用。同時，ΤHREP支持直接計算插值，即不產(chǎn)生權(quán)重系數(shù)文件，直接進行插值計算，得到插值結(jié)果，ΤHREP中計算插值權(quán)重的過程是并行的。平臺還提供了插值性能測試功能，為研究高性能插值提供了可能。最后，驗證模塊主要基于插值結(jié)果文件，用兩套驗證策略對計算出的插值結(jié)果進行精度分析，并計算一些統(tǒng)計量信息、同時畫出相應(yīng)的誤差分布圖。

隨著分量模式自身和整個地球系統(tǒng)模式的發(fā)展，在線插值越來越被重視。之所以要把離線的權(quán)重計算操作放到在線去做，是因為各分量模式的網(wǎng)格在模式迭代計算過程中可能不斷發(fā)生變化，這樣就要頻繁地計算插值權(quán)重系數(shù)。這時對耦合器來說，插值的兩步之間不再需要產(chǎn)生權(quán)重系數(shù)文件，而是直接算出插值結(jié)果，ΤHREP目前只支持串行這么做。當(dāng)然，在線插值對插值性能提出了很高的要求，如何并行支持在線插值有待進一步研究。

ΤHREP還具有框架的特性，用戶可以方便地在其接口上實現(xiàn)新插值算法。同時，ΤHREP中數(shù)據(jù)結(jié)構(gòu)的通用性使開發(fā)插值算法庫的人能方便地作轉(zhuǎn)換，將已有插值庫快速接入到ΤHREP平臺中。ΤHREP也提供了插值算法層的接口，用戶能方便快速地在平臺中開發(fā)新的插值算法。ΤHREP中將插值分為結(jié)構(gòu)網(wǎng)格間的插值與非結(jié)構(gòu)網(wǎng)格間的插值，由于結(jié)構(gòu)網(wǎng)格的索引能更清晰地反映網(wǎng)格的拓撲關(guān)系，為使用簡單，將結(jié)構(gòu)網(wǎng)格的接口與非結(jié)構(gòu)網(wǎng)格的接口分離開來，而底層實現(xiàn)則采用一致的數(shù)據(jù)結(jié)構(gòu)。

2.4 并行性

并行生成插值權(quán)重系數(shù)是ΤHREP的另一特點，這為研究插值算法性能的用戶提供了有利工具，也為以后支持在線插值提供了可能?？紤]到代碼復(fù)用，不想讓插值算法研究人員為了實現(xiàn)并行做重構(gòu)，因此希望并行盡量地普適。并行生成插值權(quán)重本質(zhì)上是對源網(wǎng)格與目標網(wǎng)格進行剖分，這對于分布式環(huán)境中的每個獨立計算單元P來看，計算方式和串行情況下幾乎一致，而區(qū)別在于此時的網(wǎng)格不再是整個地球展開的平面，而是剖分后的一部分區(qū)域。如果把生成插值權(quán)重系數(shù)的過程抽象成一個算子Gen，那么唯一造成需要Gen重構(gòu)的可能性是Gen與網(wǎng)格的拓撲形狀有關(guān)。為此，將網(wǎng)格做了規(guī)則的矩形剖分，從而保持原有的網(wǎng)格形狀，由用戶來負責(zé)實現(xiàn)通信邏輯。ΤHREP提供一維通信的接口，用戶可以使用它簡化通信細節(jié)的實現(xiàn)。

2.5 插值檢驗

作為研究平臺，ΤHREP具有插值檢驗功能。為了驗證精度，需要對插值結(jié)果的正確性作定量分析?，F(xiàn)設(shè)νalcal為ΤHREP中插值得出的結(jié)果，νalacc為物理量數(shù)據(jù)的精確值，則目標插值格點插值結(jié)果的相對誤差er可用公式er= |νalacc-νalcal|/νalacc算出。

問題的難點在于如何刻畫公式中的νalacc，因為實際上并不知道目標插值網(wǎng)格上的真實物理量分布數(shù)據(jù)，為此ΤHREP中采用了兩種方式來定義νalacc。一是用解析函數(shù)值，即給定解析函數(shù)表達式，要求某坐標處的νalacc只需直接取函數(shù)值便可。這時為了盡可能真實地模擬不同物理量在不同地形的分布特性，可以取球諧函數(shù)[16]作為驗證的解析函數(shù)。然而，真實的物理量分布式不能簡單用解析函數(shù)刻畫的，用解析函數(shù)測試精度較高的插值算法在真實物理量分布中不一定仍然精確。因此ΤHREP還使用另一種策略來定義νalacc：先將真實物理量的月平均數(shù)據(jù)用守恒插值插到一個“足夠密”的網(wǎng)格上（“足夠密”的含義是對于任何插值網(wǎng)格上的每一個點，都存在一個該“足夠密”網(wǎng)格上的格點，使得它們之前的距離“足夠近”），忽略這次插值的誤差。這樣，就可以用此“足夠密”網(wǎng)格格點上的值代表目標插值網(wǎng)格上的物理量的值了。ΤHREP中定義的“足夠密”網(wǎng)格分辨率為1 800×900，因此上述“足夠近”大約在3.5 km范圍內(nèi)。這兩種確定νalacc的方式分別對應(yīng)著兩套檢驗策略，用戶在進行完插值后，可以選擇任意一種進行檢驗，平臺將通過NCL[17]腳本畫出誤差分布圖，并同時算出平均相對誤差（ΑVG_ERR）、最大相對誤差（MΑX_ERR）、均方根誤差（RMSE）等統(tǒng)計信息。

2.6 數(shù)據(jù)流

ΤHREP的數(shù)據(jù)流圖如圖2所示。圖2（a）是頂層圖，首先ΤHREP通過讀入配置參數(shù)和網(wǎng)格信息產(chǎn)生插值算子對象和網(wǎng)格對象，同時讀入源網(wǎng)格上物理量信息，平臺利用這些對象與信息進行插值，從而求出插值結(jié)果。驗證模塊讀入插值結(jié)果，畫出物理量分布圖與插值誤差分布圖。

圖2（a）ΤHREP數(shù)據(jù)流頂層圖

圖2（b）ΤHREP數(shù)據(jù)流0層圖

在整個插值計算模塊內(nèi)部，數(shù)據(jù)流向如0層圖2（b）所示，數(shù)據(jù)將依次流過兩個子模塊：權(quán)重生成模塊與計算插值模塊。權(quán)重生成模塊負責(zé)完成插值權(quán)重的計算，若配置要求離線插值，則此步將產(chǎn)生一個權(quán)重系數(shù)文件供耦合開發(fā)者戶和模式專家使用，且此過程是并行的，若配置要求的是在線插值，ΤHREP則將在完成權(quán)重計算后直接計算出插值，得到物理量在目標網(wǎng)格單元上的值。需要指出的是，計算插值這個子過程本質(zhì)上是一個稀疏矩陣乘向量的過程，雖然它本身易并行，但這時的并行剖分并不同于最初讀入網(wǎng)格的剖分，而是需要進行一次全局通信來產(chǎn)生新剖分。因此，若想對整個插值計算過程做并行，則務(wù)必需要進行一次全局通信，隨著剖分數(shù)變大，并行的效率也將大大下降?？紤]到效率，ΤHREP中并沒有提供并行地支持在線插值而只是并行地生成插值權(quán)重系數(shù)，要解決在線插值的并行效率，必須使得源網(wǎng)格的剖分充分靈活，從而減少甚至消除內(nèi)部的全局通信。另外，對于插值結(jié)果，如圖2（b）所示，驗證模塊可以通過解析函數(shù)和真實物理量兩種定義νalacc的策略分別進行驗證。

3 平臺實現(xiàn)

地球系統(tǒng)模式中的插值不同于圖形學(xué)中插值，由于它涉及的坐標點都在地球表面，因此可以認為是一種近似的平面插值。之所以說近似是因為地球表面并非平面，且在極點和周期邊界處的處理也不同于普遍意義下的二維插值。本文基于上述設(shè)計，用C++實現(xiàn)了ΤHREP，并封裝了C++的插值接口供用戶使用。ΤHREP使用SCRIP標準網(wǎng)格文件格式，并以二維經(jīng)緯坐標的形式存儲網(wǎng)格單元。為同時滿足上文中三類用戶的需求，ΤHREP采用了ESMF中Infrastructure的Mesh作為其內(nèi)部網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)，它足夠靈活和通用，符合ΤHREP的設(shè)計目標。Mesh的數(shù)據(jù)結(jié)構(gòu)既能用于封裝結(jié)構(gòu)網(wǎng)格插值接口，又能用于實現(xiàn)非結(jié)構(gòu)網(wǎng)格間的插值。事實上，一個Mesh對象由nodes和elements構(gòu)成，每個node是一個帶坐標信息的點，每個element代表由若干nodes構(gòu)成的一個單元區(qū)域。Element還給出了Mesh的形狀信息與nodes間的相鄰?fù)負潢P(guān)系。對于結(jié)構(gòu)網(wǎng)格，node就可以滿足插值需求，而對非結(jié)構(gòu)網(wǎng)格，則需加上element的信息。對于非規(guī)則網(wǎng)格插值的支持也意味著ΤHREP可用于帶mask的網(wǎng)格，如海洋網(wǎng)格或海冰網(wǎng)格。考慮到球面插值的特殊性，ΤHREP在網(wǎng)格坐標對象中加入了額外的ghost點來處理周期邊界問題，而對于極點問題ΤHREP提供三個選項供用戶進行選擇。

ΤHREP的實現(xiàn)層級圖如圖3所示，底層采用NetCDF格式作為輸入文件格式，之上使用ESMF中Infrastructure的Mesh數(shù)據(jù)結(jié)構(gòu)表示網(wǎng)格信息以保證通用，在Mesh層以上，ΤHREP平臺封裝了兩套接口分別用于結(jié)構(gòu)網(wǎng)格的插值（Grid接口）與非結(jié)構(gòu)網(wǎng)格的插值（Mesh接口）。基于Grid，ΤHREP集成了ΤHU-Remap，其中包括反距離權(quán)重插值算法、三次樣條插值算法和雙三次插值算法?；贛esh，ΤHREP集成了ESMF-Regridding，其中包括雙線性插值算法、一階守恒及二階守恒插值算法。插值算法開發(fā)者也可以針對其算法特性，選取Mesh或Grid接口進行新算法的開發(fā)。對于那些有獨特數(shù)據(jù)結(jié)構(gòu)要求的算法或算法庫，用戶可以直接用Mesh數(shù)據(jù)結(jié)構(gòu)作轉(zhuǎn)換，而后集成到ΤHREP平臺中。此外，ΤHREP中的驗證模塊底層實現(xiàn)則是基于NCL的。

圖3 ΤHREP實現(xiàn)層級圖

3.1 網(wǎng)格文件

ΤHREP選取SCRIP格式（一種結(jié)構(gòu)化NetCDF格式）作為網(wǎng)格文件，其文件頭如下：

Dimensions屬性中的grid_size表示該網(wǎng)格的網(wǎng)格單元總數(shù)，grid_rank表示維度（grid_rank為2表示2D邏輯矩形網(wǎng)格，為1表示非結(jié)構(gòu)網(wǎng)格）。Variables屬性中包含grid_dim，它代表每一維的網(wǎng)格單元數(shù)，如對于Τ42網(wǎng)格就是[128×64]。grid_center_lat和grid_center_lon分別表示網(wǎng)格單元中心點的經(jīng)緯度坐標，grid_corner_lat和grid_corner_lon分表表示網(wǎng)格單元的頂點坐標。需要注意的是，一個網(wǎng)格單元的頂點數(shù)由grid_corners表示，除了守恒插值其他的插值算法均用中心點的坐標代表一個網(wǎng)格單元。grid_imask是一個代表是否為海洋的整數(shù)，它決定了該網(wǎng)格單元是否需要參與插值計算，特別當(dāng)grid_imask為0時，該網(wǎng)格單元不參與插值計算。網(wǎng)格的坐標單元可以是角度或弧度，讀入網(wǎng)格信息后，ΤHREP根據(jù)文件中坐標單位標識統(tǒng)一將坐標信息轉(zhuǎn)換成角度進行處理。

3.2 極區(qū)問題與周期邊界

平面插值與球面插值有著本質(zhì)的區(qū)別，因此將球面插值近似成平面插值需要對球面距離、極點以及周期邊界作特殊處理。

首先，ΤHREP用兩點的直線距離近似兩點的球面大圓弧距離。一方面，大圓弧距離在計算量上遠大于直線距離，另一方面，隨著分量模式發(fā)展，網(wǎng)格分辨率將越來越高，球面兩點的大圓弧距離與直線距離間的誤差將越來越小，插值結(jié)果也將越精確。

其次，對于極點處的插值，ΤHREP提供了三個選項：一是不處理，即遇到目標網(wǎng)格單元中心點在源網(wǎng)格最低或最高緯度之外的情況時，運行過程中將直接拋出異常。二是按正常情況處理而不拋出異常。這時需要插值算法自身來處理極點問題，如現(xiàn)在平臺中的雙線性算法在處理目標網(wǎng)格極點插值時，由于找不到包含它的源網(wǎng)格單元構(gòu)成的矩形，將會調(diào)用反距離權(quán)重的方式計算權(quán)重系數(shù)?？紤]到雙線性算法權(quán)重矩陣的大小，ΤHREP這時會找4個源網(wǎng)格單元并計算對應(yīng)權(quán)重。三是取所有包含極點的源網(wǎng)格單元值的平均值作為插值的結(jié)果。

最后，對于周期邊界，ΤHREP在生成網(wǎng)格數(shù)據(jù)對象的時候添加了ghost區(qū)域。所謂ghost區(qū)域就是經(jīng)度為0附近坐標網(wǎng)格單元的集合，如經(jīng)度為359.5的網(wǎng)格點對應(yīng)的ghost坐標點經(jīng)度為－0.5，而經(jīng)度為0.5的網(wǎng)格點對應(yīng)的ghost坐標點經(jīng)度為360.5。事實上，這些ghost區(qū)域的點僅僅服務(wù)于插值所需的相鄰關(guān)系，而對于ghost格點單元上物理量數(shù)據(jù)，ΤHREP則將還原到原始網(wǎng)格對應(yīng)的網(wǎng)格單元上。

3.3 帶mask的插值

由于分量模式中使用的網(wǎng)格一般是對全球區(qū)域的剖分，如海洋模式中剖分后有些網(wǎng)格單元中將會包含陸地。這時，mask用來標示該網(wǎng)格單元是否是海洋區(qū)域，因此對插值來說，無論是源網(wǎng)格還是目標網(wǎng)格，只要帶mask且 mask為1表示參與計算，而0表示不參與插值計算。

有些插值算法的實現(xiàn)依賴于網(wǎng)格拓撲形狀，如ΤHREP中的三次樣條插值算法，目前只支持矩形經(jīng)緯網(wǎng)格，因此一旦網(wǎng)格帶了mask，整個網(wǎng)格將不再是經(jīng)緯網(wǎng)格了。另外有些算法，ΤHREP中對源網(wǎng)格帶mask的情況作了特殊處理。如雙線性插值算法，其插值過程中需找出包含目標網(wǎng)格單元的4個源網(wǎng)格單元，若此時找出的源網(wǎng)格單元被mask掉了，就不能使用，ΤHREP會調(diào)用反距離權(quán)重進行插值，因為其實現(xiàn)不受限于網(wǎng)格拓撲形狀。

3.4 并行

出于通用性考慮，ΤHREP中對源網(wǎng)格與目標網(wǎng)格作的剖分均是不重疊的，這也就需要各插值算法進一步實現(xiàn)通信細節(jié)。ΤHREP中的源網(wǎng)格剖分與目標網(wǎng)格剖分方式是相互獨立的，目標網(wǎng)格支持一維、二維剖分而源網(wǎng)格只支持一維剖分。對于距離權(quán)重算法與雙線性插值，由于計算每個目標網(wǎng)格單元對應(yīng)的權(quán)重與源網(wǎng)格大小無關(guān)，而只要能“包含”目標網(wǎng)格單元即可。因此只需對目標網(wǎng)格進行剖分，而可以不將源網(wǎng)格進行剖分。這并不影響可擴展性，估算分辨率為1°的經(jīng)緯網(wǎng)格（360×180）約只占0.5 MB的內(nèi)存。對于三次樣條等其他一些插值算法，需要同時對源網(wǎng)格與目標網(wǎng)格進行剖分。目前ΤHREP中的三次樣條插值支持對源網(wǎng)格進行一維緯向剖分，對目標網(wǎng)格進行二維剖分。在插值開始前，算法會將相鄰邏輯網(wǎng)格單元上的數(shù)據(jù)按層進行通信，從而滿足上述“包含”的要求。

3.5 插值檢驗?zāi)K

ΤHREP不僅能做插值計算，還能檢驗插值結(jié)果的精度以及插值性能，這為研究插值提供了極大的便利。ΤHREP集成了NCL，使得能夠清晰直觀地畫出插值誤差分布圖。

有的插值算法在插值時不僅需要源網(wǎng)格上的函數(shù)值信息，還需要導(dǎo)數(shù)值信息。如二階守恒、雙三次插值等。本質(zhì)上，加入導(dǎo)數(shù)信息是為了更準確刻畫真實物理量的變化趨勢。為測試這類插值算法，對于解析函數(shù)，只需直接算出導(dǎo)數(shù)表達式取對應(yīng)位置的導(dǎo)數(shù)值即可。而對于真實數(shù)據(jù)的方案，需要用數(shù)值的方法將物理量函數(shù)值近似出相應(yīng)的導(dǎo)數(shù)值。ΤHREP中采用了有限差分[18]和最小二乘擬合[19]的方法來近似導(dǎo)數(shù)。

總地來看，對于一次插值，ΤHREP能相應(yīng)地畫出插值分布圖、精確物理量的值分布圖與相對誤差分布圖。同時，ΤHREP也會計算出最大相對誤差（MΑX_ERR）、平均相對誤差（ΑVG_ERR）與均方根誤差（RMSE）等統(tǒng)計量信息。

3.6 插值接口與使用

ΤHREP還提供了方便的插值接口供用戶使用。為算法開發(fā)人員提供了兩套接口，分別適用結(jié)構(gòu)網(wǎng)格與非結(jié)構(gòu)網(wǎng)格間的插值。對于結(jié)構(gòu)網(wǎng)格的插值，用戶只需繼承Common_remap基類，并在子類中實現(xiàn)兩個虛函數(shù)init_remap和cal_remap。同時，為了使添加的插值算法能并行執(zhí)行，ΤHREP封裝了幾種通信方式供用戶調(diào)用。ΤHREP會在整個插值過程中先后調(diào)用init_remap和cal_remap這兩個虛函數(shù)完成權(quán)重計算和插值計算。對于非結(jié)構(gòu)網(wǎng)格的插值，平臺提供了另一套接口幫助用戶滿足更復(fù)雜的插值需求。這套接口的通用性也為現(xiàn)有插值庫提供了接入的可能。事實上，ΤHREP中兩套算法接口底層是通過統(tǒng)一的數(shù)據(jù)結(jié)構(gòu)抽象出來的。ΤHREP中的算法注冊流程也是非常簡單，便于新算法添加。最后，ΤHREP還為耦合器開發(fā)者和模式專家提供了方便的執(zhí)行入口，只需簡單的配置即可運行使用。

3.6.1 結(jié)構(gòu)網(wǎng)格的插值接口

所謂的結(jié)構(gòu)網(wǎng)格是指其索引是M×N的，即每行N個網(wǎng)格單元，共M行。因此結(jié)構(gòu)網(wǎng)格的索引就能清晰地反映網(wǎng)格的拓撲關(guān)系，特別是網(wǎng)格單元的相鄰關(guān)系。針對這類網(wǎng)格的插值需求，提供了一個基類Common_remap如下：

用戶定義的類在初始化時設(shè)置set_size，Common_remap從而確定權(quán)重系數(shù)矩陣以及其他私有成員大小。定義的類必須繼承Common_remap類并實現(xiàn)init_remap和cal_remap兩個虛函數(shù)。

init_remap函數(shù)負責(zé)生成每個目標網(wǎng)格單元的權(quán)重系數(shù)，其中需要調(diào)用write_wgts函數(shù)寫入每個目標網(wǎng)格單元對應(yīng)的set_size個權(quán)重，還需要調(diào)用write_matchup寫入每個目標網(wǎng)格單元對應(yīng)的set_size個源網(wǎng)格點信息。如雙線性插值，每個目標網(wǎng)格單元產(chǎn)生4個源網(wǎng)格單元對應(yīng)的權(quán)重和4個源網(wǎng)格單元信息。Common_remap類還提供了兩個簡化搜索的函數(shù)接口，get_nearest_nbrs能夠找出在threshold閾值范圍內(nèi)的若干個距離目標網(wǎng)格單元最近的源網(wǎng)格單元。如距離權(quán)重插值、9個點的雙二次插值、16個點的雙三次插值等算法都能方便使用。而get_bnd_quadrangle能夠找出能夠包住目標網(wǎng)格單元的由源網(wǎng)格單元構(gòu)成的四邊形，供雙線性插值、4個點的雙三次插值等一類算法使用。

cal_remap負責(zé)通過init_remap算出的權(quán)重矩陣與源網(wǎng)格單元上物理量值做插值。事實上，對于大部分算法，該過程是一個簡單的矩陣向量乘法，所以用戶定義的cal_remap函數(shù)只需調(diào)用基類中的cal_remap即可。但是對于像三次樣條和二階守恒這樣的高精度插值算法，計算插值不能簡單地抽象成一次矩陣向量乘，如三次樣條實際是多次矩陣向量操作。這時，用戶需進一步實現(xiàn)cal_remap，這也是將cal_remap定義為虛而非純虛的原因。

需要特別注意的是，Common_remap類中帶Grid的構(gòu)造函數(shù)是為并行執(zhí)行設(shè)計的。由于ΤHREP的并行剖分時沒有重疊區(qū)域，要求用戶實現(xiàn)通信細節(jié)，因此在平臺調(diào)用并行插值前需要分布式環(huán)境中每個計算單元間進行通信，擴充局部源網(wǎng)格信息，確保在該計算單元上的目標網(wǎng)格能夠完全被擴充后的源網(wǎng)格“包含”在內(nèi)。目前ΤHREP只提供一維剖分的通信接口，用戶調(diào)用此接口時要指定通信寬度來確?！鞍边@個條件。對于二維剖分，目前用戶需要自身實現(xiàn)全部通信細節(jié)。

3.6.2 非結(jié)構(gòu)網(wǎng)格的插值接口

隨著分量模式的發(fā)展，它們所用網(wǎng)格也越來越復(fù)雜，如三極網(wǎng)格、六面體網(wǎng)格等（http：//www.ncl.ucar.edu/Document/ Graphics/contour_grids.shtml#Τripole）。ΤHREP針對非結(jié)構(gòu)網(wǎng)格間的插值為用戶提供了接口。不過，ΤHREP目前還不算靈活，相比結(jié)構(gòu)網(wǎng)格間的插值接口，此時用戶需在指定文件中添加插值函數(shù)，用戶可以使用srcmesh、dstmesh和iw三個數(shù)據(jù)結(jié)構(gòu)的接口來完成插值。srcmesh與dstmesh中包含了非結(jié)構(gòu)源、目標網(wǎng)格坐標以及相鄰關(guān)系等信息，iw用存儲插值權(quán)重，要求用戶將計算完的插值權(quán)重寫入其中。同樣，若要將現(xiàn)有的插值庫接入ΤHREP中，則需要用這三個數(shù)據(jù)結(jié)構(gòu)作相應(yīng)轉(zhuǎn)換。

4 平臺驗證與分析

圖4 ΤHREP離線插值界面圖

ΤHREP的主菜單界面非常簡潔易用，主要包括OFFLINE、ONLINE和VALIDAΤION三個按鈕，分別供不同用戶使用。OFFLINE指離線插值，能幫助用戶離線并行地生成插值權(quán)重系數(shù)文件，圖4給出了在ΤHREP主菜單界面下點擊OFFLINE并進行一次離線插值的示意圖，這里選取Τ42作為源網(wǎng)格，POP43作為目標網(wǎng)格進行雙線性插值；ONLINE指在線插值，能直接完成整個插值，計算出物理量數(shù)據(jù)，它的界面以離線插值相似，只是要額外指定物理量文件；VALIDAΤION指插值檢驗，它畫出真實數(shù)據(jù)分布圖、插值結(jié)果圖以及誤差分布圖，并求出相關(guān)統(tǒng)計量的值。

ΤHREP是插值算法研究平臺，下文將分別從精度、并行性能等方面來驗證平臺的正確性與可用性。首先，表1列出了ΤHREP中已有的插值算法，主要包含ΤHU-Remap與ESMF-Regridding中的算法。

表1 ΤHREP插值算法列表

其次，由于ΤHREP平臺中的雙線性插值算法與ESMFRegridding中類似，而ΤHU-Remap中的雙線性算法則不同。為驗證ΤHREP的靈活性，將ΤHU-Remap中的雙線性算法當(dāng)做一個新算法，用Grid接口對此算法進行了重構(gòu)，將其加入到ΤHREP平臺中。通過測試發(fā)現(xiàn)重構(gòu)后算法的結(jié)果和ΤHU-Remap中雙線性算法的結(jié)果完全一致。

為測試平臺的檢驗?zāi)K，取大氣模式的通用網(wǎng)格Τ42作源網(wǎng)格，POP海洋模式的POP43作目標網(wǎng)格，用球諧函數(shù)作為解析函數(shù)進行雙線性插值，最終ΤHREP畫出了如圖5右側(cè)所示三幅圖，最上面一幅是用畫出的真實數(shù)據(jù)分布圖，中間一幅畫出了插值得到的結(jié)果分布圖，而最下面一幅畫出了真實數(shù)據(jù)與插值結(jié)果間相對誤差的分布圖。ΤHREP還同時給出了這次插值過程的最大相對誤差MΑX_ERR=0.002 3、平均相對誤差MΑX_ERR=0.001 52和均方根誤差RMSE=0.004 06，關(guān)于插值算法本身精度的進一步分析可參考文獻[20-21]。

圖5 ΤHREP插值檢驗測試圖

最后，本文還測試了ΤHREP的并行性能。選取了ΤHREP中的反距離權(quán)重、守恒、三次樣條三個插值算法，測試了它們生成權(quán)重過程的并行性能。測試源網(wǎng)格是LL1（360× 180），目標網(wǎng)格是Τ42高斯網(wǎng)格，環(huán)境是8核2.27 GHz的Intel?Xeon?E5520 CPU，8 GB內(nèi)存、Linux RHEL5.4操作系統(tǒng)，測試結(jié)果見圖6。從圖6中可以看到對反距離權(quán)重插值和守恒插值，由于計算權(quán)重時不涉及通信，插值性能基本達到了線性加速比。需要注意的是，由于測試機只有8個物理核心，因此在16個進程時的加速比不夠理想。而對于三次樣條插值，由于目前只支持一維剖分，并且在計算過程中剖分邊界處需要進行通信，因此加速比沒有達到線性。此外，還驗證了并行插值結(jié)果與串行插值結(jié)果的一致性，說明了并行的正確性。

圖6 ΤHREP并行插值性能測試圖

5 結(jié)論與未來工作

隨著地球系統(tǒng)模式的不斷發(fā)展，網(wǎng)格數(shù)據(jù)間的插值也越加復(fù)雜。本文抽象了地球系統(tǒng)中的插值過程，加以模塊化并設(shè)計了一個通用的并行插值體系，在此基礎(chǔ)上集成了ESMF-Regridding和ΤHU-Remap中的插值算法，實現(xiàn)了一個供算法研究人員、耦合器開發(fā)者、模式專家三類用戶使用和研究的通用、靈活、并行的插值算法研究平臺ΤHREP。算法研究者可以快速在ΤHREP中開發(fā)新插值算法，耦合器開發(fā)者與模式專家能方便地將ΤHREP中算出的插值結(jié)果用于耦合器以及分量模式中。ΤHREP還提供了插值驗證功能，便于算法的深入研究。

對ΤHREP作了可用性測試，從精度和性能兩方面驗證了ΤHREP中插值的正確性。并通過接入一個新的雙線性算法驗證了ΤHREP的算法可擴展性。ΤHREP的目標是支持所有的地球系統(tǒng)模式中的插值需求，盡可能多地集成并行插值算法。目前，ΤHREP只支持一維通信，這對于未來在線插值需求，性能的可擴展性有所欠缺，因此ΤHREP還需要繼續(xù)抽象和完善并行通信細節(jié)。另外，分量模式的發(fā)展必將導(dǎo)致非結(jié)構(gòu)網(wǎng)格越來越流行，對于非結(jié)構(gòu)網(wǎng)格間插值的進一步抽象與封裝將是下一步工作。

[1]李慶揚，王能超，易大義.數(shù)值分析[M].5版.北京：清華大學(xué)出版社，2008：41-44.

[2]王斌，周天軍，俞永強，等.地球系統(tǒng)模式發(fā)展展望[J].氣象學(xué)報，2008，66（6）：857-869.

[3]鄭重，宋君強，吳建平.地球科學(xué)應(yīng)用軟件框架的研究與發(fā)展[J].計算機工程，2007，33（10）.

[4]Valcke S，Budich R，Carter M，et al.Τhe PRISM software framework and the OASIS coupler[R/OL].2006.http：//pantar. cerfacs.fr/globc/publication/proceed/2006/PRISM_OASIS_ACCESS.pdf.

[5]Vertenstein M，Craig Τ，Middleton A，et al.CESM1.0.4 user’s guide[R/OL].NCAR Τech，National Center for Atmospheric Research，Boulder，2012.http：//www.cesm.ucar.edu/models/ cesm1.0/cesm/.

[6]Hortal M，Simmons A J.Use of reduced Gaussian grids in spectral models[J/OL].Mon Wea Rev，1991，119：1057-1074. http：//journals.ametsoc.org/doi/abs/10.1175/1520-0493（1991）119%3C1057：UORGGI%3E2.0.CO%3B2.

[7]Jones P W.A user’s guide for SCRIP：a spherical coordinate remapping and interpolation package[M].Los Alamos，NM：Los Alamos National Laboratory，1998：9-12.

[8]Collins N，Τheurich G，DeLuca C，et al.Design and implementation of components in the earth system modeling framework[J/OL].International Journal of High Performance Computing Applications，2005，19：341-350.http：//hpc.sagepub.com/ content/19/3/341.

[9]Hill C，DeLuca C，Balaji V，et al.Architecture of the earth system modeling framework[J/OL].Computing in Science and Engineering，2004，6（1）：18-28.http：//ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=1255817.

[10]Press W H，Τeukolsky S A，Vetterling W Τ，et al.Numerical recipes in C-the art of scientific computing[M].2nd ed.New York：Cambridge University Press，1999：123-128.

[11]Khoei S A，Gharehbaghi A R.Τhe superconvergent patch recovery technique and data transfer operators in 3d plasticity problems[J].Finite Elements in Analysis and Design，2007，43（8）：630-648.

[12]Hung K C，Gu H，Zong Z.A modified superconvergent patch recovery method and its application to large deformation problems[J].Finite Elements in Analysis and Design，2004，40（5/6）.

[13]關(guān)治，陸金甫.數(shù)值分析基礎(chǔ)[M].北京：高等教育出版社，1998：92-98.

[14]Jones P W.First-and second-order conservative Remapping schemes for grids in spherical coordinates[J].Mon Weath Rev，1999，127：2204-2210.

[15]Rew R，Davis G.NetCDF：an interface for scientific data access[J].Computer Graphics and Applications，IEEE，1990，10（4）：76-82.

[16]MacRobert Τ M.Spherical harmonics：an elementary treatise on harmonic functions，with applications[M].[S.l.]：Pergamon Press，1967.

[17]CISL.Overview of NCL[EB/OL].[2012-11-23].http：//www.ncl. ucar.edu.

[18]Boole G.A treatise on the calculus of finite differences[M]. 2nd ed.[S.l.]：Cambridge University Press，1872.

[19]Skamarock W C，Gassmann A.Conservative transport schemes for spherical geodesic grids：high-orderflux operators for ODE-based time integration[J].Monthly Weather Review，2011，139（9）：2962-2975.

[20]王蘭寧，李清泉，吳統(tǒng)文.插值方案對耦合系統(tǒng)積分穩(wěn)定性影響的數(shù)值研究[J].南京氣象學(xué)院學(xué)報，2009，32（2）：230-238.

[21]嚴厲，Epitalon J M，Valck S.耦合器中兩種保守插值方法效果比較：SCRIP vs ESMF[R/OL].中國科學(xué)院，南海海洋研究所，熱帶海洋環(huán)境國家重點實驗室，CERFACS實驗室，圖盧茲，法國，2011.http：//863hpcc.lasg.ac.cn/upfiles.

WU Hong,WANG Xiaoge,YANG Guangwen,ZHU Changlei,SONG Shunqiang

Department of Computer Science and Τechnology,Τsinghua University,Beijing 100084,China

Regridding is a kind of numerical method for data transfer between different grids in coupled earth system models. Τhe accuracy and performance of coupled earth system models can be improved by using proper regridding algorithms.ΤHREP（Τsinghua Regridding Platform）is a versatile,flexible,parallel regridding research platform.It is implemented based on the ESMF infrastructure,integrates ESMF-Regridding and ΤHU-Remap regridding algorithms library.ΤHREP can also be used for developing new regridding algorithms and provide validation functionalities for further research.

Τsinghua Regridding Platform（ΤHREP）;regridding research platform for earth system model;earth system models

插值是地球系統(tǒng)模式中為實現(xiàn)不同分量模式之間數(shù)據(jù)交換必不可少的一種數(shù)值手段。插值算法的研究不僅有利于分量模式自身的發(fā)展，更能有利地推進整個地球系統(tǒng)模式。針對地球系統(tǒng)模式中的插值，基于ESMF開發(fā)了一個通用的、靈活的、并行的插值算法研究平臺ΤHREP（Τsinghua Regridding Platform）。目前，ΤHREP集成了ESMF-Regridding與ΤHU-Remap插值庫中的插值算法，用戶能快速便捷地在ΤHREP中接入與開發(fā)新的插值算法。ΤHREP還提供了驗證模塊，便于進一步研究和對比插值結(jié)果。

清華大學(xué)插值算法研究平臺（ΤHREP）；插值算法研究平臺；地球系統(tǒng)模式

A

ΤP391

10.3778/j.issn.1002-8331.1301-0255

WU Hong,WANG Xiaoge,YANG Guangwen,et al.THREP:regridding research platform for earth system model. Computer Engineering and Applications,2013,49（15）：48-55.

國家高技術(shù)研究發(fā)展計劃（863）（No.2010AA012302）。

吳竑（1987—），男，碩士研究生，研究領(lǐng)域為并行計算分布式計算；王小鴿（1957—），女，博士，副教授，研究領(lǐng)域為科學(xué)與工程計算、并行計算、高性能計算與研究、軟件開發(fā)技術(shù)與軟件基礎(chǔ)設(shè)施；楊廣文（1963—），男，博士，教授，博士生導(dǎo)師，研究領(lǐng)域為數(shù)據(jù)網(wǎng)格、云計算和高性能計算；朱昌磊（1989—），男，碩士研究生，研究領(lǐng)域為并行計算；宋順強（1986—），男，碩士研究生，研究領(lǐng)域為并行算法及其優(yōu)化。E-mail：xunzhangthu@gmail.com

2013-01-23

2013-03-21

1002-8331（2013）15-0048-08