地下水化学成分分类是进行水文地球化学研究的基础工作，通过分析区域地下水化学类型差异，可了解地下水化学成分的基本构成、形成过程及其受环境要素影响的演化规律^[1-6]。在地下水化学类型特征分析中，舒卡列夫分类是一种广泛应用的分类方法^[6]。舒卡列夫分类法根据水中6种主要的阴阳离子(K⁺+Na⁺、Ca²⁺、Mg²⁺、HCO₃^-、SO₄^2-、Cl^-)进行划分，将大于25%毫克当量百分数的离子参加分类命名，阴离子在前，阳离子在后，中间用短横线相连来对地下水化学类型进行命名^[6]。

地下水化学类型分区图可以较好地反映地下水化学成分的空间分布规律和演化特征^[6]。目前，基于舒卡列夫分类的地下水化学类型空间分区方法主要为泰森多边形法和离子毫克当量百分数的插值分区法^[6]。泰森多边形又称Voronoi图，是一种常见的离散数据空间插值方法，对于离散分布的样本点，通过建立泰森多边形可以保证每个多边形内分布一个样本点，样本点的值即代表多边形内区域的值，在水文、气象等研究中都有其应用^[6-9]。利用泰森多边形进行水化学类型分布特征分析时，多边形内的水化学类型取样本点的值即可。该方法的优点是简单易操作，缺点是空间插值没有考虑研究对象的空间连续性，当样本点数量有限或分布不均匀时，插值结果与实际情况差异较大，相邻区域内常发生水化学类型突变。

基于GIS的离子毫克当量百分数插值分区法是近年来应用较为广泛的一种水化学类型分布特征分析法^[10-11]，过程如下：①计算样本点阴阳离子的毫克当量百分数；②对阴阳离子毫克当量百分数进行空间插值，获得各个离子毫克当量百分数的区域分布数据；③将毫克当量百分数大于等于25%的区域进行提取；④对提取后的区域进行叠加分析；⑤制图综合获得水化学类型分区图。该方法的优点是考虑到了区域内水化学成分分布的空间连续性和变化趋势，通过离散样本点数据插值获得连续分布，进而计算出水化学类型分布特征。其缺点是过程复杂、效率较低，不能反映出离子间的含量大小排序关系，表达的水化学类型有限。

为了简化地下水化学类型分区的流程，客观反映地下水化学类型的空间分布特征，本文分析了常用方法的特点并提出改进的分区思路，设计并建立了基于GIS技术的地下水化学类型分区空间分析模型，通过实际样本数据进行了模型验证，对不同方法的分析结果进行了对比分析。

1 模型原理与方法 1.1 建模思路

在GIS中，栅格数据结构以编码来区分地物属性，拥有相同属性编码的连续栅格单元构成区域，在科学计算、空间分区与可视化表达方面具有简单高效的特点。根据栅格数据特点，提出以栅格单元为基本单位对6种离子毫克当量百分数空间插值数据进行水化学类型分区的分析思路，并利用GIS空间分析建模技术将分析过程流程化处理，实现地下水化学类型的自动识别与空间分区，从而解决常用方法存在的问题。模型分析过程如下：

1) 对离散的阴阳离子毫克当量百分数数据进行空间插值，以获取其空间连续分布栅格数据；

2) 计算每一个栅格单元含量大于25%的阴离子(阳离子)组合编码；

3) 按照含量大小排序，计算每一个栅格单元的阴离子(阳离子)排序编码；

4) 根据离子组合编码和离子排序编码，计算栅格单元的水化学类型编码；

5) 对水化学类型编码栅格数据进行地图符号化，得到最终的地下水化学类型分区图。

步骤1)属于空间数据准备范畴，而步骤5)属于对模型分析结果的后处理过程，空间分析模型的构建重点围绕步骤2)~4)展开。

1.2 模型原理

