申 靜, 蘇天赟, 王國宇, 劉海行, 李家鋼

(1. 中國海洋大學 信息科學與工程學院, 山東 青島 266100; 2. 國家海洋局 第一海洋研究所, 山東 青島266061; 3. 中海石油研究總院, 北京 100027)

基于Kriging算法的海底地形插值設計與實現(xiàn)

申 靜1, 蘇天赟2, 王國宇1, 劉海行2, 李家鋼3

(1. 中國海洋大學 信息科學與工程學院, 山東 青島 266100; 2. 國家海洋局 第一海洋研究所, 山東 青島266061; 3. 中海石油研究總院, 北京 100027)

針對目前海底地形構建方面存在的問題與不足, 基于普通Kriging的計算公式, 以渤海海域采集的離散高程點數據為例, 對構建海底地形高程模型的空間插值方法——Kriging方法進行了研究。針對海底地形建模的具體實現(xiàn), 重點對Kriging算法的數據分布檢驗、數據分組、球形模型擬合、網格化插值以及結果顯示等模塊進行研究。最后, 基于Visual C++ 6.0平臺對海底地形高程的插值過程進行了編程實現(xiàn)和可視化表達, 從而成功構建出渤海局部海域的海底地形特征, 為海洋科學研究和工程建設提供了參考和依據。

Kriging; 插值; 地形; Visual C++; 可視化

構建高精度的海底地形高程模型是海洋科學研究和工程開發(fā)的基礎。由于海洋特殊的地理環(huán)境, 海底地形數據的獲取需要耗費大量的人力和物力, 人們只能獲得有關海底地形的部分離散水深數據。因此, 需要借助合適的插值算法, 基于采集到的離散水深點數據, 構建出合理的海底地形高程模型, 為海洋科學研究和工程開發(fā)提供基礎數據。

Kriging插值法, 又稱空間局部估計或空間局部插值法, 是地統(tǒng)計學的主要內容之一。它是以南非礦業(yè)工程師 D.G.Krige(Kriging)名字命名的一項實用空間估計技術, 目前已被廣泛應用于礦產儲量計算、遙感數據處理、地質學、水文學、環(huán)境科學、農林科學和農田水利等諸多方面。其插值實質是利用區(qū)域化變量的原始數據和變異函數的結構特點, 對未采樣點的區(qū)域化變量的取值進行線性無偏、最優(yōu)估計。與傳統(tǒng)插值方法相比它具有兩個明顯優(yōu)勢[1]: (1)在數據網格化過程中考慮了被描述對象的空間相關性,使估計結果更科學、更接近實際情況; (2)給出了估計誤差, 使估值精度更明了。

本文主要嘗試將Kriging插值方法應用到海底地形高程模型的構建中, 實現(xiàn)海底地形高程的合理插值與可視化。

1 Kriging插值原理

1.1 插值條件和變異函數

在研究區(qū)域內, Kriging方法將所研究的對象稱為區(qū)域化變量(如孔隙度、滲透率等), 區(qū)域化變量具有兩個最顯著的特征, 即隨機性和結構性[2-3]。用Kriging方法進行插值時, 一般需要區(qū)域化變量 z(x)滿足二階平穩(wěn)假設或區(qū)域化變量 z(x)的增量[z(x) - z (x + h )]滿足本征假設[4], 即: (1)在整個研究區(qū)域內, 區(qū)域化變量 z(x)的增量[z(x) - z (x + h )]的數學期望為 0, 即 E (z(x)- z(x +h)]=0, ?x, ?h; (2)在整個研究區(qū)域內, 增量[z(x) - z (x + h )]的方差函數存在且平穩(wěn)(不依賴于 x), 即Var[z(x) - z (x + h )]=E[z(x) - z(x + h ) ]2=2γ(h),?x, ?h

式中γ(h)稱為變差函數或變異函數(variogram),因為它為方差的一半, 所以有時又稱半變差函數或半變異函數(semivariogram), 它是地統(tǒng)計分析所特有的基本工具。h為兩樣本點空間分隔距離,z(xi)和z(xi+ h )分別是區(qū)域化變量 z(x)在空間位置 xi和xi+ h 處的實測值, N(h)表示距離為 h時數據點對的數目。由變異函數公式可以進一步推出如下關系式:

