基于循環神經網絡的正交網格的自動化生成算法

黃中展，徐世明

（清華大學地球系統科學系，北京 100084）

（*通信作者電子郵箱xusm@tsinghua.edu.cn）

0 引言

動畫電影技術、天氣預報、工業制造等多個與人們生活息息相關的領域都離不開科學計算，隨著這些領域對高質量計算的需求的增大，對數值算法的要求也隨之增大。一般來說，在科學計算當中，主要采用的三類數值算法分別為有限差分方法［1］、有限體積方法［2］以及有限元方法［3］。它們廣泛地應用于各類計算當中，而這些方法的有效實施離不開高質量的離散化方案，即網格的劃分。

在不同的應用場景，網格有很多不同的類型，可從網格形狀、目標區域類型、目標區域的維度三個方面來劃分。如圖1（a）和圖1（b）分別展示了二維平面和三維空間中的正交網格樣例。

對于二維平面上的網格，除了在圖1（a）中所展示的四邊形網格之外，還有如圖2 所示的三角形網格。相較于四邊形網格，三角形網格的生成相對容易。另外，圖1（a）中的網格和圖2中的網格還有一個顯著的差異在于其目標區域的連通性上，圖2 所展示的目標區域是多連通區域。多連通區域相對于單連通區域而言，在實際生活和科學計算需要中更為普遍，但是同時網格生成難度也更大，特別是正交網格。

圖2 多連通區域的三角網格［6］Fig.2 Triangular grid of multiconnected region

對于平面上的目標區域，三角網格和四邊形網格的算法都均可以對目標區域進行離散化。三角形網格主要采用的是Delaunay 三角剖分法［7］、四叉樹法［8］以及波前法［9］，其中Delaunay 法使用最為普遍，這是由于它在數學上十分成熟，具有較好的網格質量且收斂性也能較好地保證；而四邊形網格的生成主要有格柵法［10］、拓撲分解法［11］、節點連接法［12］等。在天氣預報、氣候模擬等科學領域中，正交網格應用最為廣泛，本文主要考慮的是正交網格的自動化生成。

對于單連通區域的正交網格生成來說，利用Thompson 微分方程法［13］最為普遍，然而這類方法產生的網格在邊界處正交性保持較差，收斂性也不能保證，且具有較多的需要調節的參數，這使得當遇到較復雜的目標區域的時候，需要有較多的人工調節，不便于科學計算的實際應用。快速且自動化的生成算法在各種科學計算領域當中十分必要。

1 正交網格的自動化生成

接下來，介紹如何利用循環神經網絡和復分析知識來完成正交網格的自動化生成。首先，介紹正交網格生成所需要的基本知識，即共形映射和相關網格生成工具。然后利用這些方法將正交網格問題轉換為一個最優化問題，同時介紹利用長短期記憶（Long Short-Term Memory，LSTM）網絡［14］來降低最優化問題求解時的時間復雜度，從而能自動化地生成目標區域的正交網格。

1.1 Schwarz-Christoffel共形映射

在復分析中，斯瓦茨-克里斯托弗爾（Schwarz-Christoffel，SC）共形映射是多邊形區域中常用的保角映射［15］。SC映射是從一個簡單多邊形（即不自交的多邊形）區域到復平面中的上半平面H={ζ∈C：Imζ＞0}的一一映射方式。另外SC 映射具有顯式表達式，考慮n個頂點的簡單多邊形區域(vi，αi)，i=1，2，…，n，那么H到該區域的SC映射f可以由

給出。其中：C0和C1為常數，z1，z2，…，zn-1為復平面中實軸上的n-1 個實數，這些實數滿足如下對應關系f(∞)=vn和f(zj)=vj(j=1，2，…，n-1)。

共形映射具有一一對應且保持角度這兩個良好性質，這些性質可以用于正交網格的生成。如圖3 所示，若區域B為目標區域，即需要在區域B當中生成正交網格。可以首先利用SC 映射建立區域B到區域H的共形映射g。若可以建立一個簡單的、容易構造正交網格的區域，如區域A，同樣利用SC映射可以構建從區域A到區域H的共形映射f，注意到g具有一一映射的性質（具有可逆性），可以建立從區域A到區域B之間的共形映射f⊙g-1，而由保角性，區域A的正交網格可以保持到區域B中，由此得到目標區域B的正交網格。這就是利用共形映射構造目標區域正交網格的基本思路。

