黃中展，徐世明

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

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

0 引言

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

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

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

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

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

對(duì)于單連通區(qū)域的正交網(wǎng)格生成來(lái)說(shuō)，利用Thompson 微分方程法［13］最為普遍，然而這類方法產(chǎn)生的網(wǎng)格在邊界處正交性保持較差，收斂性也不能保證，且具有較多的需要調(diào)節(jié)的參數(shù)，這使得當(dāng)遇到較復(fù)雜的目標(biāo)區(qū)域的時(shí)候，需要有較多的人工調(diào)節(jié)，不便于科學(xué)計(jì)算的實(shí)際應(yīng)用。快速且自動(dòng)化的生成算法在各種科學(xué)計(jì)算領(lǐng)域當(dāng)中十分必要。

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

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

1.1 Schwarz-Christoffel共形映射

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

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

共形映射具有一一對(duì)應(yīng)且保持角度這兩個(gè)良好性質(zhì)，這些性質(zhì)可以用于正交網(wǎng)格的生成。如圖3 所示，若區(qū)域B為目標(biāo)區(qū)域，即需要在區(qū)域B當(dāng)中生成正交網(wǎng)格?？梢允紫壤肧C 映射建立區(qū)域B到區(qū)域H的共形映射g。若可以建立一個(gè)簡(jiǎn)單的、容易構(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)格的理論相對(duì)完備，但是在實(shí)際網(wǎng)格生成應(yīng)用以及算法上仍有很多不足。在眾多基于共形映射的網(wǎng)格生成工具中，Gridgen-c 工具是基于C 語(yǔ)言的正交網(wǎng)格生成工具，具有快遞且穩(wěn)定的網(wǎng)格生成能力。其不足主要在與需要人為給定目標(biāo)多邊形的轉(zhuǎn)角類型，利用這些人為標(biāo)定的轉(zhuǎn)角類型，從給定的目標(biāo)區(qū)域B中人為地構(gòu)造出圖3中的區(qū)域A。具體來(lái)說(shuō)，如圖4所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實(shí)驗(yàn)與模型分析

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

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

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

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

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

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

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

3 未來(lái)的工作

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

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

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

3.2 獲得更好的分類器

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

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

數(shù)據(jù)集的選擇 在1.6 節(jié)中提到了本文算法使用的訓(xùn)練集數(shù)據(jù)為GADM 數(shù)據(jù)集，該數(shù)據(jù)集的優(yōu)點(diǎn)在于數(shù)量大且多邊形的邊界形狀更接近于自然區(qū)域，符合一定的物理規(guī)律。然而，由于地理類型數(shù)據(jù)邊界形狀的多樣性，GADM 的數(shù)據(jù)集里仍然存在少部分形狀比較獨(dú)特的多邊形，如3.1 節(jié)的分析，這類多邊形關(guān)于最優(yōu)化問(wèn)題（5）往往沒(méi)有較好的解。于是這類多邊形數(shù)據(jù)在訓(xùn)練的時(shí)候便會(huì)引入較大的噪聲，影響分類器的性能。在深度學(xué)習(xí)領(lǐng)域，有很多經(jīng)典數(shù)據(jù)集，如ImageNet2012［21］、COCO［22］、CIFAR10/100［23］等，它們都經(jīng)過(guò)較好的篩選和人工標(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)練集多邊形的頂點(diǎn)數(shù) 由于最優(yōu)化問(wèn)題（5）時(shí)間復(fù)雜度的限制，在訓(xùn)練LSTM 的時(shí)候，采用的是N=10的多邊形。然而僅僅使用頂點(diǎn)數(shù)為10 的多邊形進(jìn)行訓(xùn)練的話，不一定能滿足超大規(guī)模多邊形的網(wǎng)格生成需求。如圖11 所示的是利用不同頂點(diǎn)數(shù)量多邊形訓(xùn)練得到的分類器去處理圖7（d）的結(jié)果，從結(jié)果來(lái)看，這說(shuō)明了提升訓(xùn)練集多邊形的頂點(diǎn)數(shù)對(duì)網(wǎng)格的生成有很大的幫助。另外，在求解最優(yōu)化問(wèn)題（5）的時(shí)候?qū)嶋H是遍歷滿足式（4）的所有解的組合，而這些解相互獨(dú)立，非常適合采用分布式的計(jì)算，隨著高性能計(jì)算領(lǐng)域的發(fā)展，利用多線程CPU 或者GPU 等方法能大幅度地加速式（4）的遍歷，由此能得到頂點(diǎn)數(shù)足夠高的多邊形數(shù)據(jù)集，從而提升分類器處理更大規(guī)模網(wǎng)格生成問(wèn)題的能力。

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

4 結(jié)語(yǔ)

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