基于电阻率成像法的莫高108窟墙体内水盐运移过程分析

周启友，姚炳宇，李禾澍，王彦武，陈港泉

(1.南京大学地球科学与工程学院，江苏南京 210023；2.敦煌研究院保护研究所，甘肃酒泉 736200)

0 引 言

敦煌莫高窟是中国首批全国重点文物保护单位和首批世界文化遗产，它以其精美的绘画、彩塑以及丰富的文化宝藏而闻名于世。但在长期自然环境影响下，目前莫高窟部分壁画却面临着盐害的威胁，困扰着洞窟的保护工作。这些盐害问题包括空鼓、酥碱、起甲等，其根本原因是水盐在墙体、地仗以及窟内空间运移，并在地仗层盐分浓缩结晶的结果。为此，许多学者通过多种途径探讨了洞窟墙体内的水盐运移过程。例如，闫玲等[1]通过对莫高窟第53窟内外温湿度的监测分析认为，洞窟内存在吸、放湿的调节作用，并建立了洞窟壁画的等温吸湿—放湿数理模型。王亚军等[2]也通过对敦煌莫高窟第87窟窟内温度和相对湿度的监测和分析发现，窟内温湿度特征因位置不同而有所差异，窟内上部测点平均温度高于下部测点，而上部测点相对湿度低于下部测点。李红寿等[3]更是采用封闭与开放洞窟的监测方法，推断在莫高窟北区的132窟内有来自围岩的水分进入洞窟。为了认识墙体内部的水盐运移过程，郭青林等[4]采用高密度电阻率法对莫高窟崖顶地层、窟前林带区以及降水前后崖顶电阻率的状况进行了测定和分析。杨善龙等[5]则通过挖探井取样分析说明，虽然莫高窟崖体砂砾石层中的平均含水量才1.6%，但是崖体中局部所夹的粉质黏土层却具有高达24%的含水量。郭青林等[6]通过在98窟西壁下部无壁画部位钻孔取样测定可溶盐，并在钻孔内布设高精度温湿度监测探头，在西壁下部实施由南向北的电阻率测定等研究得出结论认为，98窟西壁表面到岩体内的可溶性盐主要以硫酸盐和氯盐为主，并主要富集在表面至35 cm深的岩体内，在岩体深125 cm有凝结水产生。

显然，这些研究成果为理解洞窟墙体内的水盐运移机理提供了重要依据，但由于洞窟岩层的高度非均质性，直至目前都还未能对墙体内的水盐运移过程有清楚的把握。因此，迫切需要一种能够在三维空间上监测墙体内水盐运移过程的手段来揭示墙体内的水盐运移机理。近年来，高密度电阻率成像法以其对水盐运移过程的敏感性和动态监测的优势，在研究水盐运移方面得到了广泛的应用。例如，早在2005年，Singha等[7]就采用跨孔的高密度电阻率成像法在现场对注入盐水的扩散过程进行了成功的三维成像监测。其后，Sass等[8]应用该方法对文物保护中墙体的湿润性进行了测定，指出电阻率成像法在获得墙体的含水量分布方面具有巨大的潜力。周启友等[9]在一块砂岩上所作的饱水和排水监测实验也说明电阻率成像法在监测岩体中的水分运移时是有效的。Mol等[10]则进一步应用电阻率成像法对文物中砂岩的风化过程进行了研究，揭示了砂岩劣化的机理。Sass等[11]应用该方法对模拟降雨前后墙面上的湿润和干燥过程进行了监测，指出风化最严重的岩块表面会吸水最多也最快，但同时也干燥最快。Martinho等[12]对葡萄牙文艺复兴石质浮雕中的可溶盐分进行了电阻率成像和分布研究。Pollock[13]和Wehrer等[14]则在实验室内采用电阻率成像法监测了盐水的运移和非饱和土壤中的水分入渗与溶质传输过程。最近，Ward等[15]更是提出了根据电阻率成像法监测和卡尔曼滤波方法自动跟踪溶质迁移过程的方法。所有这些研究都充分说明高密度电阻率成像法用于对墙体中水盐运移过程的监测是完全可行的。为此，2015—2017年在莫高窟108窟后壁面上实施了1 m×1 m和2 m×2 m范围的三维高密度电阻率成像法测定，以期通过电阻率成像法的三维监测，揭示洞窟墙体内部的水盐运移过程。

1 测定与数据处理方法