1.2 Kriging方程組

設z(x)是點x的區(qū)域化變量, 且滿足二階平穩(wěn)(或本征假設)。假設在點x的臨域內共有n個實測點, 即x1,x2,… ,xn, 其樣本值為z(xi)[i= 1 ,2,… ,n]。那么,普通Kriging法的插值公式為:

z*(x)表示區(qū)域化變量z(x)的無偏最優(yōu)估計量,λi為權重系數, 表示各空間樣本點xi處的觀測值z(xi)對估計值z*(x)的貢獻程度。普通Kriging方程組為:

寫成矩陣的形式為:AB=λ

由 Kriging的求解過程可知, 協(xié)方差矩陣(變異函數) 對于Kriging插值的影響較大, 估計變異函數是一個較為復雜的過程[5], 也是 Kriging內插的第一步。由于區(qū)域化變量的變異函數不依賴于位置x, 所以為了預測, 在實際計算中常需要用理論變異函數模型來替代實驗變異函數[6]。 Kriging法常用的理論變異函數模型有球形模型、指數模型、高斯模型、冪函數模型等, 其中球形模型是地統(tǒng)計分析中應用最廣泛的理論模型[3], 公式如(5)。

2 海底地形插值與可視化設計

Kriging插值算法適用于存在空間相關性的區(qū)域化變量, 海底地形符合隨機性和結構性的特征[7], 因此本文采用普通Kriging法對海底地形進行插值, 計算過程中認為海底地形的變化為“各向同性”, 并且采用 Visual C++對海底地形插值算法進行編程實現(xiàn)和可視化表達。在VC編程過程中, 為了便于程序的編寫和維護, 需要對算法進行模塊化分解, 然后對各模塊進行單獨編程和測試。程序處理流程如圖 1所示。

圖1 編程流程Fig. 1 Program flow chart

2.1 文件數據的讀取模塊

開發(fā)讀取文件函數 InputReader::Read(string filename, int mode)完成數據的讀操作, mode控制讀入數據的維數。本文通過建立可變長的矢量三維數組vectorm_input來存放數據的X,Y坐標和屬性值Z。vector為C++標準模板類和函數庫, 它能夠存放任意類型的動態(tài)數組。

2.2 數據分布檢驗

Kriging插值算法(如普通Kriging、簡單Kriging和泛 Kriging等)應用的前提是假設數據服從正態(tài)分布。因此, 在進行分析前, 檢驗數據分布特征, 了解和認識數據具有非常重要的意義。若不符合正態(tài)分布的假設, 應對數據進行變換, 使其轉為符合正態(tài)分布的形式, 常用的變換有 Log變換和 Box-cox變換。數據的檢驗可以通過直方圖完成。直方圖指對采樣數據按一定的分級方案(如等間隔分級)進行分級, 統(tǒng)計采樣點落入各個級別中的個數或占總采樣數的百分比, 并通過條帶圖或柱狀圖表現(xiàn)出來。

2.3 數據分組

2.4 擬合球形模型模塊

開發(fā) GetSphericalModel(const ForwardIterator start,int Size)函數擬合球形模型。由公式(5)可知當。如果記x= h3, 則可以得到線性模型:2

根據參數 b0,b1, b2的定義, b0＞ 0 ,b1＞ 0,b2＜0?？墒窃趯δＰ?6)的擬合過程中, 可能會出現(xiàn)參數符號不符合要求的情況, 從而無法正確地求出球形模型參數c0，c和a。因此, 本文針對擬合結果中b0,b1,b2系數符號不符合要求的情況做了詳細的分析處理:

(2)b0＜ 0 時, 令 b0= 0 , 解 yi= b1x1i+ b2x2i,

若求出 b1＞0,b2＜0, 則可得到球形模型參數;

若求出 b1＞0,b2＞0, 則令b2=0, 重新解yi= b1x1i,若求出 b1＞ 0 則可完成擬合, 否則擬合不成功。

(3)當 b0≥ 0 ,b1＞0, 但 b2≥ 0 時, 令 b2= 0 , 解yi= b0+ b1x1i,