在使用GIS分析水化学类型时涉及到自动分类与排序问题，该过程主要通过编码来实现。

1.2.1 组合编码

舒卡列夫分类法中的49种水化学类型，是源自阳离子7种基本组合和阴离子7种基本组合的交叉构成，为便于GIS分析，对参与分类命名的阴阳离子组合进行编码，如表 1所示。为了便于表达，K⁺+Na⁺、Ca²⁺、Mg²⁺、HCO₃^-、SO₄^2-、Cl^-分别用字母N、C、M、H、S、L代表。

模型分析过程中，对插值得到的6个毫克当量百分数栅格数据，以25%为界限，按照表 1对阳离子和阴离子进行组合编码，并输出相应的组合编码栅格数据。例如，某栅格单元内阴离子HCO₃^-、SO₄^2-、Cl^-的毫克当量百分数均高于25%，则阴离子组合编码为3。

1.2.2 排序编码

获取离子的组合编码后，需要按照毫克当量百分数含量进行排序，根据排列组合可知阴阳离子共有12种基本排序类型，分别用表 2中编码表示，表中离子排序仅表示毫克当量百分数含量的高低顺序，并不代表 3种离子都符合25%的条件。

模型分析过程中，对插值得到的6个毫克当量百分数栅格数据，以栅格单元为基本单位，按照毫克当量百分数数值大小，按照表 2对阳离子和阴离子进行排序编码，并输出相应的排序编码栅格数据。例如，某栅格单元内阴离子的毫克当量百分数含量为HCO₃^- > SO₄^2- > Cl^-，则阴离子排序编码为1。

1.2.3 水化学类型编码

在获得离子组合编码和离子排序编码后，需要进一步确定其水化学类型，此时需要分别确定阴离子和阳离子的水化学类型编码，该编码应既能确定哪些离子符合25%条件, 又能表明含量排序关系。对此，在阴阳离子14种组合(见表 1)的基础上，每一种参与分类命名的离子组合用一个3位数的整型数字来表达，每位数字代表一种离子，当某一离子含量低于25%则用0代替。6种离子编码如表 3所示。至此，考虑到含量排序，参与分类命名的阴阳离子组合分别扩充到15种类型，其水化学类型编码如表 4所示。

模型分析过程中，对组合编码栅格数据和排序编码栅格数据，按照表 4分别判断出参与分类命名的阴离子和阳离子的水化学类型，并输出相应的阴阳离子水化学类型编码栅格数据。例如，某栅格单元内阴离子组合编码为3，排序编码为1，则阴离子水化学类型编码为123。

通过使用以上编码方案，对于舒卡列夫分类中的49种水化学类型细分为225种类型，可以使用下式进行编码组合：

$ C_{\mathrm{SK}}=C_{\text { type_anion }} \times 10000+9000+C_{\text { type_cation }} $

(1)

式中，C_SK表示舒卡列夫分类编码；C_{type_cation}表示阳离子水化学类型编码；C_{type_anion}表示阴离子水化学类型编码。

例如，编码1239320，123在表 4中对应HSL，320在表 4中对应NM，数字9为分隔符，其对应水化学类型即为HCO₃·SO₄·Cl-Na·Mg。

1.3 模型构建

目前，空间分析建模主要依赖于基于工作流技术的图形建模方法^[12-18]。本研究选择ESRI ArcGIS平台中ModelBuilder空间分析建模工具完成具体建模工作。在模型中，分别定义不同的图形代表输入数据、输出数据和地理处理工具，以流程图的形式组合执行空间分析操作功能^[17]。

根据建模思路，选择Map Algebra、Math等地理处理工具，利用ModelBuilder建立水化学类型分区空间分析模型，如图 1所示。其中，阴阳离子的组合、排序和类型划分使用条件函数(condition，Con)构造判定树来计算编码，水化学类型划分使用Math工具按式(1)计算编码。Con函数为：