圖3 基于SC映射的正交網格生成Fig.3 Orthogonal grid generation based on SC mapping

1.2 Gridgen-c網格生成工具

基于SC 共形映射構造正交網格的理論相對完備，但是在實際網格生成應用以及算法上仍有很多不足。在眾多基于共形映射的網格生成工具中，Gridgen-c 工具是基于C 語言的正交網格生成工具，具有快遞且穩定的網格生成能力。其不足主要在與需要人為給定目標多邊形的轉角類型，利用這些人為標定的轉角類型，從給定的目標區域B中人為地構造出圖3中的區域A。具體來說，如圖4所示。

圖4 Gridgen-c網格生成工具Fig.4 Grid generation tool Gridgen-c

給定目標區域B，和三種轉角類型。Gridgen-c 需要人為提供先驗信息，如圖4 的區域B中，只需要令以及等（對應圖4 區域A的轉角類型），就能根據這些給定的轉角類型直接構造如圖4 所示的區域A。特別地，區域A中的多邊形區域實際上只需要通過簡單的經緯網（水平、豎直兩個方向）即可構造正交網格。利用圖3中的流程便可以完成正交網格生成。注意到，如果對于只有少數頂點個數的多邊形，通過人為地給定轉角類型，可以輕松地使用Gridgen-c 得到正交網格。然而在科學計算的實際具體問題中，如計算機圖形學中需要對一定區域進行流體模擬［17］；海洋科學中需要對一定海域做溫度計算［18］等，這些具體問題所考慮的多邊形區域往往具有較多的頂點個數。若此時，直接對這樣的多邊形進行人為的頂點轉角類型標定的話：一方面將消耗大量的時間；另一方面，對于多頂點的多邊形轉角標定時，由于人的能力有限，人為的先驗信息可能無法給出一個較優的轉角標定。針對這樣的問題，下面將結合Gridgen-c和深度學習的辦法給出一種正交網格的自動化生成算法。

1.3 正交網格生成對應的最優化問題

本節首先對正交網格生成問題建模，即將利用Gridgen-c所需的條件將正交網格的生成問題轉換為一個帶線性限制條件的整數規劃問題。

不妨設對于目標N多邊形(vi，αi)，三種轉角類型的個數分別為m，p以及q。于是由簡單多邊形的內角和公式有：

從式（2）可以化簡得到-q+m=4，即：

為敘述方便，不妨設N個關于多邊形每個頂點的新的變量xi(i=1，2，…，N)，其中對所有的xi滿足：

結合式（3）即知xi中分別有q，p，m個的大小為-1，0，1。

在使用Gridgen-c 工具生成正交網格的時候，人為地為目標多邊形區域標定每一個頂點的轉角類型的時候，必須滿足式（4）的限制。換言之，對于一個目標多邊形，其各轉角類型的選擇方案必然對應式（4）的一個解。規定好轉角類型后，接著著手設計最優化目標。可以從三個角度來考慮生成的網格的品質的好壞，分別是正交性、均勻性、覆蓋性。

正交性 顧名思義即為通過數值方法得到的網格，在各個網格點中正交的程度。從復分析的理論上看，由如圖3中f⊙g-1得到的區域B中的網格應該是嚴格正交的，但是由于數值誤差，算法穩定性等問題可能會導致在實際數值算法所生產的網格的正交性不是十分嚴格，特別是在邊界上。

均勻性 生成的網格應保持足夠均勻，否則容易產生極端的網格點分布，盡管這些網格點所構成的網格的正交性可能可以保持足夠好，但由于網格不均勻所導致的尺度差異，一定程度上不利于其實際的科學計算，可能會產生較大的數值誤差等問題。

覆蓋性 這里所說的覆蓋性是指生成的網格能足夠好地覆蓋目標區域。在使用Gridgen-c 工具產生正交網格時，目標多邊形的頂點轉角的錯誤的標定容易造成所生成的網格的區域較小無法覆蓋目標區域，從而不能足夠好地滿足科學計算中的計算需求。