綜上所述，隨著教育改革的不斷發(fā)展，要清楚認識到學生才是課堂中的主體。在教學中要從以學生為主體入手，運用互動教學法來激發(fā)出學生的學習主動性，培養(yǎng)好學生的學習能力，在促進學生思維發(fā)展的基礎上來提高語文教學質量。

若求出 b0＞0,b1＞0, 則可得球形模型參數;

若b0＜0,b1＞0, 則令b0=0, 重新解yi= b1x1i。

若求出 b1＞ 0 , 同樣可以完成擬合, 否則擬合不成功。

程序針對擬合不成功的情況下, 統(tǒng)一采用yi= x1i的線性擬合。

2.5 Kriging插值模塊

2.5.1 網格化

由于讀入的數據為散亂點, 需要對其進行規(guī)則化處理, 根據讀入海底數據的橫縱坐標的最大值、最小值和用戶設置的橫縱坐標步長即網格化的橫縱坐標間隔將散亂點網格化, 并將得到的坐標值 X,Y存入三維矢量網格數組 grid中。為了統(tǒng)一識別無效網格點, 每一個網格點 grid的屬性值 Z的初始值設為–9999。

2.5.2 確定已知插值點

對于預測點, 每個輸入點都有著局部影響, 這種影響隨著距離的增加而減弱。所以在進行模型擬合時, 可以利用預測點臨近的已知點進行擬合, 根據需要, 設置搜索半徑和插值所需的最少已知點個數。如果搜索半徑范圍內已知點的個數少于待插值點所需的個數時, 可以擴大搜索半徑, 繼續(xù)搜索。

2.5.3 插值計算

由擬合球形模型模塊可知, 求變異函數時, 只要求出參與插值的兩已知點的距離和待插點與已知點間的距離就能求出矩陣 A和B。另外在求解權重系數公式(4)時, 本文用到了矩陣的LU分解法[8]來加快矩陣的求解速度。權重系數求解完畢后, 帶入普通Kriging的求解公式(3)中完成海底地形高程的插值,將結果存入grid數組中將grid默認的–9999更改為實際的插值結果。

2.5.4 交叉驗證

為了驗證插值結果的誤差, 可以利用交叉驗證法[9](對數據集中任何一個樣本i( i=1, 2, …, n; n為樣本總數), 由其余n－1個樣本的觀測值對該樣本作出的估計值來評價預測精度), 驗證參數的合理性, 并進行相應的調整。

2.6 插值結果顯示模塊

開發(fā)位圖 DIB類, 根據插值得到的高程值數組grid中的屬性值的最大值和最小值, 設置調色板的色階分類和對應的屬性值區(qū)間。對每一個網格的屬性值, 可以從調色板中獲得對應的RGB顏色索引值,進而對該網格進行渲染。遍歷每一個網格進行繪制,即可得到直觀表達海底地形特征的可視化圖件。

3 海底地形插值應用實例

本文以實測水深數據為例, 基于設計的 Kriging地形插值算法, 對渤海局部海域的海底地形進行插值計算和可視化。

在計算過程中, 本文首先對實測數據進行了平穩(wěn)假設驗證。通過對數據直方圖(圖2)進行分析, 發(fā)現(xiàn)統(tǒng)計量偏度接近于0, 峰度接近于3, 均值與中值相近, 直方圖中柱狀呈對稱分布, 類似于一個倒扣的鐘, 數據近似于正態(tài)分布, 可以進行 Kriging插值。

圖 3給出了數據分組后擬合球形模型的參數(a= 6 6395.380,c0= 6 .421,c= 1 4.683)從圖中可以看出, 隨著距離的不斷增大, 變異函數的不斷遞增, 當增加到變程a時, 變異函數趨于平穩(wěn)。

圖2 水深數據分布直方圖Fig. 2 Distribution of water depths

圖3 交叉驗證圖Fig. 3 Cross-examination

分析交叉驗證圖(圖 4), 發(fā)現(xiàn)預測誤差相對較大的數值主要分布在與周圍實測值屬性值變化較大的點上, 也就驗證了Kriging插值法的實質是一個局部估計的加權平均, 具有一定的平滑作用, 尤其可以平滑掉一些極端值[10]。從圖5插值結果可以看出, 本文采用的Kriging方法能夠比較客觀、真實地反映海底地形起伏特征。