$ \begin{array}{l} {\mathop{\rm Con}\nolimits} ( < {\rm{ condition }} > , \\ < {\rm{ true\_expression }} > , \\ < {\rm{false\_expression}} > ) \end{array} $

(2)

式中，condition表示判别条件；true_expression表示判别条件为true的表达式，为具体数值或嵌套Con函数；false_expression表示判别条件为false的表达式，为具体数值或嵌套Con函数。函数的输入、输出均为栅格数据。

1.3.1 组合编码

以图 1模型中阴离子组合编码为例，按式(2)使用Con函数进行分析的地图代数见式(3)~(6)，公式中变量含义见图 1注释。

$ \begin{array}{l} {\rm{COMB}}\_{\rm{anion }} = {\rm{Con}}\left( {{\rm{pe}}{{\rm{r}}_ - }{\rm{HCO}}3 > = 25} \right., \\ {\mathop{\rm Con}\nolimits} \left( {{\mathop{\rm per}\nolimits} \_{\rm{S}}{{\rm{O}}_4} > = 25, {\mathop{\rm Con}\nolimits} \left( {{{{\mathop{\rm per}\nolimits} }_ - }{\rm{Cl}} > = 25, 3, 2} \right)} \right., \\ {\mathop{\rm Con}\nolimits} \left( {{\mathop{\rm per}\nolimits} \_{\rm{Cl}} > = 25, 4, 1} \right)), \\ {\mathop{\rm Con}\nolimits} \left( {{\mathop{\rm per}\nolimits} \_{\rm{S}}{{\rm{O}}_4} > = 25, {\mathop{\rm Con}\nolimits} \left( {{\rm{pe}}{{\rm{r}}_ - }{\rm{Cl}} > = 25, 6, 5} \right)} \right., \\ {\mathop{\rm Con}\nolimits} \left( {{\mathop{\rm per}\nolimits} \_{\rm{Cl}} > = 25, 7, 0} \right))) \end{array} $

(3)

1.3.2 排序编码

阴离子排序编码的地图代数公式为：

$ \begin{array}{l} {\rm{sort\_anion}} = {\mathop{\rm Con}\nolimits} \left( {{\rm{pe}}{{\rm{r}}_ - }{\rm{HC}}{{\rm{O}}_3} > = {{{\mathop{\rm per}\nolimits} }_ - }{\rm{S}}{{\rm{O}}_4}} \right., \\ {\mathop{\rm Con}\nolimits} \left( {{{{\mathop{\rm per}\nolimits} }_ - }{\rm{S}}{{\rm{O}}_4} > = {{{\mathop{\rm per}\nolimits} }_ - }{\rm{Cl}}, 1, {\mathop{\rm Con}\nolimits} \left( {{\rm{pe}}{{\rm{r}}_ - }{\rm{HC}}{{\rm{O}}_3} > } \right.} \right.\\ = {{\mathop{\rm per}\nolimits} _ - }{\rm{Cl}}, 2, 5)), \\ {\mathop{\rm Con}\nolimits} \left( {{{{\mathop{\rm per}\nolimits} }_ - }{\rm{HC}}{{\rm{O}}_3} > = {{{\mathop{\rm per}\nolimits} }_ - }{\rm{Cl}}, 3, {\mathop{\rm Con}\nolimits} \left( {{\rm{pe}}{{\rm{r}}_ - }{\rm{Cl}} > } \right.} \right.\\ = {{\mathop{\rm per}\nolimits} _ - }{\rm{S}}{{\rm{O}}_4}, 6, 4))) \end{array} $

(4)

1.3.3 离子水化学类型编码

在获得COMB_anion和sort_anion基础上，进行阴离子类型编码，地图代数公式为：

