潘竹爐
地球表面上71%的區域是海洋,海洋把世界聯系到一起,然而到目前為止,人類已探測的海底只有5%,還有95%的海底是未知的[1]。海洋不僅是生命的源泉,提供了人類賴以生存的物質和空間,還蘊含豐富的生物、化學、礦產等資源。隨著人類活動空間的擴展,航海早已不局限于海面上,潛艇、水下機器人等航海新設備急切需要水下三維空間地形圖。二維電子海圖只能用水深和等深線來描述水下的地形地貌,有表述不直觀、信息不全等問題,不利于海底地形分析。利用傳統的二維矢量電子海圖提供的水深數據、等深線和海岸線等信息構建水下DEM,將使水下空間可觀測維數大大增加。三維海底地形可視化系統操作更簡便、表達更直觀、逼真性更強、信息量更大,更有利于人的接受、感知和理解。
數字高程模型(DEM)一般來說是將離散、不規則的原始采樣高程值使用某些數學手段(例如過采樣、擬合、插值等)重建出相似曲面來逼近原始的地表。常見的DEM模型根據數據結構的不同主要分為規則格網DEM(如圖1(a))和不規則三角網DEM(如圖1(b))。規則格網數據結構冗余較多,但是拓撲關系較好,數據存儲簡單;不規則三角網有較好的數學結構,數據冗余少,但拓撲關系較差,數據存儲復雜[2]。

圖1 常見數字高程模型
數字高程模型的平滑內插是基于海底地形的過渡平滑趨勢構建三維海底地形的,通過二維海圖提供的離散水深點的水深,估計擬合離散點之間的水深,從而構建出完整的水下地形。由海底的過度平滑趨勢可知,距離越近的點,相關度越高。由圖2所示,按內插點的區域可分為分塊內插、逐點內插和整體內插。但是考慮到數據量的問題,一般采用逐點內插,目前,較簡單的算法有特定區域的圓、TIN三角網插值、基于VoronoiCells的自然臨近點插值[3](例如Sibson插值、Laplace插值)、基于數據融合的插值算法[4]。

圖2 DEM平滑內插分類
分形理論隸屬于非線性科學,在生物、地球科學、材料物理等領域都有著很廣泛的應用[5]。分形理論萌芽于19世紀60年代,于20世紀七、八十年代正式演變為一門獨立的學科,分形插值是分形理論的一個重要應用。分形插值主要應用了分形的自相似原理,利用迭代算法由原始特征插值出局部細小特征。分形能表達出自然界中許多非線性現象,使用分形插值理論構建地表或許最能逼近自然界真實地形[6]。
分形理論(Fractal Theory)是一門以不規則幾何形態為研究對象的新理論、新學科。分形對象在自然界是普遍存在的,粗糙不平的地表、彎彎曲曲的海岸線、綿延不絕的山脈、裊裊的炊煙、枝繁葉茂的數木、漫天飛舞的雪花等都是分形。
Mandelbrot因為提出了“英國的海岸線有多長”的問題有了分形思想萌芽的故事幾乎成為佳話被廣為流傳。大量經典的分形是自相似的,這種分形的局部和整體具有一定比例的幾何相似性,如Sierpinski墊片(圖3(a)),Koch曲線(圖3(b))。自相似集的生成元中,每一個壓縮映射都把它變換為一個小拷貝,它在不同方向上的壓縮比是相同的。多次迭代過程最終輸出圖形與輸入變得無關,真正決定最終輸出的是反饋機器中的數學變換,Barns?ley等把這種數學變換稱為迭代函數系(Iterated Functional System),簡寫為IFS。

圖3 常見的分形圖形
分形插值由Barnsley在IFS的基礎上歸納推導出來[7]。它是分形理論的一個重要應用。自Euclid創立幾何以來,人們提出了各種插值算法來描述幾何體。然而對于現實世界的大量離散數據,除了直到20世紀60年代提出的樣條函數插值[8],幾乎沒有任何方法能完美克服幾何體的低階全局光滑問題。可見數學插值更偏愛“光滑”,“平滑”等字眼。分形插值函數[9]利用大自然中廣泛存在的精細自相似結構特征來擬合滾動性強的曲線,為解決幾何體的低階全局光滑問題開辟了新的領域。
3.2.1 三維IFS分形插值算法推導
兩個維度的數據集合 I=[a,b],J=[c,d];設平面點集 D=I×J={(x,y)∶x∈I,y∈J} ,以 △x,△y 為步長,將D按M×N格網化:

設(xm,yn){m=0,1,…,M;n=0,1,…,N}點處的高程數據為zm,n,我們的目的是構造二元分形插值函數 f∶D→R ,使

設滿足上式的X方向和Y方向仿射變換為

