吉瑋，朱文菊，李超，畢建濤

(1．國家海洋局數字海洋科學技術重點實驗室，天津300171;2．中國科學院 遙感與數字地球研究所 數字地球重點實驗室，北京100094;3．北京中科數遙信息技術有限公司，北京100101;4．山東省建設發(fā)展研究院，山東濟南250001)

海上交通、海洋調查、海洋資源開發(fā)、海洋工程建設、海洋疆界勘定和海洋環(huán)境保護等，都需要了解海洋水下地形這一重要要素。但由于海上多變的自然環(huán)境造成的較高作業(yè)難度以及多波束測深儀造價昂貴等原因，中國海洋地形測量目前仍以單波束測深儀獲取點源數據為主。這就需要對點源數據進行插值或估計而形成連續(xù)變化的曲面來對水下地形的變化進行表達。

地形高程插值方法最常見且使用較多的是克里格 (Kriging)插值法［1-2］，但是一般的Kriging方法為了體現某種趨勢而采用加權平均估計，將實際高程值較大的地方消減，實際高程值較小的地方增大，對數值造成一定的平滑，結果使高程局部結構特征的表示與實際地形有較大差距。Mandelbrot在20世紀60年代提出的分形理論作為處理地形非線性和復雜性研究的有力工具，能夠對系統(tǒng)演變過程中的復雜性、不規(guī)則性和不均勻性進行度量，且已在海岸地貌、流水地貌和喀斯特地貌等方面的研究中得到了廣泛應用［3-6］，其不僅能夠用來描述分形的復雜特征，還能特征描述分形本身的幾何支撐度量。對于高程、水流等諸多非均勻的分形現象，可以使用多重分形測度或維數的連續(xù)譜來表示［7］。目前，多重分形方法已應用到地球科學的諸多領域［8-13］，但應用于海底地形插值方面的研究還不多見。本研究中，以海州灣作為研究區(qū)域，通過數字化海圖水深數據，將多重分形理論引入到Kriging這一常見的方法中構建海底地形，并將插值結果與普通Kriging法進行比較，旨在對海洋水深這一非均勻分形現象進行更精確地預測。

1 研究區(qū)域概況與分析方法

1.1 研究區(qū)域概況

海州灣位于蘇魯隆起與蘇北之南黃海拗陷的過渡地帶，北起山東日照市嵐山區(qū)的佛手咀 (35°05'55″N，119°21'53″E)，南至江蘇連云港市高公島(34°45'25″N，119°29'45″E)，是一個瀕臨黃海的大喇叭口形開敞海灣，灣口寬42 km，岸線長86．81 km，海灣面積為 876．39 km2［14-15］。由于供沙條件、水動力條件和岸坡形態(tài)的不同，海州灣地貌特征和沖淤動態(tài)各異。自20世紀30年代建港以來，該地區(qū)各種海洋建設工程持續(xù)進行，各類近海岸海洋活動頻繁，特別是1994年攔海西大堤的建成進一步影響了海底的沉積環(huán)境［16-17］。海州灣作為江蘇沿海開發(fā)的龍頭型增長極以及連云港實施海灣型經濟的重要載體，對這一區(qū)域海底地形的構建能對相關基礎設施建設、旅游資源的利用與開發(fā)，以及漁業(yè)養(yǎng)殖等海洋經濟產業(yè)的規(guī)劃提供較為準確的科學依據。

1.2 研究數據

本研究中，海底地形構建使用的水深數據來源于2000年12月出版的日照港至灌河口海圖，比例尺為1∶120 000，采用墨卡托投影，坐標系為1954北京坐標系，高程基準為1985國家高程基準，基本等高距為40 m。海圖數據的采集單位為中國人民解放軍海軍司令部航海保證部，圖中水深數據由海軍海洋測量船使用測深儀測得。將海圖掃描之后利用ArcGIS軟件完成投影轉換與配準的預處理工作，再在ArcGIS中采用人工矢量化的方式采集水深點數據，在約4265 km2研究區(qū)范圍內共采集909個水深點，水深為0．05～37．00 m，平均水深為11．85 m。總體樣本水深有一定差異，符合本研究需求，從中隨機抽取709個作為插值樣點用于插值研究，其余200個水深點用作檢驗 (圖1)。