結合式（4）和上述正交性、均勻性和覆蓋性的考慮，可以構造如下最優化問題：

其中Lo(xN)刻畫的是網格的正交性，即

cos?anglemax和cos?anglemin表示各個網格點中夾角（銳角）余弦值的最大和最小值。容易知道，當Lo(xN)越小時，各網格點處能較好地垂直。另外Lu(xN)表示網格的均勻性，即

areamax和areamin分別表示網格最大和最小面積，當Lu(xN)越小時，網格的最大和最小的面積越接近，即網格的面積大小越均勻。而Lc(xN)表示的是網格的覆蓋性，利用式（8）定義：

其中：Area(grid)表示的是生成網格的覆蓋面積，而Area(O)表示的是目標區域的面積。注意到，平面上的正交網格生成問題實際上即多邊形區域上的正交網格生成問題，均可以通過對式（5）進行求解獲得合適的正交網格。為了進一步驗證式（5）的通用性，在第2 章中將對三種不同類型的區域驗證算法的有效性。在算力允許的情況下，直接遍歷式（4）就可以求解該最優化問題。但是實際上，最優化問題（5）是一個帶線性限制條件（式（4））的整數規劃問題，這是一個NP-hard 的問題，理論上只能通過枚舉來獲得最優解，而遍歷式（4）就至少需要O(3N)的復雜度消耗，由于是指數級復雜度，當N稍大的時候就使得問題幾乎不可解。注意到最優化問題（5）的優化目標并不具有良好的可導性和連續性，這使得一些有效的整數規劃問題的求解器［19］不能用來對其進行求解。

針對最優化問題（5）時間復雜度過高的問題，提出了利用循環神經網絡構建一個關于轉角選擇的分類器來降低時間復雜度，從而獲得較優近似解的辦法。

1.4 循環神經網絡

為了可以利用循環神經網絡來降低最優化問題（5）求解的時間復雜度，此處提出一個基本假設：

基本假設 一個頂點的轉角類型取決于其鄰近點的轉角類型，且越接近的頂點影響越大。

科學計算中所考慮的目標區域多是自然形成的邊界，如海岸線、湖泊、動畫人物等構成的多邊形。這些多邊形的邊界較為自然，若不滿足基本假設，那么它的頂點轉角相對獨立，邊界一般不太符合自然規律。若基本假設成立，那么對于目標多邊形的任意一個頂點vi，其轉角類型xi取決于其附近的多邊形的局部信息，即由vi+1，vi+2，…，vi+t以 及vi-1，vi-2，…，vi - t這些頂點序列構成的某種特征來決定。

接下來利用頂點的兩類信息來構造特征：一方面是角度信息，用A(vi)表示頂點vi的內角的余弦值，即向量和向量構成的夾角（內角）的余弦值；另一方面是長度信息，即向量和向量的模長的均值，用L(vi)來表示。但是注意到，實際上L(vi)是無界的，需要進行歸一化，即

由此，對于每一頂點vi，都可以構造如下特征

接著考慮較小規模的多邊形，如N=10，實際上可以直接對最優化問題（5）進行枚舉求解獲得最優轉角類型。也就是說對于頂點個數為10 的多邊形，可以先由式（10）獲得特征集，然后通過枚舉直接求解得到每個頂點的最優轉角，然后利用深度學習的方法來建立特征集到最優轉角的映射。

長短期記憶（Long Short Term Memory，LSTM）網絡是循環神經網絡中最著名的擴展［14］，具有很強的捕獲序列之間的信息的能力。由基本假設，多邊形任一頂點的轉角類型取決于鄰近頂點的特征構成的序列，故此處適合采用LSTM 來學習特征集到轉角類型的映射。具體而言，如圖5 所示，對于N=10時，第t步的更新計算公式為：

