邢華橋，侯妙樂，王　磊，姜曉軼

(1. 中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083； 2. 北京建筑大學(xué)，

北京 100044； 3. 國(guó)家海洋信息中心，天津 300171)

Submergence Analysis Algorithm for Large Area Based on Spherical QTM

XING Huaqiao，HOU Miaole，WANG Lei，JIANG Xiaoyi

?

基于球面QTM的大范圍有源淹沒算法研究

邢華橋1，侯妙樂2，王磊1，姜曉軼3

(1. 中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083； 2. 北京建筑大學(xué)，

北京 100044； 3. 國(guó)家海洋信息中心，天津 300171)

Submergence Analysis Algorithm for Large Area Based on Spherical QTM

XING Huaqiao，HOU Miaole，WANG Lei，JIANG Xiaoyi

摘要：面對(duì)海平面上升后沿海區(qū)域大范圍淹沒模擬的應(yīng)用需求，提出了一種基于球面QTM的大范圍有源淹沒算法，首先構(gòu)建了研究區(qū)域的QTM三角格網(wǎng)集，并利用海陸分離的思想確定了初始淹沒單元；然后基于QTM的12鄰近搜索設(shè)計(jì)了面向多目標(biāo)三角格網(wǎng)的淹沒范圍擴(kuò)張方法，將初始淹沒單元遞歸擴(kuò)張得到最終的淹沒范圍；最后應(yīng)用Visual C#語言及Microsoft DirectX三維圖形接口設(shè)計(jì)開發(fā)了相應(yīng)的試驗(yàn)原型系統(tǒng)，并利用天地圖高分辨率影像和SRTM地形數(shù)據(jù)對(duì)試驗(yàn)區(qū)域進(jìn)行了海平面上升后的大范圍淹沒模擬試驗(yàn)。結(jié)果表明，該算法具有較好的實(shí)用性，可以有效地對(duì)海平面上升后的大區(qū)域范圍進(jìn)行動(dòng)態(tài)地淹沒模擬，對(duì)輔助防災(zāi)減災(zāi)具有一定意義。

關(guān)鍵詞：球面QTM；海平面上升；海陸分離；12鄰近搜索；有源淹沒

一、引言

隨著全球氣候變暖等極端天氣的頻發(fā)，海平面異常出現(xiàn)的次數(shù)明顯增多，異常值明顯增大，近30年來全球海平面上升速率有加大的趨勢(shì)[1]，而同期中國(guó)沿海海平面的上升速率高于全球的平均值[2-3]。海平面上升帶來的最為直接和嚴(yán)重的影響是沿海地區(qū)的淹沒問題。因此，準(zhǔn)確、科學(xué)地預(yù)測(cè)海平面上升后的淹沒范圍，對(duì)海洋災(zāi)害的預(yù)測(cè)預(yù)報(bào)、災(zāi)害過程的動(dòng)態(tài)模擬，以及促進(jìn)海洋災(zāi)害機(jī)理的研究、輔助防災(zāi)減災(zāi)等都具有重要意義。

近年來，國(guó)內(nèi)外研究學(xué)者對(duì)淹沒分析的研究取得了一定的進(jìn)展，淹沒分析主要分為無源淹沒和有源淹沒[4-6]，無源淹沒將所有高程值低于給定水位的區(qū)域均記為淹沒區(qū)，有源淹沒除了考慮上述情況外，還考慮到了“流通”的淹沒情形，即只有與淹沒區(qū)域相連通且高程值低于給定水位的區(qū)域才被視為淹沒區(qū)。目前的淹沒分析算法多以研究區(qū)域的DEM數(shù)據(jù)為基礎(chǔ)[7]，利用種子蔓延算法進(jìn)行淹沒區(qū)域的模擬。如張東華等提出一種基于DEM的洪水有源淹沒算法，在GIS技術(shù)的基礎(chǔ)上應(yīng)用數(shù)字高程模型(DEM)的格網(wǎng)模型進(jìn)行洪水淹沒分析[8]；丁志雄等提出基于格網(wǎng)模型的洪水淹沒分析算法，在給定洪水水位和洪量?jī)煞N條件下，分別基于三角形格網(wǎng)模型和任意多邊形格網(wǎng)模型進(jìn)行了水淹模擬[9]；聶漢江等提出了基于庫群聯(lián)合調(diào)度和DEM的區(qū)域洪水淹沒范圍模擬方法，以聯(lián)合調(diào)度模型得到的累積成災(zāi)水量為淹沒水量，以研究區(qū)域的DEM為數(shù)據(jù)基礎(chǔ)，運(yùn)用體積法對(duì)一定體量的洪水淹沒進(jìn)行了模擬[10]。但這些研究大都是針對(duì)小范圍區(qū)域，采用傳統(tǒng)的平面數(shù)據(jù)模型進(jìn)行的洪水淹沒分析，而海平面上升帶來的淹沒通常范圍較大，傳統(tǒng)的淹沒分析方法并不能滿足這種大范圍淹沒分析的需求。

