基于對數(shù)比轉(zhuǎn)換方法的沉積物粒級組分空間預(yù)測與底質(zhì)類型制圖

劉付程, 彭 俊

(1.淮海工學(xué)院 測繪工程學(xué)院, 江蘇 連云港 222005; 2.鹽城師范學(xué)院 城市與資源環(huán)境學(xué)院, 江蘇 鹽城224002)

底質(zhì)類型圖是海洋綜合調(diào)查成果的基礎(chǔ)圖件之一, 在海洋工程、漁業(yè)、養(yǎng)殖、航海、國防以及海洋沉積環(huán)境研究等方面都有著廣泛的應(yīng)用價值[1]。利用沉積物粒度分析數(shù)據(jù)開展底質(zhì)類型空間預(yù)測制圖是一種常見的制圖方式, 其技術(shù)途徑可分為兩類, 一是通過定量估計(jì)不同粒級組分(如砂、粉砂、黏土)的空間分布來進(jìn)行底質(zhì)類型識別和制圖表達(dá)[2], 二是通過采樣站位的底質(zhì)類型識別結(jié)果, 運(yùn)用Voronoi圖、指示 kriging等方法來進(jìn)行底質(zhì)類型邊界劃分[3-4]。由于前者更符合底質(zhì)制圖的一般邏輯推理過程, 因而廣泛被人們理解和接受, 但其關(guān)鍵是沉積物粒級組分的空間預(yù)測結(jié)果要可靠, 否則底質(zhì)類型識別的準(zhǔn)確性和制圖成果的可信度將難以保證。然而在沉積物粒級組分的空間預(yù)測過程中, 人們往往忽視其作為成分?jǐn)?shù)據(jù)所具有的特殊性, 這有可能導(dǎo)致與實(shí)際情況不相符的預(yù)測結(jié)果出現(xiàn), 如同一位置處砂、粉砂和黏土組分的預(yù)測結(jié)果之和不滿足“定和”條件[5]。

成分?jǐn)?shù)據(jù)是指具有“非負(fù)”和“定和”特性的一組數(shù)據(jù), 也即各組分取值≥0且其加和結(jié)果為一常數(shù)[6-7]。成分?jǐn)?shù)據(jù)的“定和”特性使得其各組分之間存在偽相關(guān), 并由此產(chǎn)生“閉合效應(yīng)”且不滿足經(jīng)典統(tǒng)計(jì)學(xué)方法的基本假設(shè)[6-9], 因此盲目采用經(jīng)典統(tǒng)計(jì)方法來處理成分?jǐn)?shù)據(jù)有可能得出錯誤的結(jié)論[7]。沉積物粒級組分?jǐn)?shù)據(jù)屬于典型的成分?jǐn)?shù)據(jù), 在其空間預(yù)測過程中,“非負(fù)”和“定和”要求是評判預(yù)測效果的重要依據(jù)之一, 也是后續(xù)開展底質(zhì)類型識別和制圖表達(dá)的前提條件, 然而在現(xiàn)有的多數(shù)文獻(xiàn)報(bào)道中, 這一點(diǎn)往往被有意或無意地忽略了。

成分?jǐn)?shù)據(jù)的“閉合效應(yīng)”使得直接運(yùn)用成分?jǐn)?shù)據(jù)的原始值來開展組分的空間預(yù)測受到質(zhì)疑[10]。著名的統(tǒng)計(jì)學(xué)家Aitchison教授提出了運(yùn)用對數(shù)比轉(zhuǎn)換方法來消除成分?jǐn)?shù)據(jù)的“閉合效應(yīng)”[8]。該方法將成分?jǐn)?shù)據(jù)變換成其組分比值的對數(shù), 從而實(shí)現(xiàn)了成分?jǐn)?shù)據(jù)從單形空間向?qū)崝?shù)空間的映射, 并使得轉(zhuǎn)換結(jié)果近似地服從正態(tài)分布[8-9], 這為運(yùn)用經(jīng)典統(tǒng)計(jì)學(xué)方法來處理和分析成分?jǐn)?shù)據(jù)創(chuàng)造了有利條件。本文以廢黃河三角洲表層沉積物粒度分析數(shù)據(jù)為基礎(chǔ), 嘗試運(yùn)用對數(shù)比轉(zhuǎn)換和kriging插值方法來對沉積物不同粒度組分的空間分布進(jìn)行預(yù)測, 并評價預(yù)測結(jié)果的準(zhǔn)確性, 在此基礎(chǔ)上開展底質(zhì)類型識別及制圖表達(dá)。