针对莫高窟108窟的1 m×1 m范围的高密度电阻率成像法测定于2015年进行，地点是108窟西壁面的左下方处，该处没有地仗覆盖，可在壁面临时布设电极实现对墙体内部的测定(图1)。电极布置范围1 m2左右，最下边距离地面1.5 m，最左边距离后壁左边沿也是1.5 m。在1 m2的范围内，布置电极121个，形成11×11的电极网格，电极间距在水平和垂直方向均为0.10 m左右，电流和电位远电极分别布置在监测区域的左侧和右侧，距离地面2 m处的位置。电极采用58 mm直径的AgCl心电电极，设置时剥去心电电极上的透明胶，在需要设置电极的地方先涂上一小块黏土，再在黏土上设置电极。试验结果表明，这样设置的电极和黏土块之间是紧密接触的。但是，由于砂砾石层表面疏松或在砾石面上有许多浮土，致使黏土块和砂砾石层之间很难保持长时间的稳定接触，从而导致测试结果中出现了许多异常。

图1 108窟西壁面上高密度电阻率成像法的测定位置(图中数字为各点的坐标)Fig.1 Measurement positions for high-density electrical resistivity tomographic monitoring on the surface of the west side wall of Cave 108 (numbers in the figure mean the coordinates of the points)

对1 m×1 m范围测定结果的初步分析发现，该测定尺度对于认识墙体内的水盐运移过程偏小。于是，自2016年起覆盖原来1 m×1 m的范围，又进行了2 m×2 m范围的高密度电阻率成像法测定。电极设置数量和方法均与1 m×1 m的相同，但电极间距扩大了一倍，即0.2 m左右，测定范围约4 m2。

测试方法采用二极法进行，沿着X(以洞窟左壁面为0水平指向北侧)、Y(以洞窟西壁上边界为0垂直指向下方)、XY和-XY方向分别进行扫描测定(Z轴以壁面为0垂直壁面指向西侧)。由此实现了三维的高密度电阻值测定，每次扫描测定可获得1 700个数据。在1 m×1 m范围的电极位置和测定模式如图2所示，在2 m×2 m范围的测定与此相似。

图2 108窟西壁面1 m×1 m范围上的高密度电阻率成像法测定模式Fig.2 Measurement patterns used for high-density electrical resistivity tomography in 1 m×1 m region at the surface of the west side wall of Cave 108

测定仪器为ERT21D高密度电阻测定系统(南京九州勘探技术有限公司，南京211100)，该系统可以对128个电极实现二极法的任意组合测定，测定电压有24VDC、48VDC和168VDC可选，可以实现全自动的连续扫描测定。针对1 m×1 m范围的测定，在2015年5月上旬和8月初进行了多次。但由于8月初测定时未对电极做接触处理，所测数据未能成为有效资料，因此仅在5月上旬获得了2次有效的资料。在2 m×2 m范围的测定于2016年5月、8月和2017年8月进行，各获得了2次完整有效的资料。为了说明所得资料的可靠性，图3对比了各次测定的结果，所测数据基本上都有相同的变化规律，说明所得结果是比较可靠的。

图3 在108窟西壁面测得的电阻值的稳定性Fig.3 Stability of the obtained electrical resistance in 1 m×1 m (a) and 2 m×2 m (b) regions at the surface of the west side wall of Cave 108

基于所获得的有效资料，在剔除异常数据之后，对各测定组合的视电阻率值根据周启友[16]所介绍的方法按半无限空间介质的假设进行了计算，具体计算方法在此不再赘述。虽然洞窟建造于地层内部，整体上难以作为半无限空间介质来考虑，但考虑到电极所设置的壁面为临空平面，电极间距远远小于壁面的尺寸，半无限空间的假设在视电阻率的计算和分析中仍然适用。在本研究中，采用由此获得的视电阻率值进行不同时间测定结果的比较分析。

为了获得介质的真实电阻率分布，针对剔除异常后的数据，按半无限空间介质采用补修算法[17]进行了电阻率的反演计算。该算法应用有限元法进行计算，实际计算区域大于测定区域，并采用六面体微元进行三维剖分。对1 m2和4 m2的测定区域，剖分微元大小分别为5 cm×5 cm×5 cm和10 cm×10 cm×10 cm。反演迭代计算6次，以RMS误差降低到稳定值时的电阻率分布作为介质的真实电阻率分布来分析。

2 结果与讨论

2.1 在1 m×1 m范围内获得的电阻率空间分布特征

图4给出了2015年5月5日所获得的电阻率的空间分布图像。由于所得图像是三维的，通过平行于壁面(Z=0.00～0.20 m，Z=0.25～0.40 m，Z=0.45～0.60 m，Z=0.65～0.80 m)的电阻率平面图和垂直于壁面(Y=2.60～2.70 m，Y=2.85～2.95 m，Y=3.05～3.15 m，Y=3.30～3.40 m)的电阻率剖面图来表现(为符合显示习惯，图中的Y轴和Z轴值都取了负值，后文同样作此处理)。