針對(duì)這一問題，本文基于球面四元三角網(wǎng)(quaternary triangular mesh，QTM)的12鄰近搜索，提出一種面向大范圍區(qū)域的有源淹沒算法，以實(shí)現(xiàn)海平面上升后研究區(qū)域的淹沒動(dòng)態(tài)模擬。

二、QTM理論知識(shí)

QTM是由內(nèi)接于球體的正八面體經(jīng)多次遞歸剖分而成的球面擬合三角格網(wǎng)集[11]，具有符合球面的特性和本質(zhì)是柵格的特點(diǎn)，可以用QTM模擬地球表面；另外，QTM形狀和面積的近似相等能保證格網(wǎng)之間具有簡(jiǎn)潔的鄰近關(guān)系，便于格網(wǎng)間的鄰近查詢、分析等處理[12]。

1. QTM剖分與編碼

在QTM初始化剖分時(shí)，內(nèi)接于球體的正八面體頂點(diǎn)占據(jù)球面主要點(diǎn)(包括兩極)，而邊的投影與赤道、主子午線，以及90°、180°、270°子午線重合，如圖1所示。在對(duì)初始的QTM進(jìn)一步細(xì)化剖分時(shí)，用大弧平分法找到三角形邊的中點(diǎn)，通過連接它們將球面三角形分成4個(gè)小的三角形，以此類推進(jìn)行遞歸，對(duì)整個(gè)球面進(jìn)行近似均勻的剖分[13]，形成了全球的QTM離散格網(wǎng)集。

圖1　QTM初始剖分單元

QTM編碼不僅暗含著格網(wǎng)的空間位置，而且表示了格網(wǎng)的剖分層次，因此QTM編碼是進(jìn)行鄰近搜索及連通性分析的基礎(chǔ)。本文采用的編碼方案為固定方向編碼[12]，與其他編碼方案相比，該方案的所有格網(wǎng)都具有固定方向，更利于鄰近搜索。如圖2所示，編碼的第1位為八分體號(hào)，從第2位開始，0表示中間三角格網(wǎng)，1表示頂(底)三角格網(wǎng)，2表示左三角格網(wǎng)，3表示右三角格網(wǎng)，即可以用一個(gè)四進(jìn)制的Morton碼來標(biāo)示。因此，QTM的編碼M由八分碼M0加上Morton碼組成，即M=M0q1q2…qi。

圖2　固定方向編碼方案

2. QTM的鄰近關(guān)系

鄰近關(guān)系是GIS空間分析中必不可少的一類空間關(guān)系，是球面實(shí)體擴(kuò)張和空間索引的基礎(chǔ)[14-15]，在實(shí)際應(yīng)用中(如水淹分析、最短路徑分析等)具有重要的意義。如圖3所示，QTM三角格網(wǎng)的鄰近三角格網(wǎng)分為具有公共邊的邊鄰近(edge-adjacent)三角格網(wǎng)(編號(hào)為E)和具有公共點(diǎn)的角鄰近(vertex-adjacent)三角格網(wǎng)(編號(hào)為V)。將球面按QTM剖分后，每個(gè)球面三角格網(wǎng)都有3個(gè)邊鄰近三角格網(wǎng)，位于南、北極的8個(gè)頂角三角形和赤道上經(jīng)度為0、90°、180°(-180°)、-90°的16個(gè)頂角三角形均具有7個(gè)角鄰近三角形(如圖3(a)所示)，其他位置的非頂角三角形均具有9個(gè)角鄰近三角形(如圖3(b)所示)。根據(jù)QTM的鄰近關(guān)系，其鄰近搜索分為三鄰近搜索和12鄰近搜索，由于海水淹沒不具有方向性，即海水可能向所有方向蔓延，因此，利用12鄰近搜索更符合實(shí)際情況。本文主要考慮QTM格網(wǎng)的12鄰近搜索[16]。