1 材料與方法

1.1 區(qū)域概況與數(shù)據(jù)

廢黃河三角洲由 1128~1855年間黃河南泛奪淮入海所帶來的泥沙淤積而成, 屬于典型的淤泥質(zhì)海岸。1855年黃河尾閭北歸后, 因喪失了主要泥沙來源, 三角洲被廢棄轉(zhuǎn)而進(jìn)入了侵蝕調(diào)整階段。廢黃河三角洲近岸海域平均流速為 0.25 m/s, 主流向?yàn)镹NW-SSE。據(jù)1995~2005年濱海海洋站波浪觀測資料統(tǒng)計(jì), 研究區(qū)常浪向?yàn)镋NE, 強(qiáng)浪向?yàn)镹E[11]。

2008年7~8 月, 在廢黃河口兩側(cè)布設(shè)了13條沉積物采樣斷面, 斷面從潮間帶開始, 向外海延伸5~6 km。在13條斷面上共采集表層沉積物樣品188個(圖1)。沉積物樣品經(jīng)處理后采用CoulterLS-100Q型激光粒度分析儀進(jìn)行粒度分析, 獲得全部樣品的砂(Φ＝–1~4)、粉砂(Φ＝4~8)和黏土(Φ>8)組分的含量。

在全部188個樣品中, 隨機(jī)抽取50個作為驗(yàn)證數(shù)據(jù)集, 以進(jìn)行沉積物粒度組分空間預(yù)測的準(zhǔn)確性評價; 其余的 133個樣品作為插值數(shù)據(jù)集用于沉積物粒級組分的空間預(yù)測(插值)。插值數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集的空間分布見圖1。

圖1 沉積物采樣站位分布圖Fig.1 The spatial distribution of sampling sites

1.2 對數(shù)比轉(zhuǎn)換方法

成分?jǐn)?shù)據(jù)對數(shù)比轉(zhuǎn)換方法主要有加和對數(shù)比轉(zhuǎn)換(additive log-ratio)、中心化對數(shù)比轉(zhuǎn)換(centered log-ratio)和等角對數(shù)比轉(zhuǎn)換(isometric log-ratio)等多種方法[8,12]。近些年來, 有學(xué)者為解決成分?jǐn)?shù)據(jù)中的0值問題(因0值無法取對數(shù))而對加和對數(shù)比轉(zhuǎn)換方法進(jìn)行了改進(jìn), 提出了改進(jìn)的加和對數(shù)比轉(zhuǎn)換方法(modified additive log-ratio)[10]。考慮到本次沉積物粒度分析數(shù)據(jù)中, 部分樣品的砂組分出現(xiàn)0值情況, 因此本文采用改進(jìn)的加和對數(shù)比轉(zhuǎn)換方法對沉積物的砂、粉砂和黏土組分進(jìn)行轉(zhuǎn)換, 其轉(zhuǎn)換公式為:

回轉(zhuǎn)公式為:

式中,xij為第i個沉積物樣品的第j種組分的含量;yij為第i個樣品的第j種組分含量的轉(zhuǎn)換值;ηj為常數(shù),取樣本中各組分?jǐn)?shù)據(jù)除0以外最小值的一半;c為沉積物組分的類別數(shù)。

1.3 空間預(yù)測效果評價

由于沉積物的粒級組分?jǐn)?shù)據(jù)為成分?jǐn)?shù)據(jù), 因此對其空間預(yù)測效果的評價既要考慮各粒級組分空間預(yù)測結(jié)果的準(zhǔn)確性, 同時還要考慮各組分預(yù)測結(jié)果是否滿足“非負(fù)”和“定和”的要求。

