基于球坐標轉換的土壤機械組成空間插值研究

李　貞　張世文　曹　夢　黃亞捷　魏　星　黃元仿

(1.中國農業大學資源與環境學院， 北京 100193； 2.安徽理工大學地球與環境學院， 淮南 232001)

0　引言

土壤機械組成是土壤較為穩定的自然屬性，是研究有關土壤的一系列物理化學過程(如土壤改良、溶質運移、農業區劃等)的重要基礎，也是進行陸面過程及大氣過程模擬必需的重要參數。空間插值是獲取區域土壤屬性的一種重要方法，然而土壤機械組成作為一種成分數據，其空間插值不僅要滿足線性無偏、最優要求，而且要滿足定額與非負條件[1-2]，因此常常在插值前對土壤機械組成數據進行轉換。常用的轉換方法為對數比轉換方法[3-7]。雖然通過對數比轉換處理成分數據在各個學科中都有著廣泛的應用[8-9]，但是當成分數據中出現零值時，對數比轉換方法需要進行特殊處理，如李春軒等[10]在引入零值替換方法的基礎上比較了不同對數比轉換方法對成分數據空間插值效果的影響。有研究表明，球坐標轉換方法是一種新穎的成分數據處理方法[11]。趙江濤等[12]利用球坐標轉換和對數比轉換對北京市的就業情況進行預測建模，結果表明，球坐標轉換和對數比轉換均能取得很好的效果。WANG等[13]闡述了球坐標轉換方法并利用該方法對四維含有零值的成分數據進行了分析和預測，結果表明，球坐標轉換方法既可以降低數據的冗余程度，滿足成分數據的定額、非負條件，又無需考慮數據中的零值問題。土壤機械組成數據是一種三維數據，利用球坐標系表示更加直觀易懂，然而目前還沒有球坐標轉換方法對土壤機械組成數據空間插值效果的影響研究。因此，本文采用球坐標轉換方法對土壤機械組成數據進行轉換，利用轉換后的數據進行普通克里格插值，并與經對稱對數比和非對稱對數比轉換后的土壤機械組成數據的普通克里格插值結果進行比較，分析利用球坐標轉換方法進行土壤機械組成空間插值的適用性以及3種轉換方法對插值效果的影響。

1　研究區域概況與數據獲取

1.1　研究區域概況

選取寧夏銀北地區為研究區(105°51′1.83″～106°58′59.51″E，38°16′20.31″～39°23′26.99″N)。研究區地處寧夏河套灌區賀蘭山東麓、銀川平原北部和鄂爾多斯臺地西緣的高階地上，由西南向東北傾斜，主要地貌類型為山前洪積傾斜平原，黃河沖洪積平原和沖湖積平原，東西寬約51 km，南北長約130 km，面積約6 967.80 km2，行政區劃上主要包括銀川市和石嘴山市(圖1)。銀北地區屬于典型的溫帶大陸性氣候，日照充足，干旱少雨，蒸發強烈，風大多沙，多年平均降雨量為183～200 mm。土壤類型以灰鈣土和灌淤土為主，并包括少量的潮土、鹽土、風沙土、新積土、沼澤土、龜裂堿土和泥炭土，巖性構成從東到西依次為礫石、粗砂、中砂、細粉砂、砂黏土[14]。

1.2　數據獲取與處理

本研究采用網格均勻布點和分層抽樣相結合的方法在研究區內布點采樣。網格大小為2.5 km×2.5 km，在此基礎上根據土壤類型、土地利用類型、數字高程等進行分層抽樣。采樣時間為2017年4月，共取土樣184個，每個樣點用GPS(Global position system)記錄其位置，樣點分布見圖1。采用梅花取樣法，每個土樣都由取樣點附近直徑10 m范圍內5個土壤表層0～20 cm土樣混合而成。樣品經自然風干后過2 mm篩備用。土壤機械組成采用激光粒度分析儀測定，每個樣品測定3次，取其平均值。土壤顆粒分級標準采用國際制，粒徑在20～2 000 μm為砂粒，粒徑在2～20 μm為粉粒，粒徑0～2 μm為黏粒。