由图4可见，在平行后壁面的各个深度上和垂直于壁面的各个剖面上，电阻率的分布都是不均匀的。在平面图中的上部区域(Y绝对值较小一侧)，电阻率值整体偏低；而在靠近下部的区域(Y绝对值较大一侧)，电阻率值则整体偏高，并在最浅部的一层(Z=0～0.20 m)存在高低相间的条带状的高阻或低阻区域。在垂直于壁面的剖面上，平面图中下部(见图4b中的3和4)的高阻区域向深部延伸，上部(见图4b中的1和2)的低阻值区域也在深部出现。这种在西壁面的特征(下侧高上侧低)和在测定区域下部观察到的砾石透镜体相吻合，而浅部条带状的电阻率分布特征也与构成墙体的砂砾石的层状结构相一致，说明所得的电阻率图像能够反映出墙体岩层的结构。

图4 在1 m×1 m区域内所得电阻率的空间分布(2015年5月5日)Fig.4 Spatial distributions of resistivity obtained in 1 m×1 m region (May 5，2015)

由于高电阻率区域应该对应粗粒介质的、胶结致密的或水分和盐分含量都较低的区域，而低电阻率区域则应该是细粒介质的、胶结疏松的或水分和盐分含量较高的区域。因此，所得的电阻率分布图像在反映墙体岩层结构总体格局的基础上，在同样条件下还可以反映出岩层在岩性、胶结程度、含水量、盐分等方面的差异。为此，1 m×1 m区域的上侧应该是介质颗粒细小、水分和盐分都最有可能集中的区域。

2.2 在2 m×2 m范围内获得的电阻率空间分布特征及其动态变化

图5～7分别给出了2016年5月21日、2016年8月24日和2017年8月9日在2 m×2 m范围内获得的电阻率空间分布，图中虚线区域是2015年的1 m×1 m的测定范围。平行于壁面的平面图的深度范围是Z=0.00～0.35 m，Z=0.40～0.70 m，Z=0.80～1.10 m和Z=1.20～1.50 m。垂直于壁面的剖面图在Y方向的范围是Y=2.10～2.30 m，Y=2.50～2.70 m，Y=2.90～3.10 m和Y=3.30～3.50 m。

图5 在2 m×2 m区域内所得电阻率的空间分布(2016年5月21日)Fig.5 Spatial distributions of resistivity obtained in 2 m×2 m region (May 21，2016)

由图可见，无论在平面上还是在剖面上，墙体电阻率值的空间分布在各个测定时间都是极不均匀的，都存在着范围较大的低阻条带区域，这和墙体岩层的层状结构一致。然而，尽管低阻条带区域在3次测定结果中都具有相同的空间分布结构，但它们在具体的分布范围和形状上却有着很大的不同。如果电阻率值仅仅是岩层结构的反映，那么低阻条带区域的分布在3次测定结果中应该是基本一致的，因为岩层结构不可能在短期内出现较大的变化。如果没有其他因素的影响，低阻条带区域在3次测定结果中分布的差异就只能是介质内部在水分和盐分含量上的变化所引起的。由于在3次测定间仅有较小的温度变化，且温度的变化也不可能仅引起局部电阻率的变化。低阻条带区域在分布范围和形状上的变化是墙体水分和盐分含量发生变化的结果。

基于上述推理，便可以很好地解释所得电阻率空间分布在不同时间的变化过程。例如，在2016年5月21日的图像上，有一较大的低阻区域出现在1 m×1 m测定范围的上侧，分布形状呈块状特征。而在2016年8月24日的图像上，该低阻区域却阻值增高，中心位置发生转移，分布形状更具横向特征。到2017年8月9日，该低阻区域更是在阻值的大小和分布形状上都发生了变化。2016年5月21日相对集中的低阻值区域，应该是外部水分入渗在1 m×1 m区域附近，一时性的水分集中和电阻率降低所造成的。而随后，通道内水分向周围孔隙细小的区域扩散并在那里蒸发和浓缩盐分，岩层电阻率值降低。而原来的通道内部则由于水分的减少而含水量降低，岩层电阻率值升高。正是由于这种水分的再分配、蒸发和盐分浓缩过程造成了低阻区域的转移和空间分布形态的变化。因此，不同时间获得的低阻区域在位置、范围和形状上的不同就很正常了。

非均质多孔介质中的饱和与非饱和过程完美地解释了所得电阻率的空间分布特征及其在时间上的变化过程，为认识108窟墙体内的水盐运移过程提供了重要依据。

2.3 视电阻率空间分布的时间变化过程