本文對不同粒級組分空間預(yù)測結(jié)果準(zhǔn)確性的評價采用平均絕對誤差(mean absolute error, MAE)、均方根誤差(root mean square error, RMSE)和一致性指標(biāo)(index of agreement,d)來進(jìn)行。平均絕對誤差和均方根誤差越小、一致性指標(biāo)越大表示預(yù)測效果越好[13]。各指標(biāo)的計(jì)算公式如下:

式中,n為驗(yàn)證數(shù)據(jù)集中的樣本數(shù);Xi、Zi分別表示驗(yàn)證數(shù)據(jù)集中的第i個樣本中某粒級組分的實(shí)測值和預(yù)測值;表示相應(yīng)粒級組分實(shí)測值的平均值。

對于“非負(fù)”檢驗(yàn)可直接由各粒級組分的預(yù)測結(jié)果來判別, 而對“定和”的檢驗(yàn)則是通過疊加各組分預(yù)測結(jié)果的柵格圖, 再逐個柵格判斷其 3組分預(yù)測值的加和是否為100%。

2 結(jié)果與分析

2.1 沉積物不同粒級組分含量的描述性統(tǒng)計(jì)

從表1可以看出, 研究海域沉積物組分中粉砂的平均含量最高, 為 55.86%, 超出了黏土和砂的平均含量之和, 反映了廢黃河三角洲沉積物的粉砂質(zhì)特性。3種組分的變異系數(shù)差別較大, 表明其空間分布的均勻性存在著顯著的差異。砂的變異系數(shù)達(dá)到了 1.216, 說明其空間分布的異質(zhì)性較強(qiáng); 粉砂的變異系數(shù)最小, 反映其空間分布均勻性要比砂和黏土更為明顯。表1中插值數(shù)據(jù)集與全體數(shù)據(jù)集有較為相似的統(tǒng)計(jì)參數(shù), 可將其視為全體數(shù)據(jù)集的一個理想子集。由插值數(shù)據(jù)集開展沉積物粒級組分的空間分布預(yù)測結(jié)果能在一定程度上反映全體數(shù)據(jù)集所體現(xiàn)的底質(zhì)空間分布特征。

表1 沉積物粒級組分含量的統(tǒng)計(jì)特征Tab.1 Statistical characteristics of different grain size compositions of the sampling sediments

2.2 沉積物粒級組分含量的空間分布預(yù)測

沉積物粒級組分的空間分布預(yù)測通常是通過空間插值來實(shí)現(xiàn)的。由于沉積物不同粒級組分的空間分布與位置有關(guān), 并表現(xiàn)出一定的空間自相關(guān)性,因此地統(tǒng)計(jì)學(xué)中的kriging插值方法常被認(rèn)為是一種理想的空間預(yù)測方法[14]。Kriging插值方法本質(zhì)上是一種加權(quán)平均方法, 其權(quán)重是在滿足最優(yōu)無偏估計(jì)條件下通過半方差函數(shù)求得的[15]。本文采用地統(tǒng)計(jì)學(xué)軟件GS+7.0對經(jīng)加和對數(shù)比轉(zhuǎn)換后的砂、粉砂和黏土數(shù)據(jù)進(jìn)行半方差計(jì)算和理論模型擬合, 模型參數(shù)由交叉驗(yàn)證法來確定。圖2給出了相應(yīng)的擬合理論模型及參數(shù)值(圖2a, b, c)。

圖2 基于不同數(shù)據(jù)處理方法的沉積物粒級組分理論半方差模型Fig.2 The fitted models for smivariograms of sediment grain size compositions based on data processed with different methods

圖2同時也給出了基于砂、粉砂和黏土原始數(shù)據(jù)的半方差散點(diǎn)圖及其擬合理論模型參數(shù)(圖2d, e,f)?？紤]到砂和粉砂原始數(shù)據(jù)的統(tǒng)計(jì)分布呈偏態(tài)(表1), 其半方差函數(shù)可能不夠穩(wěn)健[15], 有可能會影響到砂和粉砂的最終插值結(jié)果, 因此, 對砂和粉砂原始數(shù)據(jù)進(jìn)行了去除趨勢化處理, 即將插值數(shù)據(jù)集中各樣點(diǎn)砂和粉砂的實(shí)測值減去其各自的二階趨勢面值(選擇二階趨勢面是因?yàn)榻?jīng)檢驗(yàn)發(fā)現(xiàn)砂和粉砂含量的空間分布存在二階趨勢), 從而獲得砂和粉砂的二階趨勢面殘差。由于砂和粉砂二階趨勢面殘差的統(tǒng)計(jì)分布均近似為正態(tài), 由此計(jì)算得到的半方差函數(shù)更為穩(wěn)健。圖2g, h為砂和粉砂二階趨勢面殘差的半方差函數(shù)圖。