其中，W、U和b是LSTM 各個門的參數。ct為第t步中LSTM 的記憶單元，ht為其隱含層的輸出，首先初始化ct和ht，即h0=0以及c0=0。接著，以頂點vi為例，LSTM 的第1步的輸入y1為，接著第2步的輸入y2為，類似地，第3 步的輸入y3為，第4 步輸入為，如此類推，一直到第10步。由此，對于N=10 的情況，可以由LSTM 訓練得到一個能判斷目標多邊形轉角類型的分類器，即給出各頂點在三種轉角類型的概率。注意到，由基本假設，實際上一個頂點的轉角類型只和其附近的頂點信息有關，這意味著對于一個M邊形，其中M≥N，需要判斷其中一個頂點的轉角類型，可以只利用上其鄰近的頂點信息，如10 個頂點的信息，由此可以用上在N=10時訓練出來的分類器來解決M邊形的網格生成問題。

圖5 LSTM結構Fig.5 Structure of LSTM

但是如果僅僅依靠分類器來進行轉角類型的判斷，由于分類器存在誤差，得到的解不一定能滿足式（4），即使能滿足，得到的網格質量不一定能足夠地好（取決于轉角分類的質量）。面對這樣的問題，可以引入一定量的枚舉，將LSTM 作為一種降低時間復雜度的工具。

1.5 降低時間復雜度

由于訓練數據或者算法的不足，導致LSTM 分類器的性能存在一定的誤差，直接使用LSTM 來求解最優化問題（5）很可能會有無解或者解的質量不高的問題。引入一定量的枚舉能夠緩解這兩個問題。

如圖6 所示，采用一種簡單的策略，對于一個N多邊形的一個頂點，首先可以利用由10 多邊形訓練得到的LSTM 分類器來對其頂點類型進行判斷，即得到三種轉角類型{0，-1，1}的概率P，給定閾值P0。用兩個例子來說明圖6 的算法流程，如P0=0.6，當轉角概率P={0.7，0.1，0.2}時，此時P最大的概率max(P)為0.7，大于P0，此時可以直接確定轉角類型，即此時轉角類型為0。如當P={0.5，0.3，0.2}時，此時最大的概率為0.5，比P0小，此時分類器沒有足夠的把握判斷目標多邊形的該頂點的轉角類型，于是，可以不考慮概率最小的類型，即類型1，僅僅只考慮較大概率的兩種轉角類型。

圖6 利用LSTM降低時間復雜度Fig.6 Reducing time complexity by LSTM

一般來說，最優化問題（5）的求解至少需要的O(3N)時間復雜度。若利用圖6 所示的策略，有L個頂點可以直接確定（即轉角類型的最大概率大于閾值P0），此時時間復雜度降為O(3N-L)，而剩下的頂點都不考慮概率最小的情況，此時時間復雜度進一步降低為O(2N-L)，一般來說當L稍大的時候，該計算消耗是可以承受的，L的大小取決于LSTM 的分類能力和閾值P0的選取。由此，通過圖6的策略，可以完成的N多邊形（其中N≥10）的正交網格的自動化生成。

1.6 訓練數據

實際上，為了得到足夠好的LSTM 分類器，需要較大數量的多邊形樣本，其中這些多邊形的邊界還需要滿足一定的物理意義。也就是說，訓練集中的多邊形不能有過于復雜的邊界，否則一方面可能最優化問題（5）沒有足夠好的解，另一方面較差的解會給分類器的訓練引入較大的噪聲，降低分類能力，不利于計算復雜度的降低。注意到GADM［20］數據庫中有全球各國各行政單位的地理邊界數據。該數據集數據量大，且這些地理邊界多具有較自然的邊界，較為符合基本假設的要求，恰好可以用作此處LSTM 分類器的訓練。為了得到合適的訓練集，至少做如下3個預處理。

1）保證所考慮的多邊形是簡單多邊形。如果訓練數據中含有非簡單多邊形容易產生較大的數值錯誤，對模型的訓練注入較大的噪聲。

2）保證考慮的多邊形是單連通區域。在GADM數據集中包含很多群島等地理邊界，這些邊界本身就構成了多連通的區域，這些數據應該被預先剔除。

3）由于LSTM 分類器的訓練需要的多邊形的頂點個數為10，所以對于GADM 數據集中的多邊形，可以對其頂點做采樣，獲得10 邊形。然后根據式（10）得到對應的特征，且利用式（5）直接枚舉出最優解，作為訓練數據。

2 實驗與模型分析