1.3 分析方法

1．3．1 Kriging插值法 Kriging插值法又稱空間局部插值法，是以變異函數理論和結構分析為基礎，在一定的區(qū)域內對區(qū)域化變量進行無偏最優(yōu)估計的一種方法，是地質學研究的主要方法之一。Kriging法是根據待插值點與臨近實測點的空間位置，對待插值點的高程值進行線性無偏最優(yōu)估計，通過生成一個關于高程的Kriging插值圖來表達研究區(qū)域的原始地形。本研究中，采用普通Kriging方法進行插值，具體的原理與插值步驟參見文獻 ［2］，具體計算過程在ArcGIS軟件中完成。

1．3．2 多重分形插值法 多重分形 (Multifractal)是指用多個維數來描述非均勻的復雜幾何體，以此來全面刻畫其特征。在地形變化復雜的地區(qū)或者相對較大的區(qū)域，高程異常的變化不規(guī)則，用規(guī)則的曲面取代實際上不規(guī)則的似大地水準面，必然會導致大地水準面的不規(guī)則變化，導致擬合后高程異常的精度降低、誤差較大。多重分形可通過高程值計算出高程局部異常的奇異性指數，通過奇異性指數得到多重分形Kriging插值的結果。

局部奇異性分析方法實際上是將樣點高程值段的密度在分形空間中進行度量，以確定分形密度和分形維數，奇異性所度量的是場值隨量度范圍大小的變化規(guī)律。奇異性指數 (Z)表示基于對場值的某種滑動加權平均，計算公式如下:

其中:Ω(x0，ε)為圍繞中心點x0半徑ε的小滑動窗口;ω(‖x0-x‖)為對在Ω (x0，ε)中與中心點x0相隔‖x0-x‖距離的任意點x的加權函數，它往往與距離呈負相關。加權函數的選擇不僅與距離有關，與場的空間性自相關，還與處理目的有關。

成秋明［8］給出的多重分形方法將滑動平均關系表達為

其中，α(x0)為x0點處的局部奇異性指數?？梢钥闯觯?(2)中不僅包含了空間相關性的成分，還有度量奇異性指數。如果α(x0)=2時，通過該方法所計算的加權平均值與通常的加權平均無異;但當處于高程值差異較大地段而且局部地區(qū)具有奇異性，即α(x0)＜2時，通過該方法所得的結果將高于通常的加權平均結果;相反，當處于高程值差異較低地段，即α(x0)＞2時，通過該方法所得的結果將低于通常的加權平均結果。由此可見，該方法有利于加強峰谷值，對于具有奇異性的空間結構來說，傳統(tǒng)的插值方法不能很好地估計出峰谷值。指數α對應分形空間維數，用分形空間維數與正常的歐氏空間維數的差Δα=2-α即可表示分形密度與正常密度的空間維數的差異。

2 結果與分析

2.1 普通Kriging法的插值結果

利用探索性數據分析工具得知插值樣點集(圖1中黑色實心點)偏態(tài)值為-0．157，可認為樣點呈近正態(tài)分布。再對數據進行試驗半方差函數的計算和球面模型擬合 (步長取3500 m，組數為12組)，得到研究區(qū)域的海底地形插值結果。從普通Kriging法的插值結果 (圖2)可以看出，海州灣水下地形起伏不大，呈明顯條帶狀分布。提取柵格中檢查點位置的水深值，并與檢查點的初始值進行線性擬合 (圖3)。結果表明:水深為-15～0 m范圍內的點比較集中，說明插值結果與實際值比較接近;-25～-15 m范圍內的插值結果與實際值有一定偏離，但插值結果與實際值均分布在斜率k=1的直線附近，表明總體估計較為準確。

2.2 多重分形Kriging法的插值結果

圖1 樣本點分布示意圖Fig.1 The distribution of sampling points

圖2 普通Kriging法的插值結果Fig.2 Map of interpolation by Kriging method

圖3 基于普通Kriging法插值的交叉驗證圖Fig.3 Cross－examination by Kriging method