土壤機械組成數據的描述性統計、正態性檢驗采用SPSS軟件，土壤質地三角圖采用Origin完成，地統計學分析和空間分布圖制作以及交叉檢驗利用ArcGIS完成。

2　研究方法

2.1　土壤機械組成數據轉換方法

球坐標轉換(Spherical coordinate transform, SCT)是將成分數據由直角坐標系轉換到球坐標系。對于空間上第i個點上具有p種成分的成分ui=(ui,1,ui,2,…,ui,p),由于定和限制即ui,1+ui,2+…+ui,p=1,可對其各分量開根號,u′i,j=(ui,j)0.5(j=1,2,…,p)做簡單的非線性變換，此時有(u′i,1)2+(u′i,2)2+…+(u′i,p)2=1,則u′i=(u′i,1,u′i,2,…,u′i,p)∈Rp分布在一個半徑為1的p維超球面上，因此，可將ui=(ui,1,ui,2,…,ui,p)(i=1,2, …,n)從直角坐標系變換到球坐標系(r,θi,2, … ,θi,p)∈Θp，具體映射公式為

(1)

由球坐標系轉回直角坐標系為[13]

(2)

式中ui,j——第i個樣點上第j種成分的質量分數

θi,j——第i個樣點球坐標系上第j個分量，0<θi,j≤π/2，j=2,3, …,p

p——成分種類數，取3m——樣點數

2.1.2對數比轉換方法

AITCHISON[15-16]提出成分數據的對數比轉換方法，將成分數據變換成其主分的比值對數(稱對數比)，解決了成分數據統計分析中的閉合效應和非正態分布問題， PAWLOWSKY等[1]將對數比方法與地質統計學方法相結合，提出了成分數據的區域化統計方法，解決了成分數據插值的定額、非負、誤差最小和無偏估計的要求。常用的對數比轉換方法為非對稱對數比轉換(Asymmetry logratio transform, ALT)，有學者在此基礎上又進行了改進，提出了對稱對數比轉換(Symmetry logratio transform, SLT)[10,17]。非對稱對數比轉換及其轉回公式為

(3)

(4)

對稱對數比轉換及其轉回公式為

財務公司產業鏈金融服務對于企業集團以及產業鏈相關的企業都意義重大，能夠幫助各企業整體運營資源并優化發展計劃。財務公司需要明確在企業集團和產業鏈金融服務中的職能地位，憑借自身的優勢來不斷優化服務細節。在新時期財務公司的產業鏈金融服務應當更好地利用信息化手段并加強對潛在管理問題和風險的處理能力，保證相關金融服務的科學性和有效性，推動企業集團與其他企業更高水平地發展。

(5)

(6)

式中φi,j——ui,j經非對稱對數比轉換后的數值

βi,j——ui,j經對稱對數比轉換后的數值

ηj——常數，通常取所有樣點第j種成分除0以外最小質量分數的一半

2.2　土壤機械組成空間插值方法

地統計學利用采樣點的空間自相關性預測未知點的值，是區域制圖的一種重要方法。由于本文側重于不同轉換方法對土壤機械組成空間插值精度的影響，因此選擇最常用的普通克里格(Ordinary Kriging,OK)插值方法進行空間插值，該方法認為未知點的值可以通過已知點值的線性組合得到[18-19]。

2.3　預測精度檢驗

本文采用常用的交叉驗證方法進行插值精度檢驗。選取平均絕對誤差MAE、均方根誤差RMSE、一致性指標d、相關系數R評價不同轉換方法的插值精度，其中，平均絕對誤差和均方根誤差越小，插值精度越高，一致性指標和相關系數越大，插值精度越高[20-23]，計算方法為

(7)

(8)

(9)

其中

(10)

(11)

式中pi——第i個點的預測值

oi——第i個點的實測值

n——驗證點的個數，文中為采樣點數

3　結果與分析

3.1　土壤機械組成的描述性統計分析

