基于球坐标转换的土壤机械组成空间插值研究

李　贞　张世文　曹　梦　黄亚捷　魏　星　黄元仿

(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.