$ \begin{array}{l} {\rm{type\_anion = Con(COMB\_anion = = 1, }}\\ {\rm{100, }}\\ {\rm{Con(COMB\_anion = = 2, Con}}\left( {{\rm{sort\_anion = }}} \right.\\ {\rm{ = 1, 120, Con}}\left( {{\rm{sort\_anion = = 3, 210, 0}}} \right){\rm{), }}\\ {\rm{Con(COMB\_anion = = 3, Con}}\left( {{\rm{sort\_anion = }}} \right.\\ {\rm{ = 1, 123, Con(sort\_anion = = 2, 132, }}\\ {\rm{Con}}\left( {{\rm{sort\_anion = }}} \right.\left( {{\rm{sort\_anion = }}} \right.{\rm{ = 3, 213, Con}}\left( {{\rm{sort\_anion}}} \right.\\ {\rm{ == 4, 231, Con(sort\_anion = = 5, 312, Con(sort}}\\ {\rm{\_anion = = 6, 321, 0)))))), }}\\ {\rm{Con(COMB\_anion = = 4, Con}}\left( {{\rm{sort\_anion}}} \right. = \\ {\rm{ = 2, 130, Con}}\left( {{\rm{sort\_anion = = 5, 310, 0}}} \right){\rm{), Con}}\\ {\rm{(COMB\_anion = = 5, 200, }}\\ {\rm{Con(COMB\_anion = = 6, Con(sort\_anion = }}\\ {\rm{ = 4, 230, }}\\ {\rm{Con}}\left( {{\rm{sort\_anion = = 6, 320, 0}}} \right){\rm{), }}\\ {\rm{Con}}\left( {{\rm{COMB\_anion = = 7, 300, 0}}} \right){\rm{))))))}} \end{array} $

(5)

1.3.4 水化学类型编码

阳离子COMB_cation、sort_cation、type_cation的地图代数公式与阴离子类似，在完成相应的编码计算后，按式(1)进行水化学类型划分，地图代数公式如下：

$ \begin{array}{c}{\text { hydrochemtype }=\text { type _} \text { anion } \times 10000+} \\ {9000+\text { type_cation }}\end{array} $

(6)

图 2所示为ModelBuilder在建立的水化学类型分区空间分析模型基础上自动生成的可执行对话框，在指定各模型参数对应的数据路径后，即可按照模型逻辑完成相应的水化学类型分析。

2 模型实例验证 2.1 样本数据

以陕西省泾惠渠灌区为研究区进行验证，灌区北依黄土塬，西至泾河，南至渭河，东至石川河，为连续分布的完整水文地质单元(见图 3)。灌区地质构造属于鄂尔多斯台向斜南缘渭河断陷盆地的中段北侧，是新生代以来的构造下陷区，为第四纪沉积物所覆盖^[19]。灌区内第四系地层比较发育，自西向东、由北向南，沉积厚度递增，岩性颗粒由粗变细。灌区浅层地下水含水层主要为分布最广泛的第四系全新统冲积层(alQ₄)，局部地区为第四系上更新统风积层(eolQ₃)，主要由沙质黏土、亚砂土、粉细砂和砂砾石组成。

选取灌区20个地下水水质分析数据作为样本数据，数据资料见表 5。采样点用GPS精准定位，采样过程按照《水质样品的保存和管理技术规定》(HJ 493-2009)和《环境水质监测质量保证手册》进行。采样点水化学类型使用Excel分析确定，利用ArcGIS的Geostatistical Analyst工具分别调用反距离加权(inverse distance weighted，IDW)、径向基函数(radial basis function，RBF)和克里金插值(Kriging)模型对采样数据进行空间插值，比对各插值模型的预测误差分布情况^{[18, 20]}，最终选择Kriging模型进行空间插值，获得6种离子毫克当量百分数的栅格数据。

2.2 结果与分析 2.2.1 插值结果

6种离子毫克当量百分数插值结果如图 4所示。阴离子沿地下水流向形成条带形分布，毫克当量百分数平均值排序为HCO₃^->SO₄^2->Cl^-，其中HCO₃^-毫克当量百分数在全区范围超过25%，是控制区域地下水化学类型的主要阴离子。阳离子分布无明显规律，毫克当量百分数平均值排序为K⁺+Na⁺>Mg²⁺>Ca²⁺，其中Ca²⁺毫克当量百分数在全区范围低于25%。

2.2.2 模型验证分析

1) 基于泰森多边形的水化学类型分区