由于銀北地區特殊的氣候、地質條件，使得該地區荒漠化和鹽漬化同時并存，因此本研究中部分樣點粉粒和黏粒的含量為零。由表1可知，整個研究區砂粒含量的分布區間為12.55%～100%，粉粒含量的分布區間為0～73.42%，黏粒含量的分布區間為0～22.32%。從土壤顆粒的變異系數可以看出，黏粒的變異性最大，達到75.84%，其次為粉粒53.42%，砂粒的變異系數最小，為28.38%。按照變異系數的劃分標準[24-25]，研究區3種土壤顆粒均處于中等強度變異水平。根據土壤質地三角圖(圖2)可以看出該地區主要為砂土及砂壤土，還有少量的粉砂粘壤，這與張秀珍等[26]的研究結果基本一致。

表1　研究區土壤機械組成統計特征Tab.1　Statistical characteristics of soil mechanical composition in study area

圖2　研究區土壤質地三角圖Fig.2　Ternary diagram for all samples

3.2　土壤機械組成空間插值

通過地統計學分析得到3種轉換方法轉換后數據半變異函數的擬合模型及相關參數，見表2。由表2可知，經對數比轉換后的砂粒、粉粒、黏粒分別符合球狀模型、高斯模型和指數模型。由于球坐標轉換方法的降維特性，因此經球坐標轉換后只有兩列數據參與分析和插值，其半變異函數分別符合球狀模型和高斯模型。由于非對稱對數比轉換無法處理原數據中的零值，因此僅有177個樣點參與地統計分析和插值。球坐標轉換中砂粒、粉粒、黏粒的排列方式影響轉換后數據分布，經比較，本文中以黏粒、砂粒和粉粒的排列方式為最優。

C0/(C0+C1)常用來表示變量空間自相關性強度，低于25%和高于75%表明變量分別具有較強和較弱的空間自依賴性[27]，由表2可知，除了經非對稱對數比轉換的黏粒和對稱對數比轉換的黏粒具有較強的空間自相關性以外，其他數據均具有中等程度空間自相關性，表明該地區土壤砂粒、粉粒的空間分布是由母質、氣候、地形等結構性因素以及耕地、施肥、排灌、管理措施等隨機因素共同作用的結果，而黏粒的分布受結構性因素影響更多一些。張世文等[4]對北京市土壤機械組成空間插值的研究結果也表明，對稱對數比轉換后的砂粒、粉粒呈中等空間自相關性，轉換后的黏粒具有較強的空間自相關性。

利用球坐標轉換、非對稱對數比轉換和對稱對數比轉換數據所獲得的土壤機械組成空間分布如圖3所示。3種轉換方法所獲得的砂粒、粉粒和黏粒在空間分布上顯示出相似的分布趨勢，如砂粒高值區域主要沿賀蘭山東麓和鄂爾多斯臺地由東南向西北方向呈帶狀分布；而粉粒的高值區域主要位于銀北地區中間地帶，由南向北沿黃河呈帶狀分布，這里地勢平坦，分布著大量肥沃的耕地；黏粒的高值區域則顯得較為零散，主要位于研究區東南角以及零星散布于研究區中間地帶。

表2　轉換后土壤機械組成數據半變異函數擬合模型及相關參數Tab.2　Semivariogram of transformed soil mechanical composition and parameters

圖3　基于不同轉換方法的土壤機械組成普通克里格插值圖Fig.3　OK prediction maps of soil mechanical composition for different transformation methods

3種轉換方法所得研究區土壤砂粒、粉粒、黏粒空間分布圖經柵格計算加和為1，表明球坐標轉換方法和對數比轉換方法一樣可以滿足土壤機械組成空間插值的定額、非負、線性無偏以及最優要求。

3.3　土壤機械組成空間插值精度檢驗