为了进一步探讨墙体电阻率空间分布在时间上的变化过程，直接对比了2015年5月、2016年5月、8月和2017年8月4次测定所得的视电阻率空间分布(图8)。图8中，圆点越大，颜色越接近红色，代表视电阻率值越高。可以看出，5月份在靠近洞窟壁面尚有不少高的视电阻率区域出现，而到了8月份在靠近洞窟壁面却不再有这些高的视电阻率区域，近表层完全为低的视电阻率区域所占据。假如温度、电极接触效果以及其他因素对测定结果的影响可以忽略的话，这正好说明后壁中的水分和盐分至少从5月份到8月份这一阶段有着一个由深部向浅部逐渐聚集的过程。且实际向浅部运移和聚集的时间段有可能更长，并可能与来自洞窟外部的降雨入渗过程密切相关。

2.4 洞窟墙体内的水盐运移机理

综上所述，莫高108窟墙体内的电阻率空间分布是非均质的，低阻值区域的空间位置、范围和形状也是随时间变化的。说明在构成墙体的非均质砂砾石层内存在着重要的水分饱和与非饱和过程(再分配过程)，正是这一不断重复的饱和与非饱和过程促进了墙体内水分的运移和盐分在一些特定区域的聚集。具体过程如下。

如图9所示，当墙体内因为降雨入渗等原因有较多的水分进入时，这些携带盐分的水分首先通过并集中在介质颗粒较粗、孔隙较大的区域或通道(图9a中A区)，水体流动按近似饱和流的方式进行(如入渗水量不足以达到饱和流，则水分以非饱和流的形式靠近细颗粒一侧流动)，水体优先选择大孔隙的通道并快速流动(图9a中带箭头的实线)。随着降雨入渗的减少或停止，以及由于周围颗粒较细、孔隙较小的区域的毛细吸力的作用，通道内的水分含量减少，水分向周围颗粒细小的区域扩散，水体流动主要在墙体内颗粒细小的区域上缓慢进行(图9a中B区)。当整个墙体内的水分进一步减少时，集中于颗粒细小区域的水体之间不再连续，并在三维空间上出现相互独立的现象，水分向颗粒更细的区域收缩或集中，形成独立的水分相对集中的区域(图9a中C区)，称之为水块[18]。由于颗粒细小的区域毛细吸力远大于颗粒较粗的区域，水分在此长时间滞留，面向窟内空间以及孔隙度较大的区域蒸发，水块内盐分逐渐浓缩并结晶(图9b中C区)。从这些水块上蒸发的水汽或进入窟内空间，或在孔隙较大的孔隙网络内运移，受墙体内温度和压力梯度的控制(图9b中带箭头的虚线)。当第二次降雨等引起的水分再次来临时，上述过程再次重复，盐分再次向颗粒细小的区域输送并因随后的蒸发而在那里不断浓缩，盐害由此产生。

图6 在2 m×2 m区域内所得电阻率的空间分布(2016年8月24日)Fig.6 Spatial distributions of resistivity obtained in 2 m×2 m region (Aug 24，2016)

图7 在2 m×2 m区域内所得电阻率的空间分布(2017年8月9日)Fig.7 Spatial distributions of resistivity obtained in 2 m×2 m region ( Aug 9，2017)

图8 2015年至2017年所得视电阻率空间分布的对比Fig.8 Comparison of the apparent resistivity distributions obtained from 2015 to 2017

图9 洞窟墙体内的水盐运移过程示意图Fig.9 Schematic diagram of the water and salt transport processes in the wall of the cave

上述盐分浓缩并结晶的关键是水分的间歇性补充及其随后的蒸发作用。莫高窟年降雨量不大，但降雨集中，墙体内存在长期的反复入渗作用，而戈壁气候环境又很容易促进入渗之后水分的蒸发，为此认为正是上述饱和与非饱和过程造成了洞窟壁画的盐害。由此，根据上述水盐运移机理，可以作出预测：1)洞窟盐害出现的区域应该主要分布在那些介质颗粒细小的特定区域。2)颗粒细小的区域水分含量高，盐分浓度大；反之，如粗大的砾石和孔隙度较大的区域则少。杨善龙等[5]的研究结果以及壁画盐害的实际分布区域正好验证了这一预测。

3 结 论

通过在108窟西壁面上的电极设置，成功实施了1 m×1 m和2 m×2 m范围的三维高密度电阻率成像，获得了5月份和8月份两个时间点上洞窟墙体电阻率的三维分布图像。所得结果表明，无论在1 m×1 m，还是在2 m×2 m的范围内，墙体内电阻率的空间分布都是非均匀的、动态变化的，低阻值区域在空间位置、分布范围和形状上都会随时间而变化。说明在构成墙体的非均质砂砾石层内存在着重要的水分饱和与非饱和过程(再分配过程)，正是这一不断重复的饱和与非饱和过程促进了墙体内水分的运移和盐分在一些特定区域的聚集。

致 谢：在野外工作期间，得到了敦煌研究院保护研究所郭青林和刘洲的热情帮助，在此表示感谢。