一種近海底質類型圖生成的非參數方法

劉付程，彭俊，張瑞，董春來，張存勇

（1. 淮海工學院 測繪工程學院，江蘇 連云港 222005；2. 華東師范大學 河口海岸學國家重點實驗室，上海 200062）

一種近海底質類型圖生成的非參數方法

劉付程1，彭俊2，張瑞1，董春來1，張存勇1

（1. 淮海工學院 測繪工程學院，江蘇 連云港 222005；2. 華東師范大學 河口海岸學國家重點實驗室，上海 200062）

近海底質類型圖在近海工程和經濟活動中有著廣泛的應用價值。針對傳統制圖方法中存在的問題，提出了一種基于非參數指示Kriging的底質類型圖生成方法。該方法能夠有效地規避制圖過程中的主觀性，且對取樣數據的平穩性和統計分布沒有特殊要求，并能對制圖結果的不確定性進行定量評價。該方法在連云港南部海域的應用實踐表明，在相同的條件下，該方法可獲得比傳統方法更為精確地的制圖結果，具有一定的實用價值。

沉積物類型；指示Kriging；制圖

1 引 言

海洋底質類型是指依據沉積物不同粒級組分的組成情況來對沉積物進行的一種類型劃分[1]。海洋底質類型圖是海洋沉積調查的重要成果圖件之一，它在現代海洋沉積環境分析、海洋工程建設以及浮標拋設與船舶錨泊地的選擇等方面都有著廣泛的實用價值[2,3]。

傳統的底質類型圖生成方法是根據取樣點的沉積物分類結果，結合水下地形信息和專家知識而繪制的。該方法很好地利用了已知信息和專家知識，但繪制過程中主觀性強，且工作效率低[2,4]。近些年來，隨著計算機技術和空間分析理論的發展，基于數據驅動的底質類型圖生成方法也得到了快速發展。張立華等[3]討論了運用 Voronoi圖生成技術來繪制底質類型圖的方法，楊康等[4]運用距離倒數權重（IDW）和克里格（Kriging）插值方法繪制沉積物不同粒級組分的空間分布圖，然后再采用柵格疊合技術制作底質類型圖。盡管這些方法在很大程度上規避了制圖過程的主觀性，但在實際應用過程中仍然還存在不同程度的局限性，如Voronoi圖法和 IDW 法對沉積物粒度空間分布的自相關性沒有給予足夠的重視，而 Kriging法則要求空間現象具備二階平穩或服從某種統計分布假設，這些都對制圖結果的可信度和方法的適用范圍產生不利影響。此外，目前大多數成圖方法都未能給出制圖結果不確定性的定量評價信息，這也在一定程度上制約了成果圖的有效應用。本文以連云港南部海域為例，嘗試運用非參數地質統計學中的指示克立格（Indicator Kriging）方法來對近海沉積物類型進行制圖，并對方法的適用性和制圖結果精確性作出評價。

2 非參數法制圖的基本原理

2.1 基本思想

對于研究海域的M個表層沉積物樣本，可依據一定的沉積物分類和命名方法（如Shepard分類法)將其劃分為N種類型。若這N種沉積物在該海域的某一個非取樣點處發生或出現的概率分別為P1、P2、…、PN，則存在如下的一個合理推斷：當P取最大值時，它所對應的沉積物類型即為該非取樣點處的沉積物類型。基于這一推斷，制圖的關鍵就轉變為如何估計不同類型沉積物在所有非取樣點處發生的概率，然后再通過比較概率值的大小來確定非取樣點處的沉積物類型，從而實現制圖。

若將沉積物取樣看成是對指定類型沉積物空間分布概率的一次實現，則相對于指定類型的沉積物而言，所有采樣點的概率值要么是 1（即該采樣點為指定類型沉積物），要么是0（即該采樣點為其它類型沉積物），據此可以采用非參數地質統計學中的指示 Kriging方法來估計指定類型沉積物在所有非取樣點處出現的概率。

2.2 指示克立格（Indicator Kriging）