研究區樣點砂粒、粉粒與黏粒的實測值與3種轉換方法獲得的普通克里格預測值之間的平均絕對誤差MAE、均方根誤差RMSE、一致性指標d、相關系數R、實測值與預測值的散點圖見圖4。由圖4可知，對于砂粒和粉粒的預測，球坐標轉換方法的MAE、RMSE均低于對稱對數比轉換和非對稱對數比轉換方法，一致性指標d和相關系數R均高于其他兩種轉換方法，表明球坐標轉換方法對于砂粒和粉粒的空間預測效果相比對稱對數比和非對稱對數比轉換方法而言具有一定提高作用。而對于黏粒的空間預測，非對稱對數比轉換方法的MAE和RMSE均小于對稱對數比轉換和球坐標轉換，一致性指標d和相關系數R均高于對稱對數比轉換和球坐標轉換，表現出較好的預測效果。由實測值和預測值之間的散點圖也可以看出，球坐標轉換方法對砂粒和粉粒的空間預測效果略高于非對稱對數比轉換，而對于黏粒的空間預測，非對稱對數比轉換方法則表現出較好的預測效果。

圖4　不同轉換方法的交叉驗證分析Fig.4　Cross-validation analyses of different transformation methods

4　討論

研究了3種轉換方法對土壤機械組成普通克里格插值精度的影響，結果表明轉換方法的確對土壤機械組成的空間插值精度有一定的影響，這和前人研究結果一致[5-6,10]。研究結果顯示，整體而言，3種轉換方法對土壤機械組成空間插值的插值精度由大到小為黏粒、粉粒、砂粒，這和WANG等[28]基于對數比轉換的中國黑河流域土壤質地插值結果以及張世文等[5]基于對數比轉換的北京市土壤質地插值結果相同，均表現為黏粒的插值精度較高，這可能與黏粒含量較小有關；單就球坐標轉換方法而言，其對砂粒和粉粒的預測效果較好，而砂粒、粉粒、黏粒含量定和為1，砂粒、粉粒較接近真實值，從理論上而言黏粒的空間插值精度也應該較高，出現球坐標轉換對黏粒插值精度略低的原因可能與用于插值的經球坐標轉換后的數據有關。因此以后可針對不同轉換方法對土壤機械組成空間插值精度的影響機制做進一步研究。

空間結構影響插值精度，一般來說樣點間具有較強的空間自依賴性，其插值精度就比較高[29]。本研究中經非對稱對數比轉換的黏粒的C0/(C0+C1)小于對稱對數比轉換的黏粒，而且兩者均小于經球坐標轉換后的兩列數據的C0/(C0+C1)，表明經對數比轉換后的黏粒較經球坐標轉換后的兩列數據具有較強的空間自相關性，與此同時插值結果也表明經非對稱對數比轉換的黏粒的插值精度較高，預測值與實測值的相關系數達到0.981，一致性指標達到0.974，平均絕對誤差和均方根誤差分別為0.004和0.007，接近于0，對稱對數比轉換后的黏粒插值精度僅次于非對稱對數比轉換，而球坐標轉換后的黏粒預測效果最差。對于砂粒和粉粒而言，經球坐標轉換后數據的C0/(C0+C1)較對稱對數比轉換和非對稱對數比轉換后數據的C0/(C0+C1)均有一定程度的下降，而球坐標轉換對砂粒和粉粒的預測效果均有一定的提升，尤其相對于非對稱對數比轉換方法。因此，可以認為3種轉換方法對土壤機械組成普通克里格插值精度不同可能是由于轉換后用于插值的數據的C0/(C0+C1)不同所影響的。

插值方法是影響插值效果的重要因素[27,30]，有研究表明普通克里格插值方法計算復雜，需要的參數較多，因此其插值精度較低[30-31]，以后可以開展一些3種轉換方法對其他插值方法插值效果影響的研究。

5　結論

(1)以寧夏銀北地區土壤機械組成數據為例，驗證了球坐標轉換方法和對數比轉換方法同樣可以滿足成分數據空間插值的定額、非負、線性無偏以及最優要求，且球坐標轉換方法無需考慮零值，計算量少，是一種簡便且直觀的轉換方法。

(2)比較了球坐標轉換和對稱對數比轉換以及非對稱對數比轉換方法對土壤機械組成普通克里格插值效果的影響，結果表明，球坐標轉換方法對于砂粒和粉粒的預測精度最高，非對稱對數比轉換方法對于黏粒的空間預測效果最好。