本章將通過四個樣例來觀察生成的正交網格的質量。首先考慮的是兩個簡單形狀的圖形，即如圖7（a）所示的每一個轉角都是90°或者270°的圖形。這樣的圖形能通過簡單的縱橫劃分（經緯網）來獲得正交網格。

接著，圖7（b）所示的圖形是正16 邊形，實際上它比較接近一個圓形。實際上，由于其對稱性，它的每一個點的轉角類型應該都是要相同的，然而由于式（4）的限制，它的各點轉角類型不可能完全一樣，所以這樣的目標多邊形區域本身不存在像圖7（a）這樣足夠好的解。從實驗結果來看，圖7（a）所示的圖形生成的正交網格恰好是通過縱橫劃分獲得的，達到7（a）所示圖形的網格生成的最優解。如果直接利用Gridgen-c工具和人工選擇轉角類型，最好的方式也是將90°的轉角設為轉角類型1，將270°的轉角設為轉角類型-1，也就是說本文的方法在此樣例上達到了最優解。另一方面，盡管正多邊形這樣的圖形在生成正交網格的問題上存在客觀的困難，但從圖7（b）所示的正16 邊形來看，生成的網格已經足夠地好，與直接人工進行轉角類型選擇所生成的正交網格一致。

圖7 正交網格生成樣例Fig.7 Examples of orthogonal grid generation

接下來是考察真實目標區域，首先是如圖7（c）所示的動畫形象區域。在計算機圖形學領域中，動畫領域的模擬和計算應用是重要的課題。而圖7（d）所示的是真實地理區域，以非洲大陸為例。在海洋科學和大氣科學等領域，如天氣預報、溫度預測、模擬等問題上需要高質量的正交網格劃分。與圖7（a）和圖7（b）相比，這兩個樣例所示的目標區域更為復雜，生成難度更大。從生成的正交網格結果來看，網格貼體性和正交性均保持較好水平，能滿足科學計算的需求。

圖8首先展示了圖7（c）和圖7（d）兩種具有較復雜邊界的圖形在不同的L的取值下，最優化問題（5）的可行解的個數。另外圖8中baseline 表示的是若沒有使用本文算法的情況下，對于40-L邊形公式（4）的可行解個數（圖7（c）和圖7（d）的目標多邊形均為40邊形）。根據1.5節的分析，本文的算法可以將時間復雜度從至少O(3N)降為大約至少O(2N-L)，從圖8 可以看出，最優化問題（5）需要枚舉的可行解個數已經大幅度減少，當L為31時分別減少了88.42%和91.16%的可行解數量，說明了算法的有效性。

正交網格是天氣預報、氣候模擬等科學應用中最為重要的網格。表1 展示的是這些領域中部分先進的正交網格生成工作。自動化網格生成是減少人力成本的關鍵，如圖9所示。

在Linux16.04，CPU 為i7-8700 環境下對比SCtoolbox 和本文方法。注意到SCtoolbox［28］也是開源且自動化生成網格方法，對于圖9（a）、（b）的 簡單圖形來說，本文方法與SCtoolbox 生成網格相同，但本文方法生成速度有明顯優勢。對于邊界復雜的圖9（c）、（d）而言：一方面SCtoolbox 在面對復雜邊界的時候效果較弱，所示樣例的網格點集中在圖形底部，沒能較好地進行網格劃分；另一方面，SCtoolbox在復雜圖形樣例中同樣需要更多的時間。Delft3D［30］是優秀的商業軟件，在實際工程任務中被廣泛使用，但其生成網格的過程中仍然需要一定的人工先驗信息。注意到大多數工作都是由C/C++或者Matlab編寫，而本文方法由Python實現，這使得本文方法具有很強的可擴展性，為科學計算領域提供更大的便利，特別是在人工智能科學計算的應用上。

3 未來的工作

3.1 更多的轉角類型

實際上，基于Gridgen-c 工具的網格生成并不一定能保證得到很好的正交網格，即最優化問題（5）的最優解并不一定能滿足科學計算的需要，其主要的原因在于Gridgen-c 工具本身，即它只有三種轉角類型，分別刻畫了180°、90°和270°的轉角，然而這并不能對所有多邊形都有效，如圖10所示。

