基于對數比轉換方法的沉積物粒級組分空間預測與底質類型制圖

劉付程, 彭 俊

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

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

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

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

1 材料與方法

1.1 區域概況與數據

廢黃河三角洲由 1128~1855年間黃河南泛奪淮入海所帶來的泥沙淤積而成, 屬于典型的淤泥質海岸。1855年黃河尾閭北歸后, 因喪失了主要泥沙來源, 三角洲被廢棄轉而進入了侵蝕調整階段。廢黃河三角洲近岸海域平均流速為 0.25 m/s, 主流向為NNW-SSE。據1995~2005年濱海海洋站波浪觀測資料統計, 研究區常浪向為ENE, 強浪向為NE[11]。

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

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

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

1.2 對數比轉換方法

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

回轉公式為:

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

1.3 空間預測效果評價

由于沉積物的粒級組分數據為成分數據, 因此對其空間預測效果的評價既要考慮各粒級組分空間預測結果的準確性, 同時還要考慮各組分預測結果是否滿足“非負”和“定和”的要求。

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

式中,n為驗證數據集中的樣本數;Xi、Zi分別表示驗證數據集中的第i個樣本中某粒級組分的實測值和預測值;表示相應粒級組分實測值的平均值。

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

2 結果與分析

2.1 沉積物不同粒級組分含量的描述性統計

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

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

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

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

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

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

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

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

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

2.3 沉積物粒級組分含量的預測結果評價

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

表2 不同數據處理方法下的預測結果準確性評價Tab.2 Evaluation of prediction accuracy at validation sampling stations based on different data processing methods

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

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

2.4 基于沉積物粒級組分預測結果的底質類型制圖

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

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

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

3 結論

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

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

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

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

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

[4] 劉付程, 彭俊, 張瑞, 等.一種近海底質類型圖生成的非參數方法[J].海洋通報, 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] 周蒂.對數比統計分析及粒度數據中沉積水動力環境信息的萃取[J].沉積學報, 1996, 14(增刊):149-157.

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

[9] 周蒂.地質成分數據統計分析—困難與探索[J].地球科學—中國地質大學學報, 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] 陸勤.廢黃河三角洲淤泥質海岸穩定性研究[D].上海: 華東師范大學, 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].海洋科學, 2010, 34(7): 54-59.

[15] 王政權.地統計學及其在生態學中的應用[M].北京:科學出版社, 1999: 59-149.