1PAWLOWSKY V, OLEA R A, DAVIS J C. Estimation of regionalized compositions: a comparison of three methods[J].Mathematical Geology, 1995, 27(1):105-128.

2DE GRUIJTER J J, WALVOORT D J J, VAN GAMS P F M. Continuous soil maps—a fuzzy set approach to bridge the gap between aggregation levels of process and distribution models [J]. Geoderma, 1997, 77(2): 169-195.

3ZHANG Shiwen, KONG Weifang, HUANG Yuanfang, et al. Spatial prediction of topsoil texture in a mountain-plain transition zone using unvariate and multivariate methods based on symmetry logratio transformation[J].Intelligent Automation & Soft Computing, 2014, 20(1): 115-129.

4ZHANG Shiwen, SHEN Chongyang, CHEN Xiaoyang, et al. Spatial interpolation of soil texture using compositional Kriging and regression Kriging with consideration of the characteristics of compositional data and environment variables [J]. Journal of Integrative Agriculture, 2013, 12(9): 1673-1683.

5張世文,王勝濤,劉娜,等. 土壤質地空間預測方法比較[J]. 農業工程學報，2011,27(1): 332-339.

ZHANG Shiwen, WANG Shengtao, LIU Na, et al. Comparison of spatial prediction method for soil texture [J]. Transactions of the CSAE, 2011, 27(1):332-339. (in Chinese)

6檀滿枝,陳杰.土壤模糊隸屬度不同數據轉換方法及其對空間插值結果的影響[J]. 生態學報, 2009, 29(6)：3147-3153.

TAN Manzhi, CHEN Jie. Different transforms of fuzzy membership values of sampled soils and theirs influences on resulted interpolation prediction [J]. Acta Ecologica Sinica ,2009, 29(6): 3147-3153. (in Chinese)

7ODEH I O A, TODD A J, TRIANTAFILIS J. Spatial prediction of soil particle-size fraction as compositional data [J]. Soil Science, 2003, 168(7): 501-515.

8PAWLOWSKY-GLAHN V, EGOZCUE J J. Spatial analysis of compositional data: a historical review [J]. Journal of Geochemical Exploration, 2016, 164: 28-32.

9SMITH P F, RENNER R M, HASLETT S J. Compositional data in neuroscience: if you’ve got it, log it! [J]. Journal of Neuroscience Methods, 2016, 271: 154-159.

10李春軒,羅毅,包安明,等. 基于對數比轉換的成分數據空間插值研究[J].中國農業科學, 2012, 45(4): 648-655.

LI Chunxuan, LUO Yi, BAO Anming, et al. Study on spatial interpolation of compositional data based on log-ratio transformation [J]. Scientia Agricultura Sinica, 2012, 45(4): 648-655. (in Chinese)

11PAWLOWSKY-GLAHN V, EGOZCUE J J. Spatial analysis of compositional data: a historical review [J]. Journal of Geochemical Exploration, 2016, 164: 28-32.

12趙江濤,黃薇,王惠文. 兩種成分數據預測建模方法的比較研究[J].北京航空航天大學學報：社會科學版，2003, 16(2): 37-40.

ZHAO Jiangtao, HUANG Wei, WANG Huiwen. Contrastive study on two forecast modeling methods of compositional data [J].Journal of Beijing University of Aeronautics and Astronautics： Social Sciences Edition, 2003, 16(2): 37-40. (in Chinese)

13WANG Huiwen, LIU Qiang, HENRY M K, et al. A hyperspherical transformation forecasting model for compositional data [J]. European Journal of Operational Research, 2007, 179(2): 459-468.

14趙文娟. 寧夏銀北地區鹽漬土水鹽運移數值模擬研究[D].銀川:寧夏大學, 2014.

ZHAO Wenjuan. The study of numerical simulation application in soil water-salt flow in Yinbei region of Ningxia Province [D].Yinchuan: Ningxia University, 2014.(in Chinese)

15AITCHISON J. The statistical analysis of compositional data [J]. Journal of the Royal Statistical Society: Series B, 1982, 44(2): 139-177.

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