在普通Kriging插值方法的基礎上，為準確反映海州灣海底地形，度量水深分布的局部奇異性，引入多重分形理論。首先采用不同格網尺寸等級，滑動計算多重分形測度范圍內高程分布所包圍形體的體積，計算格網虛擬局部奇異性指數;然后回歸計算雙對數曲線的曲率，得到奇異性指數［即分形維數α(x0)］。具體計算過程在Matlab軟件中完成。當曲率大于指定的閾值 (在此設為0．9)時，雙對數關系成立;否則，取α(x0)=2。結合海州灣海岸線，得到了基于多重分形理論的Kriging插值結果。該插值結果 (圖4)與普通Kriging插值結果趨勢相同，都呈明顯的條帶狀分布。將檢查點位置的水深值與檢查點的初始值進行線性擬合得到交叉驗證圖 (圖5)。

圖4 多重分形Kriging法的插值結果Fig.4 Map of interpolation by multifractal Kriging method

圖5 多重分形Kriging法插值的交叉驗證圖Fig.5 Cross－examination of multifractal Kriging method

2.3 兩種插值方法的比較

將兩種插值方法的結果與驗證數據集進行對比。從表1可以看出:實際值中最大值為-0．10 m，最小值為-29．50 m;利用普通Kriging法插值后的最大值為-0．09 m，最小值為-31．16 m，可見，其對數據集的拉伸較為明顯;基于多重分形理論的Kriging法插值后的最大值為-0．10 m，最小值為 -3 0．0 2 m，該方法對數據集的拉伸較小。從標準差與均方根誤差上看，基于多重分形的Kriging方法均比普通Kriging方法小，說明其得到的插值結果更為準確。

表1 兩種插值方法的結果比較Tab.1 Comparison of two methods

另外，對比圖2與圖4的插值結果，發(fā)現圖2較圖4平滑，沒有保留異常點的水深信息。由于后者彌補了半變異函數局部平滑的缺陷，其能夠有效地勾勒局部細節(jié)，度量地表高程局部隆起和下陷的局部奇異特征，更好地反映了海州灣的局部地形。

3 海州灣水下地形的三維可視化與分析

為了更形象地刻畫海州灣海底的地形，利用Phong簡單光照模型，基于多重分形理論的Kriging方法的插值結果，通過調整光照 (設定光照方向為水平向135°，垂直向為45°)、顏色與選擇觀察方向，參考自然地物可視化的相關研究［18-22］，利用IDL語言開發(fā)繪制了海州灣水下三維地形圖(圖6)，從宏觀上較好地展示了海州灣的整體水下地形;除此之外，海灣里的幾個暗礁 (圖6中紫色矩形區(qū))也較好地呈現了出來。

由于古黃河的南徙，其攜帶的大量泥沙沉積不斷使陸地向前推進，最終造成原來的孤島云臺山(圖6中A點)與大陸相連［14，23］，從而形成了從A到B的“大喇叭口”［24］形開敞海灣。海岸地貌類型主要有海蝕地貌和海積地貌。如圖6所示，海域寬闊但水深較淺，低潮線以下多為水下暗坡，地勢向東北方向傾斜，平均比降在0．37‰。在灣口附近和口外淺海，水深為10～27 m，是一個起伏和緩的沖刷面，地勢向東偏北方向傾斜。利用2012年獲取的研究區(qū)域THEOS衛(wèi)星融合影像 (圖7)，對比圖6中暗礁的坐標，確定其對應著灣口東北部的基巖型島嶼:達山島、車牛山島和牛背島，周圍還存在著達東礁和牛嘴礁等。雖然采集的水深數據無法模擬出海平面以上的島嶼，但通過多重分形Kriging方法準確地擬合出了海底暗礁的位置。

4 結論

利用研究區(qū)域的水深數據點，分別使用普通Kriging插值法和基于多重分形的Kriging插值方法進行了海州灣的水下地形構建與三維可視化。研究表明:(1)普通Kriging插值法通過半變異函數從預測點周圍的觀測值中生成權系數進行預測，可以較好地反映整個區(qū)域內的數據空間特征，但插值結果較平滑，丟失了局部信息;(2)多重分形理論由于采用空間自相似性過濾，改進了普通Kriging方法中半變異函數局部平滑的缺點，考慮了空間數據分布的局部奇異性，可以獲得更好的空間插值效果。對于研究區(qū)域海州灣，采用多重分形Kriging插值方法在宏觀上較細膩地反映了平坦的海底地形，在微觀上也展示了達山島和車牛山島等周圍的暗礁。證明該方法可廣泛用于模擬復雜的幾何形體，尤其適用于復雜粗糙地形的重建和對數據異常的敏感度量場的重建。