指示克立格是條件普通Kriging的非參數形式，最早由Journel提出[5]。該方法可對非取樣點處的條件概率分布函數（Conditional Cumulative Distribu--tion Function, CCDF）進行估計，且該方法不依賴空間現象的平穩性，也不要求區域化變量服從某種分布假設，因此該方法在很多領域得到了廣泛地應用[6,7]。限于篇幅，本文只結合不同類型沉積物空間分布概率的估算來對其原理作一簡要的介紹，更詳細的內容，讀者可以參考文獻[5,8]。

假設在某海域進行沉積物取樣，經實驗室分析后獲得各樣品不同粒級組分的質量分數。若Z為基于指定分類方法（如Shepard分類法）的某一類型沉積物不同粒級組分的組合條件(即砂、粉砂和粘土的質量分數)，則在位置x處可定義如下的階梯函數(即指示函數)：

式中z(x)為位置x處的沉積物粒度組成情況。z(x)∈Z表示位置x處的粒度組成符合指定類型沉積物的粒度組成要求。在給定的Z條件下，z(x)∈Z時(即I(x,z)=1)的概率為

則指示函數I(x，Z)的期望值為

上式表明，位置x處的沉積物為指定類型沉積物的概率等于指示變量的平均值。當I(x+h, Z)和I(x, Z)為被矢量h（又稱步長）分割的兩個位置x+h和x處的指示函數值時，則指示變異函數γI(h, Z)可定義為

根據(4)式可求得指示變異函數相對于步長h的散點圖，該散點圖可用球狀、指數、高斯等理論模型來擬合，從而可以求得指示變異函數的擬合理論模型。據此可進一步采用普通Kriging方法來估計未知點 x處的 I(x,Z)值，也即指定類型沉積物在點 x處出現的概率，即

式中F*(Z)為滿足粒度組成條件Z的指定類型沉積物，在待估點x處的概率F的Kriging估計；n是待估點x鄰域中的樣點數；xk是鄰域中的第k個樣點的位置；λk是第k個樣點的權系數，它是在滿足最優無偏估計條件下，由指示變異函數來確定[5,8]。

3 試驗與分析

3.1 試驗數據

試驗區域為連云港港口以南、灌河口以北、低潮線和大約-12 m等深線之間的海域。2005年9月在該海域采集了106個表層沉積物樣品，樣點位置以GPS定位，采用垂直岸線的非網格化樣點布設方式，樣點分布見圖 1。沉積物樣品經實驗室分析獲得相應的粒度組成數據，再依據Shepard分類法對其進行類型劃分和命名[1]。為便于對制圖結果進行評價，所有樣本數據被隨機地分成2組，一組為插值數據集，用于制圖過程，共有 96個樣本數據；另一組為驗證數據集，用于檢驗制圖結果準確性，共有10個數據，其空間分布見圖1。

圖 1 采樣點分布圖及沉積物分類結果Fig. 1 Sampling sites and their sediment types based on Shepard textural classification system

3.2 結果與分析

3.2.1 不同類型沉積物指示變異函數的計算根據Shepard分類法，研究海域所有樣本數據可劃分為砂、粉砂質砂、砂-粉砂-粘土、粉砂、砂質粉砂和粘土質粉砂等6種類型沉積物（圖1）。據此可分別按照6種類型沉積物的粒度組成要求，對所有樣本進行指示變換，得到相應的指示變量值1和0，在此基礎上運用地統計學軟件 GS+7.0分別計算 6種類型沉積物的指示變異函數值，再經理論模型擬合后得到相應的理論模型模型參數（圖2）。

圖 2 不同類型沉積物的指示變異函數及其擬合模型和參數Fig. 2 Semi-variance of different types of sediment and their fitted models and parameters

從圖2可以看出，6種類型沉積物的指示變異函數均可以用球狀模型進行擬合，表明其空間變異均存在明顯的結構性特征。因此可以進一步利用Kriging方法對6種類型沉積物的指使變量（即出現概率）進行空間估計。