首先，右上、右下、左下三個轉角應該選擇90°（即轉角類型1），而左上兩個角有對稱性，它們應該要有一樣的轉角類型，但是可以驗證不管它們選擇怎樣的轉角類型都不能滿足式（4）。這說明利用Gridgen-c工具不能很好解決圖10所示的區域的網格生成問題。再比如一個頂點數足夠多的正多邊形（接近圓周），其頂點的選擇也是會面臨一定的困難（圖7（b））。這意味著應該考慮更多的轉角類型，如針對圖10 的樣例，可以以45°為間隔的轉角，即180°、135°、90°、45°、225°、270°以及315°這7 類轉角，這樣對于更加復雜的目標多邊形能得到更加適合的轉角方案。不過這樣，一方面，圖3 區域A中的正交網格劃分將不再適用，需要更好的劃分方式；另一方面，由于考慮了7 類轉角，所以最優化問題（5）的時間復雜度至少為O(7N)，此時即使有圖6 的策略仍然很難讓時間復雜度降為可以承受的范圍，所以，轉角類型的設定需要有更多的研究和考慮，是未來工作的一個重點。

3.2 獲得更好的分類器

分類器性能的好壞是本文提出算法的關鍵。根據圖6 及其分析，如果能直接確定下來的轉角類型較少，那么意味著并不能大幅度地降低計算量。同時，若被直接確定的轉角類型由于分類器的性能較弱判斷失誤，那么對后續關于最優化問題（5）的求解會帶來嚴重的干擾。為了獲得更好的分類器，可以從3個角度來思考。

生成的算法 隨著機器學習和深度學習領域的不斷進步，對于序列類型數據的處理和相關分類問題的算法將越來越強，未來可以不斷地將LSTM 替換成相應的算法來完善分類器的性能。

數據集的選擇 在1.6 節中提到了本文算法使用的訓練集數據為GADM 數據集，該數據集的優點在于數量大且多邊形的邊界形狀更接近于自然區域，符合一定的物理規律。然而，由于地理類型數據邊界形狀的多樣性，GADM 的數據集里仍然存在少部分形狀比較獨特的多邊形，如3.1 節的分析，這類多邊形關于最優化問題（5）往往沒有較好的解。于是這類多邊形數據在訓練的時候便會引入較大的噪聲，影響分類器的性能。在深度學習領域，有很多經典數據集，如ImageNet2012［21］、COCO［22］、CIFAR10/100［23］等，它們都經過較好的篩選和人工標記。而由1.6 節，本文所使用的數據集僅僅進行了適量的預處理，這使得數據集的質量并不能有足夠好的保證。這意味著，網格生成領域也需要有高質量的標準數據集，而這需要較大的人力和物力。

增加訓練集多邊形的頂點數 由于最優化問題（5）時間復雜度的限制，在訓練LSTM 的時候，采用的是N=10的多邊形。然而僅僅使用頂點數為10 的多邊形進行訓練的話，不一定能滿足超大規模多邊形的網格生成需求。如圖11 所示的是利用不同頂點數量多邊形訓練得到的分類器去處理圖7（d）的結果，從結果來看，這說明了提升訓練集多邊形的頂點數對網格的生成有很大的幫助。另外，在求解最優化問題（5）的時候實際是遍歷滿足式（4）的所有解的組合，而這些解相互獨立，非常適合采用分布式的計算，隨著高性能計算領域的發展，利用多線程CPU 或者GPU 等方法能大幅度地加速式（4）的遍歷，由此能得到頂點數足夠高的多邊形數據集，從而提升分類器處理更大規模網格生成問題的能力。

圖11 訓練集中多邊形頂點數的影響Fig.11 Influence of polygon vertex number in training set

4 結語

本文首次對正交網格生成問題在數學上將SC 映射和整數規劃問題建立聯系，構建了一套較為完整的分析框架。由所提出的基本假設，利用循環神經網絡LSTM 和目標多邊形自身的信息，可以大幅度地減少網格生成所需計算量。由于對于平面上正交網格的生成問題可以歸結為多邊形區域的網格生成問題，本文所提出的方法具有一定通用性，實驗表明，在簡單圖形區域，動畫圖形區域以及地理邊界區域等多類型區域上均能自動化地產生高質量正交網格。