利用ArcGIS的Create Thiessen Polygons工具对采样数据进行分析，生成水化学类型分区图如图 5所示。研究区内面积比例较大的水化学类型有4种，分别是HCO₃·SO₄·Cl-Na·Mg、SO₄·Cl·HCO₃-Na·Mg、HCO₃·SO₄-Na·Mg和SO₄·HCO₃-Na·Mg，研究区西南部阴离子以HCO₃^-为主，东北部阴离子以SO₄^2-为主，全区主要阳离子为K⁺+Na⁺、Mg²⁺。

2) 离子毫克当量百分数插值分区

根据毫克当量百分数插值分区流程，利用ArcGIS进行水化学类型分析(见图 6)。由图 6可知，研究区内有4种舒卡列夫水化学类型，分别是HCO₃·SO₄·Cl-Na·Mg、HCO₃·SO₄-Na·Mg、HCO₃·SO₄-Na和HCO₃-Na·Mg，由于舒卡列夫水化学类型未进行含量排序，无法判断主要的阴阳离子。

3) 基于空间分析模型的水化学类型分区

调用空间分析模型进行水化学类型分析，生成水化学类型分区图(见图 7)。由图 7可知，研究区内面积比例较高的水化学类型为HCO₃·SO₄·Cl-Na·Mg、SO₄·HCO₃·Cl-Na·Mg和HCO₃·SO₄-Na·Mg，研究区大部分地区阴离子以HCO_3-为主，东北部局部地区阴离子以SO_42-为主，全区主要阳离子为K⁺+Na⁺、Mg²⁺。

4) 对比分析。

分别对泰森多边形、插值分区和空间分析模型的分析结果进行统计，数据见表 6。

由表 6可知，按照面积比例，3种方法均可以确定HCO₃·SO₄·Cl-Na·Mg为主要水化学类型，其他类型存在一定的差异。泰森多边形方法中占23.13%、11.44%的SO₄·Cl·HCO₃-Na·Mg和SO₄·HCO₃-Na·Mg，在空间分析模型分析中分别为1.90%、2.59%。这是由于泰森多边形建立受采样点分布控制，当采样点密度较小时，所建立的泰森多边形单元面积较大，与水化学类型分布的实际情况存在差异。通过插值数据计算阴离子含量比较图可知(见图 8)，SO₄^2-≥HCO₃^-且HCO₃ < Cl^-区域主要集中于采样点2、8、10、11、14附近，对应SO₄·Cl·HCO₃-Na·Mg区域。由图 4可知，Cl^-毫克当量百分数小于25%区域主要集中于采样点9、19附近，对应SO₄·HCO₃-Na·Mg区域。

泰森多边形方法中占有13.08%、6.30%的HCO₃·SO₄-Na·Mg和SO₄·HCO₃·Cl-Na·Mg，在空间分析模型分析中分别为22.57%、26.48%。HCO₃·SO₄-Na·Mg分布差异主要集中在采样点12、18和20附近，SO₄·HCO₃·Cl-Na·Mg分布差异主要集中在采样点2、8、9、10、11、14附近。从空间模型角度分析，基于地质统计学的空间插值本身具有空间预测的功能，2种水化学类型空间分布是根据已有数据的空间分布特征进行的合理预测，结合图 4、图 7可知，2种水化学类型空间分布差异区域应该是水化学成分随地下水流动演化形成。根据区内阴离子含量普遍呈现HCO₃^- > SO₄^2- > Cl^-这一基本特征来比较，空间分析模型的分析结果更加客观。