3.2.2 不同類型沉積物空間分布的條件概率估計圖3是基于圖2給出的不同類型沉積物指示變異函數的擬合理論模型，經GS+7.0中的Krigng插值方法得到的不同類型沉積物空間分布的條件概率圖。從圖3可以看出，不同類型沉積物空間分布的條件概率存在明顯的空間互補性，也即它們條件概率高值的分布在空間上是相互錯位的，這表明不同類型沉積物分別占住了研究海域的不同空間位置，這顯然對制圖是十分有利的。

圖3 不同類型沉積物空間分布的條件概率分布圖Fig. 3 Conditional probability distribution map of different types of sediments

3.2.3 底質類型圖的生成及其不確定性評價依據不同類型沉積物空間分布的條件概率圖，可采用硬化的方法獲得研究海域的底質類型圖。硬化是指將每一位置最大概率值所代表的沉積物類型作為該位置的沉積物類型。該過程可通過地理信息系統軟件ArcGIS9.2中的地圖代數運算功能來實現。圖4即硬化后的研究海域底質類型圖，該圖清楚地給出了研究海域不同類型沉積物的空間分布情況。

圖 4 連云港南部海域底質類型圖Fig. 4 Map for sediment types in Lianyungang southern sea area

由于硬化過程在接受一種類型沉積物的同時，也忽略了其它沉積物在該點發生的可能性，因此硬化過程產生一種忽略不確定性。忽略不確定性是與該位置沉積物類型的過渡性相關，當不同類型的沉積物在某點出現的概率越接近時，其忽略不確定性也就越大，制圖結果的可靠性也就越低，因此計算每點的忽略不確定性可生成相應的不確定性圖，該圖可以用于評價制圖結果的可靠性。忽略不確定性可以通過熵來表征[9]，熵表達了沉積物出現概率集中趨向某一特定類型沉積物的程度。熵的計算公式如下：

式中，Pikj是第k種沉積物在柵格單元(i, j)處出現的概率；l是研究海域已知的沉積物類型，本文中 l為6，Hij是點(i, j)處的熵，其取值范圍從0到1。當Hij為0時，表示點(i, j)處的沉積物完全屬于某個類型，在硬化過程中不會產生對其它類型的忽略。當Hij為1時，表示不同類型沉積物在該點處出現的概率相同，即沒有一種類型沉積物可以較好地代表該點的沉積物，因此該點劃分為任何一種類型的沉積物都會產生最大的忽略不確定性。根據式(6)，運用ArcGIS9.2的柵格運算功能可獲得研究海域的熵分布圖（圖5），也即忽略不確定性大小分布圖。

圖 5 底質類型圖的忽略不確定性Fig. 5 Ignored uncertainty of the sediment type map

從圖5可以看出，研究海域忽略不確定性在研究區的東南部和西北部相對較小，表明該海域制圖結果的可靠性高；而忽略不確定性在中部呈現斑塊狀高值分布，說明該斑塊區制圖結果的可靠性相對較低。結合圖1可以發現，忽略不確定性的高值分布區主要分布在具有多種類型沉積物的交界區，而低值區則主要分布在同一種類型沉積物的集中連片分布區域。這現象符合一般的推理結論，表明不確定性圖可以用于評價底質類型圖的可靠性。

3.2.4 制圖結果的精確性評價 為比較和評價制圖結果的精確性，表1給出了驗證數據集中所有樣點處沉積物類型的實測結果和制圖結果，表1同時還給出了按照文獻[3]和[4]中方法的制圖結果。

從表1可以看出，本文所提出的非參數方法只有5號樣點與實測結果不一致，其它9個樣點的制圖結果都與實測結果相同；而IDW和Kriging方法的制圖結果均有4個樣點（2、3、4、5號樣點）與實測結果不一致，表明非參數方法在制圖結果的精確性方面要好于IDW和Kriging方法。盡管Voronoi圖法的制圖結果也只有1號樣點與實測結果不一致（表1），但該方法對樣點的位置分布十分敏感，任何一個樣點位置的微小偏移都會造成制圖結果的顯著變化，因此該方法的穩健性不夠好；相反，非參數方法是基于指示變量空間變異的結構性特征來進行空間估計，它對位置敏感性相對較低，因此其適用性要較Voronoi圖法更為寬泛。