圖3　QTM的鄰近三角格網(wǎng)

三、基于QTM 12鄰近搜索的有源淹沒算法

1. 總體思路

首先構(gòu)建研究區(qū)域的一定層次的QTM格網(wǎng)集，基于海陸分層的思想確定初始淹沒單元，并以此為基礎(chǔ)進(jìn)行12鄰近搜索獲取其鄰近格網(wǎng)，如果鄰近格網(wǎng)的高程低于給定的上升水位，即HQTM

圖4　基于QTM 12鄰近搜索的有源淹沒算法流程

2. 基于海陸分離思想的初始淹沒單元判定

確定初始淹沒區(qū)域是海平面上升分析與模擬的首要問題，海平面上升后的淹沒過程是從陸地與海洋相接處開始的。因此，本文利用海陸分離的思想，將研究區(qū)域的QTM按照3個(gè)頂點(diǎn)的高程劃分為陸地和海洋，并確定陸地邊緣格網(wǎng)。具體方法為：對(duì)于每個(gè)格網(wǎng)，如果其3個(gè)頂點(diǎn)均位于海洋(頂點(diǎn)的高程值均為-9999 m，如圖5(a)所示)則為海洋格網(wǎng)；如果其3個(gè)頂點(diǎn)均位于陸地(頂點(diǎn)的高程值均大于-9999 m，如圖5(b)所示)則為陸地格網(wǎng)；如果其3個(gè)頂點(diǎn)并不完全位于陸地或海洋，那么將該格網(wǎng)視為陸地邊緣格網(wǎng)(如圖5(c)所示)，所有陸地邊緣格網(wǎng)構(gòu)成的集合即為初始淹沒單元。

圖5　QTM三角格網(wǎng)陸地與海洋的判別

3. 面向多目標(biāo)三角格網(wǎng)的淹沒范圍擴(kuò)張

由淹沒單元通過鄰近搜索擴(kuò)張得到其鄰近格網(wǎng)時(shí)，淹沒單元往往不是單獨(dú)的一個(gè)QTM格網(wǎng)，而是多個(gè)連續(xù)的格網(wǎng)，這就對(duì)鄰近搜索提出了特殊的需求，即需要對(duì)多個(gè)三角格網(wǎng)同時(shí)進(jìn)行擴(kuò)張，而面向單目標(biāo)三角格網(wǎng)的擴(kuò)張方法會(huì)造成格網(wǎng)的重復(fù)擴(kuò)張，從而增加了系統(tǒng)的計(jì)算量。本文在淹沒范圍的擴(kuò)展過程中，剔除擴(kuò)張得到的鄰近格網(wǎng)中的重復(fù)格網(wǎng)，同時(shí)為進(jìn)一步減小工作量，每一次擴(kuò)張只將上一次擴(kuò)張得到的被淹沒格網(wǎng)向外擴(kuò)張，而不再處理已經(jīng)判斷過的格網(wǎng)。面向多目標(biāo)三角格網(wǎng)的淹沒范圍擴(kuò)張過程如圖6所示。該擴(kuò)張方法的過程具體如下：

1) 利用12鄰近搜索算法依次計(jì)算初始淹沒單元S中所有格網(wǎng)的鄰近格網(wǎng)，并從中剔除重復(fù)格網(wǎng)，將不重復(fù)的鄰近格網(wǎng)集合記作G。

2) 利用雙線性插值算法計(jì)算格網(wǎng)集G中每個(gè)格網(wǎng)中心點(diǎn)的高程。

3) 將格網(wǎng)集G中每個(gè)格網(wǎng)的高程與當(dāng)前水位數(shù)據(jù)進(jìn)行對(duì)比，高程低于當(dāng)前水位的格網(wǎng)為被淹沒的格網(wǎng)，記作G′。

4) 以G′代替步驟1)中的初始淹沒單元S，執(zhí)行步驟1)—步驟3)，直至G中所有格網(wǎng)的高程均高于當(dāng)前水位(即G′為空)或整個(gè)試驗(yàn)區(qū)的格網(wǎng)搜索完畢。

圖6　面向多目標(biāo)三角格網(wǎng)的淹沒范圍擴(kuò)張過程

四、試驗(yàn)與分析