且滿足邊界條件:
由式(3)和(4)可解得

則可得平面點集x和y方向上得仿射變換為

設高程z方向的仿射變換為

式中hm,n為高程z方向上的壓縮倍率,它是決定分形插值后曲面的粗糙起伏程度的自由參數,需滿足 0≤hm,n<1(否則IFS不收斂),一般稱它為垂直比例因子。
而且可得邊界條件:

聯立式(7)和(8)可得

令 Wm,n(x,y,z)=(?m(x),ψn(y),Fm,n(x,y,z));(m=1,2,…,M;n=1,2,…,N),那么就定義了一個迭代函數系統。 Fm,n(x,y,z)為分形插值函數fm,n(x,y)的隱函數。
3.2.2 分形插值曲面的計算步驟
根據上面推導的公式,分形插值曲面的算法步驟如下:
1)輸入待插值的格網化后的點云(xm,yn,zm,n);m=0,1,2,…,M;n=0,1,2,…,N ;
2)輸入自由參數 hm,n;
3)根據式(9)計算參數 gm,n,em,n,fm,n,km,n;
4)根 據 式(6)和(7)計 算 :?m(xi),ψn(yj),Fm,n(xi,yj,zi,j) ,m=1,2,…,M ;n=1,2,…,N ;i=0,1,…,M ;j=0,1,…N ;
5)經過步驟4)原來的 (M+1)×(N+1)格網變成((M-1)2+1)×((N-1)2+1)格網。數據量大大增加,這就是一次分形插值曲面。若不夠密集,可以繼續進行二次迭代。
為了驗證測試算法的插值效果,選取一個上海長興島附近的二維電子海圖作為數據源,如圖4所示,圖中左邊陸地是橫沙島,右邊是長興島,上面的黑點就是離散點的水深位置分布。

圖4 離散水深點在谷歌地球上的分布
為了方便構建DEM,在陸地部分加入SRTM[10]高程數據,然后將補充數據后的平面離散點集使用平滑內插算法—移動擬合法生成的DEM如圖5所示。

圖5 平滑內插算法輸出的DEM
原始二維海圖離散的水深數據共有1790組(以一個高程點數據為一組),以平均點距的1.62倍格網化成29×24格網化,使用Matlab實現分形插值曲面算法得到785×530的格網數據,如圖6所示,選取的自由參數為0.3,此時由于自由參數比較大,所以地形起伏比較大,但是基本上保持了實際的陸地海洋邊界,如果自由參數選的足夠小,則分形后的地形變得很光滑,也就退化為傳統插值算法。
由圖5和圖6可以看出,分形插值前后的主要變化是
1)離散數據集比較小的情況下,數據量大大增加,高程點變得更加密集;
2)分形插值后局部和整體具有一定比例的幾何相似性,這就是分形的自相似性。
選取自由參數為0.1的分形插值后的三維高程數據,采用DEM可視化技術,導入MultiGen Creator[11]地形建模工具生成的視圖如圖7所示。
三維海底地形圖在在軍事仿真、水下作業、游戲開發和虛擬現實等方面都有十分廣闊的應用前景。本文綜合利用離散的水深信息構建數字高程模型,推導了數字高程模型的三維分形內插算法,并通過Matlab建模比較了傳統平滑內插算法和分形插值算法,為三維水下地形圖提供了一種行之有效的實現方法。
[1]http://www.noaa.gov/ocean.html[EB/OL].
[2]沈海峰,李慶陽.電子海圖高程數據的離散化研究[C]//廣州:第四屆廣東海事高級論壇論文集,2012:94-97.
[3]申浩,田峰敏,趙玉新.基于VoronoiCells插值的三維海底地形圖生成方法[J].系統仿真學報,2006,18(2):444-446,450.
[4]趙凱,胡大斌,肖建波.基于數據融合的DEM插值與可視化[J].計算機仿真,2012,29(7):144-146,213.
[5]朱華,姬翠翠.分形理論及其應用[M].北京:科學出版社,2011.1,27-47.
[6]邱秋香.三維海底地形仿真技術的研究與實現[D].哈爾濱:哈爾濱工程大學,2011:9.
[7]盧福丹,陳元枝,蔡續.三維海底地形仿真的研究[J].計算機應用與軟件,2014,31(1):77-80.
[8]陳颙,陳凌.分形幾何學[M].北京:地震出版社,2005:5-6,35-44.
[9]孫洪泉.分形幾何與分形插值[M].北京:科學出版社,2011:99-219.
[10]http://en.wikipedia.org/wiki/SRTM[EB/OL].
[11]王乘,周均清,李利軍.Creator可視化仿真建模技術[M].武漢:華中科技大學出版社,2005:261,356.