表 1 沉積物類型的制圖結果與真實結果的比較Tab. 1 Comparison between the mapping types and actual types of sediment

4 結 語

基于指示 Kriging的近海底質類型圖生成的非參數方法，能夠有效地規避制圖過程中的主觀性，且對沉積物取樣數據的平穩性和統計分布沒有特殊要求，并能對制圖結果的不確定性進行定量評價。該方法在連云港南部海域底質類型制圖的應用實踐表明，在相同的條件下，它可獲得比傳統方法更為精確地制圖結果，具有較強的適用性和應用價值。

需要指出的是，沉積物類型的空間分布受水下地形地貌的影響顯著，在本文的非參數成圖方法體系中，地形地貌因素未被考慮，因此該方法在水下地形變化較為復雜海域的應用效果還需要進一步檢驗。如何將沉積物類型的空間分布與水下地形地貌及水動力環境變化之間的關系納入到非參數成圖方法體系中，是值得進一步探討的問題。

[1] 國家海洋局908專項辦公室編. 海洋底質調查技術規范 [M]. 北京: 海洋出版社, 2006: 33.

[2] 劉錫清. 最新中國近海陸架底質類型圖 [J]. 海洋地質與第四紀地質. 1992, 12(4): 11-20.

[3] 張立華, 崔高嵩, 張建軍, 等. 一種海底底質與地形的信息疊置可視化方法及應用 [J]. 測繪科學, 2007, 32(4):111-112, 115.

[4] 楊康, 張永戰. 基于柵格疊合的沉積物底質圖生成方法 [J]. 第四紀研究, 2007, 27(5): 889-895.

[5] Journel A G. Non-parametric estimation of spatial distribution [J].Mathematical Geology, 1983, 15(3): 445-468.

[6] 劉瑞民, 王學軍, 張巍. 天津表土 PAHs區域環境風險評價研究[J]. 環境科學, 2008, 29(6): 1719-1723.

[7] 徐英, 陳亞新, 王俊生, 等. 農田土壤水分和鹽分空間分布的指示克立格分析評價 [J]. 水科學進展, 2006, 17(4): 477-482.

[8] 侯景儒. 指示克立格法的理論及方法 [J]. 地質與勘探, 1990, 26(3): 28-38.

[9] Zhu A X. Measuring uncertainty in class assignment for natural resource maps under a similarity model [J]. Photogrammetric Engineering & Remote Sensing, 1997,63:1195-1202.

A non-parametric method to generate type and distribution of coastal surface sediment map based on indicator kriging

LIU Fu-cheng1, PENG Jun2, ZHANG Rui1, DONG Chun-lai1, ZHANG Cun-yong1

(1. School of Geodesy and Geomatics Engineering, Huaihai Institute of Technology, Lianyungang 222005, China;
2. State Key Lab of Estuarine and Coastal Research, East China Normal University, Shanghai 200062, China)

Coastal surface sediment type map has been widely used in marine economic and engineering activities, but the traditional mapping methods had some shortcomings due to their intrinsic assumption. In this paper, a non-parametric method of indicator kriging has been proposed for generating types and distribution of coastal surface sediment. The method can effectively avoid mapping subjectivity, and has no special requirements for the sample data to meet second-order stationary or normal distribution, and can also provide useful information on the quantitative evaluation of mapping uncertainty. The application of the method in the southern sea area of Lianyungang showed that much more mapping precision could be obtained compared with the traditional methods such as IDW, kriging and Voronoi under the same condition, so the method had some promotion and practical values.

type of sediment; indicator kriging; mapping

P736.21

A

1001-6932(2011)05-0551-06

2010-12-03；

2011-04-26

江蘇省高校自然科學基金項目（07KJD170012）；江蘇省測繪科研項目(JSCHKY201005)；淮海工學院自然科學基金項目(Z2008009)。

劉付程（1971—），男，安徽安慶人，副教授，博士，主要從事海岸帶環境演變和GIS應用研究。電子郵箱：iliufucheng@126.com。