為了對(duì)算法進(jìn)行驗(yàn)證，本文以C# 3.0為開發(fā)語言，DirectX為三維開發(fā)包，從底層開發(fā)了試驗(yàn)原型系統(tǒng)。試驗(yàn)硬件環(huán)境為Pentium(R) Dual-core CPU T4200 @2.00 GHz，2.00 GB內(nèi)存，采用SRTM 地形數(shù)據(jù)作為DEM數(shù)據(jù)源，天地圖高分辨率影像為影像數(shù)據(jù)源，在給定淹沒水位的情形下對(duì)海平面上升后研究區(qū)域進(jìn)行動(dòng)態(tài)淹沒模擬。試驗(yàn)中首先構(gòu)建了一定層次的QTM三角格網(wǎng)集作為研究區(qū)域的模擬圖層，然后利用海陸分離思想提取試驗(yàn)區(qū)域的初始淹沒單元，如圖7所示。

圖7　試驗(yàn)區(qū)域的初始淹沒單元

進(jìn)而利用本文提出的算法在給定淹沒水位下獲取所有被淹沒的格網(wǎng)存入淹沒格網(wǎng)集，并將淹沒格網(wǎng)集中所有格網(wǎng)進(jìn)行著色和渲染，然后與無源淹沒算法的淹沒效果進(jìn)行了對(duì)比。圖8為淹沒水位為5 m的試驗(yàn)區(qū)域淹沒范圍的模擬結(jié)果，圖9為利用無緣淹沒算法得到的試驗(yàn)區(qū)域淹沒范圍。通過淹沒模擬對(duì)比圖可以看到，利用本文算法得到的淹沒區(qū)域具有連通性，更加符合淹沒的實(shí)際情況。

為了進(jìn)一步比較本文算法與無源淹沒算法的區(qū)別，并驗(yàn)證本文算法在大區(qū)域范圍進(jìn)行淹沒模擬的有效性，在相同的試驗(yàn)硬件環(huán)境下，應(yīng)用不同的淹沒水位分別在不同面積的大范圍試驗(yàn)區(qū)域進(jìn)行有源淹沒和無源淹沒模擬，并比較了不同淹沒算法的淹沒面積，見表1。試驗(yàn)結(jié)果表明，隨著淹沒水位的升高，淹沒范圍逐漸增大；在利用相同淹沒水位對(duì)同一的試驗(yàn)區(qū)域進(jìn)行模擬時(shí)，本文算法得到的淹沒范圍相對(duì)于無源淹沒算法較小，這是由于本文算法考慮到了淹沒區(qū)域的連通性，只有高程低于給定水位且與初始淹沒單元相連通的區(qū)域才被視為最終的淹沒區(qū)域。此外，與傳統(tǒng)的有源淹沒算法相比，本文算法所應(yīng)用的試驗(yàn)區(qū)域達(dá)到了數(shù)萬平方千米的級(jí)別，可以有效地應(yīng)用于更大范圍的淹沒分析中， 并且能將淹沒結(jié)果動(dòng)態(tài)地在三維球面上進(jìn)行模擬和展示。

圖8　本文算法得到的淹沒模擬圖

圖9　無源淹沒算法得到的淹沒模擬圖

km2

五、結(jié)論與展望

本文針對(duì)海平面上升后的沿海區(qū)域淹沒模擬問題，提出了基于球面QTM 12鄰近搜索的大范圍有源淹沒算法，并設(shè)計(jì)開發(fā)了相應(yīng)的原型試驗(yàn)系統(tǒng)。試驗(yàn)證明該方法可以用于球面大范圍海平面上升后的水淹動(dòng)態(tài)模擬，能夠?yàn)榉罏?zāi)減災(zāi)提供決策依據(jù)。但該算法也有一定的局限性，海平面上升是個(gè)漸變的過程，其年變化范圍在毫米量級(jí)，而本文假定的海平面上升水位實(shí)際是疊加了平均潮位和極值潮位的最終上升水位。如何根據(jù)真實(shí)的海平面上升高度模擬淹沒范圍，并對(duì)海平面上升后給淹沒區(qū)域帶來的社會(huì)、經(jīng)濟(jì)影響等進(jìn)行分析將是下一步的研究重點(diǎn)。

參考文獻(xiàn)：

