基于循環(huán)神經(jīng)網(wǎng)絡(luò)的正交網(wǎng)格的自動化生成算法

黃中展，徐世明

（清華大學(xué)地球系統(tǒng)科學(xué)系，北京 100084）

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

0 引言

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

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

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

圖2 多連通區(qū)域的三角網(wǎng)格［6］Fig.2 Triangular grid of multiconnected region

對于平面上的目標(biāo)區(qū)域，三角網(wǎng)格和四邊形網(wǎng)格的算法都均可以對目標(biāo)區(qū)域進(jìn)行離散化。三角形網(wǎng)格主要采用的是Delaunay 三角剖分法［7］、四叉樹法［8］以及波前法［9］，其中Delaunay 法使用最為普遍，這是由于它在數(shù)學(xué)上十分成熟，具有較好的網(wǎng)格質(zhì)量且收斂性也能較好地保證；而四邊形網(wǎng)格的生成主要有格柵法［10］、拓?fù)浞纸夥ǎ?1］、節(jié)點連接法［12］等。在天氣預(yù)報、氣候模擬等科學(xué)領(lǐng)域中，正交網(wǎng)格應(yīng)用最為廣泛，本文主要考慮的是正交網(wǎng)格的自動化生成。

對于單連通區(qū)域的正交網(wǎng)格生成來說，利用Thompson 微分方程法［13］最為普遍，然而這類方法產(chǎn)生的網(wǎng)格在邊界處正交性保持較差，收斂性也不能保證，且具有較多的需要調(diào)節(jié)的參數(shù)，這使得當(dāng)遇到較復(fù)雜的目標(biāo)區(qū)域的時候，需要有較多的人工調(diào)節(jié)，不便于科學(xué)計算的實際應(yīng)用?？焖偾易詣踊纳伤惴ㄔ诟鞣N科學(xué)計算領(lǐng)域當(dāng)中十分必要。

1 正交網(wǎng)格的自動化生成

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

1.1 Schwarz-Christoffel共形映射

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

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

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

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

1.2 Gridgen-c網(wǎng)格生成工具

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

圖4 Gridgen-c網(wǎng)格生成工具Fig.4 Grid generation tool Gridgen-c

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

1.3 正交網(wǎng)格生成對應(yīng)的最優(yōu)化問題

本節(jié)首先對正交網(wǎng)格生成問題建模，即將利用Gridgen-c所需的條件將正交網(wǎng)格的生成問題轉(zhuǎn)換為一個帶線性限制條件的整數(shù)規(guī)劃問題。

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

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

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

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

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

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

均勻性 生成的網(wǎng)格應(yīng)保持足夠均勻，否則容易產(chǎn)生極端的網(wǎng)格點分布，盡管這些網(wǎng)格點所構(gòu)成的網(wǎng)格的正交性可能可以保持足夠好，但由于網(wǎng)格不均勻所導(dǎo)致的尺度差異，一定程度上不利于其實際的科學(xué)計算，可能會產(chǎn)生較大的數(shù)值誤差等問題。

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

結(jié)合式（4）和上述正交性、均勻性和覆蓋性的考慮，可以構(gòu)造如下最優(yōu)化問題：

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

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

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

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

針對最優(yōu)化問題（5）時間復(fù)雜度過高的問題，提出了利用循環(huán)神經(jīng)網(wǎng)絡(luò)構(gòu)建一個關(guān)于轉(zhuǎn)角選擇的分類器來降低時間復(fù)雜度，從而獲得較優(yōu)近似解的辦法。

1.4 循環(huán)神經(jīng)網(wǎng)絡(luò)

為了可以利用循環(huán)神經(jīng)網(wǎng)絡(luò)來降低最優(yōu)化問題（5）求解的時間復(fù)雜度，此處提出一個基本假設(shè)：

基本假設(shè) 一個頂點的轉(zhuǎn)角類型取決于其鄰近點的轉(zhuǎn)角類型，且越接近的頂點影響越大。

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

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

由此，對于每一頂點vi，都可以構(gòu)造如下特征

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

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

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

圖5 LSTM結(jié)構(gòu)Fig.5 Structure of LSTM

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

1.5 降低時間復(fù)雜度

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

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

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

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

1.6 訓(xùn)練數(shù)據(jù)

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

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

2）保證考慮的多邊形是單連通區(qū)域。在GADM數(shù)據(jù)集中包含很多群島等地理邊界，這些邊界本身就構(gòu)成了多連通的區(qū)域，這些數(shù)據(jù)應(yīng)該被預(yù)先剔除。

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

2 實驗與模型分析

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

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

圖7 正交網(wǎng)格生成樣例Fig.7 Examples of orthogonal grid generation

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

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

正交網(wǎng)格是天氣預(yù)報、氣候模擬等科學(xué)應(yīng)用中最為重要的網(wǎng)格。表1 展示的是這些領(lǐng)域中部分先進(jìn)的正交網(wǎng)格生成工作。自動化網(wǎng)格生成是減少人力成本的關(guān)鍵，如圖9所示。

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

3 未來的工作

3.1 更多的轉(zhuǎn)角類型

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

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

3.2 獲得更好的分類器

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

生成的算法 隨著機(jī)器學(xué)習(xí)和深度學(xué)習(xí)領(lǐng)域的不斷進(jìn)步，對于序列類型數(shù)據(jù)的處理和相關(guān)分類問題的算法將越來越強(qiáng)，未來可以不斷地將LSTM 替換成相應(yīng)的算法來完善分類器的性能。

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

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

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

4 結(jié)語

本文首次對正交網(wǎng)格生成問題在數(shù)學(xué)上將SC 映射和整數(shù)規(guī)劃問題建立聯(lián)系，構(gòu)建了一套較為完整的分析框架。由所提出的基本假設(shè)，利用循環(huán)神經(jīng)網(wǎng)絡(luò)LSTM 和目標(biāo)多邊形自身的信息，可以大幅度地減少網(wǎng)格生成所需計算量。由于對于平面上正交網(wǎng)格的生成問題可以歸結(jié)為多邊形區(qū)域的網(wǎng)格生成問題，本文所提出的方法具有一定通用性，實驗表明，在簡單圖形區(qū)域，動畫圖形區(qū)域以及地理邊界區(qū)域等多類型區(qū)域上均能自動化地產(chǎn)生高質(zhì)量正交網(wǎng)格。