圖4 擬合球形模型圖Fig. 4 Spherical modeling

圖5 渤海局部水深地形Fig. 5 Local water depth terrain of the Bohai Sea

4 結語

通過算法設計、開發(fā)和應用實踐, 本文設計的Kriging插值程序能夠真實地反映出空間對象的變化特征。但相對于其他插值方法(如距離反比加權法) ,Kriging法復雜得多,需要根據具體的插值對象設置變異函數的參數和距離h的分類, 計算步驟繁瑣, 影響了插值速度。文章進行 Kriging算法編程過程中,對算法進行了優(yōu)化, 計算效率得到了提高。但是程序對插值點數進行限制是以犧牲精度來提高效率, 并且當處理的數據量比較大時, 容易出現(xiàn)內存溢出和運算效率降低的現(xiàn)象, 今后需要重點研究和解決Kriging算法在這方面存在的問題。

[1] 包世泰, 廖衍旋, 胡月明, 等. 基于 Kriging的地形高程插值[J]. 地理與地理信息科學, 2007, 23(3):28-32.

[2] 趙俊蘭. Kriging法在GIS空間數據內插中的應用[J].有色金屬(礦山部分), 1998, 3: 35-38.

[3] 李曉軍, 王長虹, 朱合華. Kriging 插值方法在地層模型生成中的應用[J]. 巖土力學, 2009, 30(1):157-162.

[4] 王靖波, 潘懋, 張緒定. 基于 Kriging方法的空間散亂點插值[J]. 計算機輔助設計與圖形學學報, 1999,11(6): 525-529.

[5] 顏輝武, 祝國瑞, 徐智勇, 等. 基于 Kriging 水文地質層的三維建模與體視化[J]. 武漢大學學報(信息科學版), 2004, 29(7): 611-614.

[6] 鮮乾坤, 萬寧, 賀全兵, 等. 基于 Kriging插值的 3D地形可視化[J]. 四川理工學院學報(自然科學版),2009, 22(6): 47-53.

[7] 李國敏, 陳崇希. 克立格法在地下水數值模型中的應用[J]. 地質科技情報, 1995, 14(3): 95-99.

[8] 紀坤, 陳建平. LU分解分塊算法的研究與實現(xiàn)[J].南京師大學報(自然科學版), 2008, 31(計算機專輯):55-60.

[9] 吳學文, 晏路明. 普通 Kriging 法的參數設置及變異函數模型選擇方法——以福建省一月均溫空間內插為例[J]. 地球信息科學2007, 9(3): 104-108.

[10] 尉英華, 郭品文, 劉洪濱. 利用插值法建立歷史旱澇格點資料的可行性[J]. 氣象與減災研究, 2007, 30(3):1-6.

Submarine topography visualization based on Kriging algorithm

SHEN Jing1, SU Tian-yun2, WANG Guo-yu1, LIU Hai-xing2, LI Jia-gang3
(1. Ocean University of China, the College of Information Science and Engineering, Qingdao 266100, China;2. First Institute of Oceanography, State Oceanic Administration, Qingdao 266061, China; 3. Petroleum Research Institute, Beijing 100027, China)

Oct.,24,2011

Kriging; interpolation; topography; Visual C++; visualization

On the basis of Kriging interpolation algorithm, we constructed the sea floor elevation model using discrete bathymetry points. In the model constructing, key modules including normal distribution verification, data binning, spherical model fitting, grid interpolation and grid visualization were designed. The visualization program and key modules were developed on the Visual C++ 6.0 platform. Part of the Bohai Sea elevation model was constructed using discrete bathymetry survey points.

P208; TP391

A

1000-3096(2012)05-0024-05

2011-10-24;

2012-03-06

“十二五”國家科技重大專項子任務(2011ZX05056-001-01)

申靜(1987-), 女, 山東濰坊人, 中國海洋大學在讀研究生,E-mail: shenjingshenwei@163.com

(本文編輯：劉珊珊)