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

（1. 淮海工學(xué)院 測(cè)繪工程學(xué)院，江蘇 連云港 222005；2. 華東師范大學(xué) 河口海岸學(xué)國家重點(diǎn)實(shí)驗(yàn)室，上海 200062）

一種近海底質(zhì)類型圖生成的非參數(shù)方法

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

（1. 淮海工學(xué)院 測(cè)繪工程學(xué)院，江蘇 連云港 222005；2. 華東師范大學(xué) 河口海岸學(xué)國家重點(diǎn)實(shí)驗(yàn)室，上海 200062）

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

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

1 引 言

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

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

2 非參數(shù)法制圖的基本原理

2.1 基本思想

對(duì)于研究海域的M個(gè)表層沉積物樣本，可依據(jù)一定的沉積物分類和命名方法（如Shepard分類法)將其劃分為N種類型。若這N種沉積物在該海域的某一個(gè)非取樣點(diǎn)處發(fā)生或出現(xiàn)的概率分別為P1、P2、…、PN，則存在如下的一個(gè)合理推斷：當(dāng)P取最大值時(shí)，它所對(duì)應(yīng)的沉積物類型即為該非取樣點(diǎn)處的沉積物類型?；谶@一推斷，制圖的關(guān)鍵就轉(zhuǎn)變?yōu)槿绾喂烙?jì)不同類型沉積物在所有非取樣點(diǎn)處發(fā)生的概率，然后再通過比較概率值的大小來確定非取樣點(diǎn)處的沉積物類型，從而實(shí)現(xiàn)制圖。

若將沉積物取樣看成是對(duì)指定類型沉積物空間分布概率的一次實(shí)現(xiàn)，則相對(duì)于指定類型的沉積物而言，所有采樣點(diǎn)的概率值要么是 1（即該采樣點(diǎn)為指定類型沉積物），要么是0（即該采樣點(diǎn)為其它類型沉積物），據(jù)此可以采用非參數(shù)地質(zhì)統(tǒng)計(jì)學(xué)中的指示 Kriging方法來估計(jì)指定類型沉積物在所有非取樣點(diǎn)處出現(xiàn)的概率。

2.2 指示克立格（Indicator Kriging）

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

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

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

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

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

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

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

3 試驗(yàn)與分析

3.1 試驗(yàn)數(shù)據(jù)

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

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

3.2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

3.2.4 制圖結(jié)果的精確性評(píng)價(jià) 為比較和評(píng)價(jià)制圖結(jié)果的精確性，表1給出了驗(yàn)證數(shù)據(jù)集中所有樣點(diǎn)處沉積物類型的實(shí)測(cè)結(jié)果和制圖結(jié)果，表1同時(shí)還給出了按照文獻(xiàn)[3]和[4]中方法的制圖結(jié)果。

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

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

4 結(jié) 語

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

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

[1] 國家海洋局908專項(xiàng)辦公室編. 海洋底質(zhì)調(diào)查技術(shù)規(guī)范 [M]. 北京: 海洋出版社, 2006: 33.

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

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

[4] 楊康, 張永戰(zhàn). 基于柵格疊合的沉積物底質(zhì)圖生成方法 [J]. 第四紀(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] 劉瑞民, 王學(xué)軍, 張巍. 天津表土 PAHs區(qū)域環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)研究[J]. 環(huán)境科學(xué), 2008, 29(6): 1719-1723.

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

[8] 侯景儒. 指示克立格法的理論及方法 [J]. 地質(zhì)與勘探, 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

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

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