從圖2可以看出, 經(jīng)不同數(shù)據(jù)處理方法得到的砂、粉砂和黏土的變量半方差均可用指數(shù)模型或球狀模型來擬合, 表明變量的空間變異具有明顯的結(jié)構(gòu)性特征[15], 可進(jìn)一步運(yùn)用 ArcGIS的普通 kriging方法來對其空間分布進(jìn)行插值, 插值柵格邊長為30 m×30 m。為便于對不同數(shù)據(jù)處理方法下的預(yù)測結(jié)果進(jìn)行比較, 本文運(yùn)用 ArcGIS的柵格運(yùn)算功能, 將所有插值結(jié)果回轉(zhuǎn)到砂、粉砂和黏土的原始數(shù)據(jù)尺度, 即根據(jù)公式(2)將基于對數(shù)比轉(zhuǎn)換數(shù)據(jù)的插值結(jié)果恢復(fù)到其相應(yīng)的體積百分比數(shù)據(jù), 得到基于對數(shù)比轉(zhuǎn)換kriging方法的砂、粉砂和黏土空間預(yù)測結(jié)果分布圖(圖3a, b, c); 將砂、粉砂二階趨勢面殘差的kriging插值結(jié)果與其對應(yīng)的二階趨勢面進(jìn)行加和,得到基于趨勢面殘差kriging方法的砂、粉砂空間預(yù)測結(jié)果分布圖(圖3g, h)。圖3d, e, f是基于砂、粉砂和黏土原始數(shù)據(jù)的kriging插值結(jié)果分布圖。

從圖3可以看出, 基于不同數(shù)據(jù)處理方法所獲得的相應(yīng)粒級組分的空間分布格局基本一致, 砂的高值區(qū)主要分布在廢黃河口南側(cè)近岸區(qū)域, 呈與潮流方向一致的條帶狀分布; 粉砂的高值區(qū)主要分布在研究區(qū)域的西北部及東部; 黏土的低值區(qū)與砂的高值區(qū)基本對應(yīng), 其高值區(qū)主要分布在研究區(qū)域的北部近岸和東北部?；趯?shù)比轉(zhuǎn)換kriging方法的預(yù)測結(jié)果對 3組分空間分布的刻畫更為精細(xì), 比較圖3a, b, c和表1可以發(fā)現(xiàn), 其預(yù)測結(jié)果的上、下限值與插值數(shù)據(jù)集基本接近。基于原始數(shù)據(jù)的 kriging預(yù)測結(jié)果表現(xiàn)出顯著的平滑效應(yīng), 其值域范圍較插值數(shù)據(jù)集顯著縮小(圖3d, e, f和表1)。趨勢面殘差kriging方法對砂組分的預(yù)測結(jié)果出現(xiàn)負(fù)值(主要出現(xiàn)在東北部采樣點(diǎn)外側(cè)區(qū)域), 而對粉砂的預(yù)測結(jié)果超出了實(shí)測數(shù)據(jù)值域范圍的上限值, 與實(shí)測數(shù)據(jù)不符(圖3g, h和表1)。

圖3 基于不同數(shù)據(jù)處理方法的沉積物粒級組分預(yù)測圖Fig.3 Prediction maps of sediment grain size compositions based on different data processing methods

2.3 沉積物粒級組分含量的預(yù)測結(jié)果評價