17檀滿枝,密術曉,李開麗,等. 不同插值方法對成分數據空間預測結果的影響[J]. 土壤, 2009, 41(6): 998-1003.

TAN Manzhi, MI Shuxiao, LI Kaili, et al. Influences of different interpolation methods on spatial prediction of compositional data [J]. Soils, 2009, 41(6):998-1003. (in Chinese)

18瞿明凱.幾種地統計學方法在縣域土壤空間信息處理上的應用與研究[D].武漢：華中農業大學, 2012.

QU Mingkai. Application and study of several geostatistical methods in soil spatial information processing at county scale [D]. Wuhan：Huazhong Agricultural University, 2012.(in Chinese)

19LIU Ruimin, CHEN Yaxin, SUN Chengchun, et al. Uncertainty analysis of total phosphorus spatial-temporal variations in the Yangtze River Estuary using different interpolation methods [J]. Marine Pollution Bulletin, 2014, 86(2): 68-75.

20GUMIERE S J, LAFOND J A, HALLEMA D W, et al. Mapping soil hydraulic conductivity and matric potential for water management of cranberry: characterization and spatial interpolation methods [J].Biosystems Engineering, 2014, 128: 29-40.

21XIE Yunfeng, CHEN Tongbin, LEI Mei, et al. Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: accuracy and uncertainty analysis [J].Chemosphere, 2011, 82(3): 468-476.

22葉回春,黃珊瑜,張世文,等. 土壤有機碳空間變異性對采樣密度的響應研究[J/OL]. 農業機械學報，2014，45(12): 215-222.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20141232&flag=1.DOI:10.6041/j.issn.1000-1298.2014.12.032.

YE Huichun, HUANG Shanyu, ZHANG Shiwen, et al. Spatial variability response of soil organic carbon to sampling density change [J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(12): 215-222. (in Chinese)

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

24張娜,張棟良,屈忠義,等. 內蒙古河套灌區區域土壤質地空間變異分析——以解放閘灌域為例[J]. 干旱區資源與環境, 2015, 29(12):155-163.

ZHANG Na, ZHANG Dongliang, QU Zhongyi, et al. The spatial variation of soil texture in Hetao Irrigation District in Inner Mongolia [J].Journal of Arid Land Resource and Environment, 2015, 29(12): 155-163. (in Chinese)

25王衛華,李建波,王鑠,等. 土壤熱特性參數空間變異性與擬合方法研究[J/OL]. 農業機械學報,2015,46(4):120-125. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20150418&flag=1.DOI:10.6041/j.issn.1000-1298.2015.04.018.

WANG Weihua, LI Jianbo, WANG Shuo, et al. Spatial variability of soil thermal parameters and its fitting method [J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2015,46(4):120-125.(in Chinese)

26張秀珍,劉秉儒,詹碩仁. 寧夏境內12種主要土壤類型分布區域與剖面特征[J].寧夏農林科技, 2011, 52(9):48-50,63.

27MIRZAEI R, SAKIZADEH M. Comparison of interpolation methods for the estimation of groundwater contamination in Andimeshk-Shush Plain, Southwest of Iran [J]. Environmental Science and Pollution Research, 2016, 23(3): 2758-2769.

28WANG Zong, SHI Wenjiao. Mapping soil particle-size fractions: a comparison of compositional Kriging and log-ratio Kriging [J]. Journal of Hydrology, 2017, 546: 526-541.

29KRAVCHENKO A N. Influence of spatial structure on accuracy of interpolation methods [J]. Soil Science Society of America Journal, 2003, 67(5): 1564-1571.

30LIU Ruimin, CHEN Yaxin, SUN Chengchun, et al. Uncertainty analysis of total phosphorus spatial-temporal variations in the Yangtze River Estuary using different interpolation methods [J]. Marine Pollution Bulletin, 2014, 86(2): 68-75.

31QU Linglu, XIAO Huayun, ZHENG Nengjian, et al. Comparison of four methods for spatial interpolation of estimated atmospheric nitrogen deposition in South China [J].Environmental Science & Pollution Research, 2017, 24(3): 2578-2588.