[1]李響，劉克修， 董軍興，等. 中國(guó)沿海海平面上升風(fēng)險(xiǎn)區(qū)劃方法研究[J].海洋預(yù)報(bào)， 2014，31(2)：41-49.

[2]周子鑫.我國(guó)海平面上升研究進(jìn)展與前瞻[J]. 海洋地質(zhì)動(dòng)態(tài)，2008， 24(10)：14-18.

[3]王慧，劉克修，范文靜，等. 2012年中國(guó)沿海海平面上升顯著成因分析[J].海洋學(xué)報(bào)， 2014，36(5)：8-17.

[4]劉仁義，劉南. 基于GIS的復(fù)雜地形洪水淹沒區(qū)計(jì)算方法[J].地理學(xué)報(bào)，2001，56(1)：1-6.

[5]劉小生，陳英俊，黃玉生.基于GIS技術(shù)的洪水淹沒區(qū)確定[J]. 測(cè)繪科學(xué)，2007，32(5)：136-137.

[6]王建軍，吳北平，朱萬里，等.基于ArcGIS的風(fēng)暴潮災(zāi)害疏散系統(tǒng)的研究與開發(fā)[J].測(cè)繪通報(bào)，2013(9)：91-94.[15]國(guó)家測(cè)繪局．GB/T13990—921∶5000、1∶10000地形圖航空攝影測(cè)量?jī)?nèi)業(yè)規(guī)范[S].北京：中國(guó)標(biāo)準(zhǔn)出版社，1992.

[7]韓麗蓉，張濤.基于RS/GIS的大南川水庫潰壩災(zāi)害模擬方法研究[J].測(cè)繪通報(bào)，2014(1)：97-99.

[8]張東華，劉榮，張?jiān)佇?，?一種基于DEM的洪水有源淹沒算法的設(shè)計(jì)與實(shí)現(xiàn)[J].東華理工大學(xué)學(xué)報(bào):自然科學(xué)版，2009，32(3)：181-184.

[9]丁志雄，李紀(jì)人，李琳.基于GIS 格網(wǎng)模型的洪水淹沒分析方法[J].水利學(xué)報(bào)，2004(6)：56-60.

[10]聶漢江，劉彬，何立新，等. 基于庫群聯(lián)合調(diào)度和DEM的區(qū)域洪水淹沒范圍模擬[J].水電能源科學(xué)，2014，32(5)：43-46.

[11]DUTTON G. Encoding and Handling Geospatial Data with Hierarchical Triangular Meshes [C]∥Proceeding of 7th International Symposium on Spatial Data Handling. Netherlands: Delft University of Technology, 1996. 34-43.

[12]SAHR K, WHITE D，KIMERLING A. Geodesic Discrete Global Grid Systems[J]. Cartography & Geographical Information Science, 2003,30(2): 121-134.

[13]LEE M，SAMET H. Navigating through Triangle Meshes Implemented as Linear Quadtree[J]. ACM Transactions on Graphics, 2000,19(2): 79-121.

[14]陳軍，趙仁亮.GIS 空間關(guān)系的基本問題與研究進(jìn)展[J]. 測(cè)繪學(xué)報(bào)，1999，28(2)：96-102.

[15]侯妙樂，邢華橋，趙學(xué)勝，等.球面四元三角網(wǎng)的復(fù)雜拓?fù)潢P(guān)系計(jì)算[J].武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2012，37(4)：468-471.

[16]趙學(xué)勝，侯妙樂，白建軍.全球離散格網(wǎng)的空間數(shù)字建模[M].北京：測(cè)繪出版社，2007：40-46.引文格式： 顏金彪，胡最，禹信. CORS與靜態(tài)相對(duì)定位技術(shù)下的坐標(biāo)轉(zhuǎn)換研究[J].測(cè)繪通報(bào)，2015(12)：54-56.DOI:10.13474/j.cnki.11-2246.2015.377

作者簡(jiǎn)介：邢華橋(1988—)，男，博士生，研究方向?yàn)槿螂x散格網(wǎng)、服務(wù)計(jì)算。E-mail：xinghuaqiao@126.com

基金項(xiàng)目：國(guó)家自然科學(xué)基金(41171304；40701152；41171306)

收稿日期：2015-02-26； 修回日期： 2015-10-27

中圖分類號(hào)：P208

文獻(xiàn)標(biāo)識(shí)碼：B

文章編號(hào)：0494-0911(2015)12-0046-04