表2是根據(jù)驗(yàn)證數(shù)據(jù)集中各樣點(diǎn)的預(yù)測值與實(shí)測值, 運(yùn)用式(3)~式(5)計(jì)算得到的 MAE、RMSE、d等評價指標(biāo)的值。從表2可以看出, 砂、粉砂和黏土組分在不同預(yù)測方法下的MAE、RMSE、d值雖然差異不大, 但3組分的MAE、RMSE均表現(xiàn)為對數(shù)比轉(zhuǎn)換 kriging法<趨勢面殘差 kriging法<原始數(shù)據(jù)kriging法, 而一致性評價指標(biāo)d的排序正好相反, 這表明僅從各粒級組分空間預(yù)測的準(zhǔn)確性來說, 對數(shù)比轉(zhuǎn)換kriging法要優(yōu)于另兩種方法。

表2 不同數(shù)據(jù)處理方法下的預(yù)測結(jié)果準(zhǔn)確性評價Tab.2 Evaluation of prediction accuracy at validation sampling stations based on different data processing methods

將同一方法的砂、粉砂和黏土的預(yù)測結(jié)果柵格圖進(jìn)行疊加, 獲得3組分“加和”分布圖, 可用來判別不同空間位置處 3組分預(yù)測結(jié)果是否滿足“定和”要求。圖4給出了不同預(yù)測方法的3組分“加和”分布圖,其中基于趨勢面kriging方法的“加和”分布圖中的黏土組分采用的是基于原始數(shù)據(jù)kriging方法的空間插值結(jié)果。從圖4可以看出, 基于對數(shù)比轉(zhuǎn)換 kriging方法的預(yù)測結(jié)果在所有空間位置處均滿足“定和”為100%的要求; 而另兩種方法只有不到 10%的區(qū)域面積滿足這一條件, 其絕大部分位置處的 3組分預(yù)測結(jié)果之和均大于或小于100%, 表明其預(yù)測結(jié)果只具有空間格局表達(dá)的參考價值, 不能作為定量分析的依據(jù)。

圖4 沉積物粒度組分預(yù)測結(jié)果加和分布圖Fig.4 Spatial distribution of the sum of predicted results of sand, slit and clay

2.4 基于沉積物粒級組分預(yù)測結(jié)果的底質(zhì)類型制圖

考慮到只有對數(shù)比轉(zhuǎn)換kriging方法對砂、粉砂和黏土的預(yù)測結(jié)果滿足“非負(fù)”和“定和”要求, 因此本文根據(jù)其預(yù)測結(jié)果, 按照Shepard沉積物分類和命名方案, 運(yùn)用 ArcGIS的柵格分析功能, 逐個柵格判別其所屬沉積物類型并最終繪制出研究區(qū)域的底質(zhì)類型分布圖(圖5)。從圖中可以看出, 研究海域存在7種底質(zhì)類型, 與實(shí)測數(shù)據(jù)的分類結(jié)果一致, 其中黏土質(zhì)粉砂的分布范圍最廣, 占插值區(qū)域面積的75%左右, 主要分布在廢黃河口北側(cè)近岸和5 m等深線以深海域; 其次是砂質(zhì)粉砂, 占插值區(qū)域面積的11%,主要分布在廢黃河口門外的鄰近海域; 粉砂質(zhì)砂占總面積的9%左右, 主要分布在廢黃河口南側(cè)近岸區(qū)域; 其他類型沉積物呈斑塊狀分布且面積均不大。

圖5 廢黃河口海域底質(zhì)類型分布圖Fig.5 Sediment type map of the coast of the abandoned Yellow River Delta

為進(jìn)一步評價制圖的精度, 將制圖結(jié)果與相應(yīng)實(shí)測樣本的沉積物類型進(jìn)行比較, 發(fā)現(xiàn)在全部 188個采樣站位處, 共有 174個站位的底質(zhì)類型與制圖結(jié)果是一致的, 即底質(zhì)制圖的總體精度達(dá)到了 92.6%,其中插值數(shù)據(jù)集和驗(yàn)證數(shù)據(jù)集中的一致性樣本數(shù)分別為 135個和 39個, 占其各自樣本總數(shù)的97.8%和78.0%, 表明底質(zhì)類型預(yù)測制圖的總體效果較好。

3 結(jié)論