圖6 研究區(qū)地形模擬結果Fig.6 Terrain simulation of the observed area

圖7 海州灣THEOS融合影像Fig.7 THEOS fusion image of Haizhou Gulf

［1］申靜，蘇天赟，王國宇，等．基于Kriging算法的海底地形插值設計與實現［J］．海洋科學，2012，36(5):24-28．

［2］包世泰，廖衍旋，胡月明，等．基于 Kriging的地形高程插值［J］．地理與地理信息科學，2007，23(2):28-32．

［3］張捷，包浩生．分形理論及其在地貌學中的應用——分形地貌學研究綜述及展望［J］．地理研究，1994，13(3):104-112．

［4］高義，蘇奮振，周成虎，等．基于分形的中國大陸海岸線尺度效應研究［J］．地理學報，2011，66(3):331-339．

［5］金德生，陳浩，郭慶伍．河道縱剖面分形——分線性形態(tài)特征［J］．地理學報，1997，52(2):154-162．

［6］熊波，王建力，張?zhí)煳模鍠|南巖溶山區(qū)耕地利用變化下土壤顆粒體積分形特征研究［J］．中國巖溶，2011，30(3):295-301．

［7］劉鵬舉，趙仁亮，朱金兆，等．保持地貌特征的數字高程模型生成方法研究［J］．中國礦業(yè)大學學報，2006，35(4):523-526．

［8］成秋明．多重分形與地質統(tǒng)計學方法用于勘查地球化學異?？臻g結構和奇異性分析［J］．地球科學，2001，26(2):161-166．

［9］Lee C K．Multifractal characteristics in air pollutant concentration time series［J］．Water，Air and Soil Pollution，2002，135(1/4):389-409．

［10］李錳，朱令人，龍海英．不同類型地貌的各向異性分形與多重分形特征研究［J］．地球學報，2003，24(3):237-242．

［11］曹漢強，朱光喜，李旭濤，等．多重分形及其在地形特征分析中的應用［J］．北京航空航天大學學報，2004，30(12):1182-1185．

［12］崔靈周，李占斌，郭彥彪，等．基于分形信息維數的流域地貌形態(tài)與侵蝕產沙關系［J］．土壤學報，2007，3(2):197-203．

［13］Wang D，Fu B J，Lu K S，et al．Multifractal analysis of landuse pattern in space and time:a case study in the Loess Plateau of China［J］．Ecological Complexity，2010(7):487-493．

［14］中國海灣志編纂委員會．中國海灣志．第四分冊［M］．北京:海洋出版社，1993:354-380．

［15］左書華，龐啟秀，楊華，等．海州灣海域懸沙分布特征及運動規(guī)律分析［J］．山東科技大學學報:自然科學版，2013，32(1):10-17．

［16］王寶燦，虞志英，劉蒼字，等．海州灣岸灘演變過程和泥沙流動向［J］．海洋學報，1980，2(1):79-96．

［17］張存勇，馮秀麗，陳斌林．海州灣南部近岸柱狀樣沉積物重金屬污染［J］．海洋地質與第四紀地質，2008，28(5):37-43．

［18］薛安，馬藹乃，李天宏．基于OpenGL實現真實感地形表現的研究［J］．中國圖象圖形學報，2001，6(2):800-805．

［19］孫博文．分形算法與程序設計——Delphi實現［M］．北京:科學出版社，2005:257-331．

［20］王夢，金文標．基于函數迭代系統(tǒng)的3-D分形插值算法［J］．計算機應用，2006，26(11):2701-2703．

［21］黃天云，張傳武．分形插值算法在分形自然景物模擬中的應用［J］．計算機工程與設計，2007，28(16):3394-3397．

［22］張濤，徐曉蘇，王其，等．基于分形插值的三維海底地圖生成算法［J］．中國慣性科學學報，2008，16(2):171-173．

［23］黃志強．全新世云臺山的海陸變遷［J］．徐州師范學院學報:自然科學版，1992，10(2):28-33．

［24］耿秀山．黃渤海地貌特征及形成因素探討［J］．地理學報，1981，36(4):423-434．