与实测数据的一致性方面，分析表 5和图 7可知，空间分析模型在各采样点位置的水化学类型与实测数据分析结果完全一致。插值分区法可确定的水化学类型少于其他2种方法，且无法区分地下水中起到控制作用的阴阳离子类型，其主要原因为插值分区法无法对离子含量进行排序。

本实例分析中，3种方法的水化学类型分区结果整体趋势一致，在表达水化学类型空间分布特征方面存在一定差异。泰森多边形分析结果能真实反映各采样点实测水化学类型，基本表达了相邻样本点间的水化学类型演化，但存在相邻区域水化学类型突变现象，例如采样点2、4、7、12、14等，均出现Cl^-含量比例的异常升高。对于插值分区和空间分析模型，其结果不仅能够真实反映各采样点实测水化学类型，且可以明显演示出水化学类型的演化路径，在空间连续性与差异性的表达上优于泰森多边形方法，但插值分区法无法体现阴阳离子的含量关系。此外，从可视化效果方面来看，基于空间分析模型方法要优于泰森多边形方法，生成的图形具有较好的光滑度，图面清晰简洁，有利于进一步的分析解读。

3 讨论与分析

早期学者对运用GIS进行水化学类型分析开展了相应的研究，但过程复杂且自动程度较低，也无法获得离子含量大小排序关系^[6]。本文基于空间分析建模技术，运用ArcGIS ModelBuilder建立的水化学类型分区空间分析模型，将GIS进行水化学类型分区的流程进一步集成化处理，分析过程更加简便高效，可表达的水化学类型完整。

本文所建立的水化学类型分区空间分析模型，是以空间统计学为理论基础的空间分析方法，充分考虑了地下水化学场空间连续性、空间变异、空间结构性特征，分析结果能够更加客观地反映出地下水化学成分的空间演化特征。

水化学类型分区空间分析模型对水化学类型分区的准确性主要受空间插值数据质量影响。由于地下水中不同水化学成分的空间相关性存在较大差异，在布设采样点时，应通过历史资料进行地统计分析，获得不同水化学成分半变异函数的变程值，以最小变程值作为采样间隔，在经济成本控制范围内尽量加大采样点密度。在进行空间插值时，应充分考虑研究区内地下水化学成分的形成与分布特征，当采样点密度较大且地下水化学形成机制简单，地下水化学成分空间变异有限时，不同插值方法的结果差异较小；当地下水化学成分存在较大的空间变异特征时，全局空间自相关性较差，在建立模型时需要对其趋势特征进行探索分析，选择适当的模型进行插值。

对水化学类型分区空间分析模型的精度评价参考遥感影像分类，以已有水化学类型分区图件或采样点实测分析值为参考数据，通过ArcGIS分别统计在参考数据和空间分析模型结果图件相应位置的水化学类型，建立混淆矩阵并计算总体分类精度、生产者精度、用户精度和Kappa系数进行评价。本实例分析中，因无参考图件则以采样点作为参考数据进行精度评价，各水化学类型的总体分类精度、生产者精度、用户精度均为100%，相应Kappa系数为1.0，即空间分析模型在采样点位置处水化学类型划分无错误。

4 结语

本文针对水化学类型分区中常用的泰森多边形和毫克当量百分数插值方法特点，基于GIS空间分析建模技术，利用ArcGIS的ModelBuilder建立了地下水化学类型空间分析模型。利用实测水质分析数据对模型进行了验证，分析结果显示，地下水化学类型分区空间分析模型能够较好地进行地下水化学类型的判定与自动分区，与传统分析方法相比，更好地反映了水化学类型的空间演化特征，地图可视化效果更好。

本文模型建立在空间统计学原理基础上，采样点密度和插值方法对模型分析结果的准确性影响较大，在实际应用中需要控制数据质量，以增强模型的适用性。