(1) 沉積物粒級組分?jǐn)?shù)據(jù)屬于成分?jǐn)?shù)據(jù), 具有非負(fù)和定和特性?；诔练e物粒級組分原始數(shù)據(jù)的kriging預(yù)測方法難以保證預(yù)測結(jié)果的非負(fù)和定和要求, 其預(yù)測結(jié)果只能在一定程度上反映組分的空間分布格局, 不能作為定量分析的依據(jù)。而基于對數(shù)比轉(zhuǎn)換 kriging方法, 不但能確保預(yù)測結(jié)果符合成分?jǐn)?shù)據(jù)的基本要求, 而且還有著較高的組分預(yù)測準(zhǔn)確度,其預(yù)測結(jié)果可進(jìn)一步用于底質(zhì)類型識別和制圖表達(dá)等定量分析過程。

(2) 基于對數(shù)比轉(zhuǎn)換 kriging方法的砂、粉砂和黏土預(yù)測結(jié)果, 開展了廢黃河三角洲海域的底質(zhì)類型識別和制圖表達(dá)。制圖結(jié)果表明, 研究海域存在7種底質(zhì)類型, 與實(shí)測數(shù)據(jù)的分類結(jié)果一致, 其中以黏土質(zhì)粉砂分布最廣, 占插值區(qū)域的 75%左右, 且主要分布在廢黃河口北側(cè)近岸和5 m等深線以深海域; 砂質(zhì)粉砂主要分布在廢黃河口門外的鄰近海域,粉砂質(zhì)砂主要分布在廢黃河口南側(cè)近岸區(qū)域, 其他類型底質(zhì)分布面積較小且呈斑塊狀零星分布。基于實(shí)測數(shù)據(jù)的評價結(jié)果表明, 底質(zhì)制圖的總體精度達(dá)到了92.6%, 符合制圖要求。

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

[2] 楊康, 張永戰(zhàn).基于柵格疊合的沉積物底質(zhì)圖生成方法[J].第四紀(jì)地質(zhì), 2007, 27(5): 889-895.

[3] 王濤, 陳惠榮, 王少帥, 等.基于 Voronoi的海洋底質(zhì)區(qū)域劃界方法研究[J].測繪與空間地理信息, 2011,34(2): 242-243.

[4] 劉付程, 彭俊, 張瑞, 等.一種近海底質(zhì)類型圖生成的非參數(shù)方法[J].海洋通報(bào), 2011, 30(5): 551-556.

[5] Jerosch K.Geostatistical mapping and spatial variability of surficial sediment types on the Beaufort Shelf based on grain size data[J].Journal of Marine System, 2012, doi: 10.1016/j.jmarsys.2012.02.013.

[6] Aitchison J.The statistical analysis of compositional data[J].J Royal Stat Soc B, 1982, 44(2): 139-177.

[7] 周蒂.對數(shù)比統(tǒng)計(jì)分析及粒度數(shù)據(jù)中沉積水動力環(huán)境信息的萃取[J].沉積學(xué)報(bào), 1996, 14(增刊):149-157.

[8] Aitchison J.The statistical analysis of compositional data[M].London: Chapman and Hall, 1986: 58-61.

[9] 周蒂.地質(zhì)成分?jǐn)?shù)據(jù)統(tǒng)計(jì)分析—困難與探索[J].地球科學(xué)—中國地質(zhì)大學(xué)學(xué)報(bào), 1998, 23(2): 147-152.

[10] Odeh I O A, Todd A J, Triantafilis J.Spatial prediction of soil particle-size fractions as compositional data[J].Soil Science, 2003, 168 (7): 501-514.

[11] 陸勤.廢黃河三角洲淤泥質(zhì)海岸穩(wěn)定性研究[D].上海: 華東師范大學(xué), 2011.

[12] Egozcue J J, Pawlowsky-Glahn V, Mateu-Figueras G, et al.Isometric logratio transformations for compositional data analysis[J].Mathmatical Geology, 2003, 35(3):279-300.

[13] Willomtt C J.Some comments on the evaluation of model performance[J].Bulletin American Meteorological Society, 1982, 63(11): 1309-1313.

[14] 劉付程, 張存勇, 彭俊.海州灣表層沉積物粒度的空間變異特征[J].海洋科學(xué), 2010, 34(7): 54-59.

[15] 王政權(quán).地統(tǒng)計(jì)學(xué)及其在生態(tài)學(xué)中的應(yīng)用[M].北京:科學(xué)出版社, 1999: 59-149.