面向高压直流输电接地极设计的土壤电阻率遥感估算

刘 斌，王 振，谭 荣 荣，李 志 斌，赵 天 杰

(1.北京洛斯达科技发展有限公司，北京 100120；2.国家基础地理信息中心，北京 100830；3.中国科学院空天信息创新研究院，遥感科学国家重点实验室，北京 100101)

0 引言

特高压直流输电接地极的选择关系着整个接地系统的安全运行。由于接地极入地电流大，为满足温升、跨步电压等需求，要求接地极极址的可用面积大、土壤导电(热)性能好、热容率高、表层土壤厚等[1]，因此，接地极极址土壤的物理和化学性质均影响直流输电接地极的接地性能[2-4]。其中，土壤电阻率是影响土壤导电性能的主要因素，进而影响整个接地系统的安全[5]，如何高效、准确地获取土壤电阻率的空间分布对于特高压直流输电接地极极址的选择至关重要。

目前，部分学者在土壤电阻率和导电性方面进行了一系列研究。例如：孙宇瑞发现在不同土壤含水量条件下，土壤电导率与土壤含水量相关关系不同，含水量在15%～20%之间时，两者呈显著的线性关系[6]；刘春泉等研究表明，土壤含水量和气温是影响土壤电阻率的主要因素[7]；段旭等发现土壤电阻率与土壤的总孔隙度、体积含水量均具有良好的相关性[8]；Lech等研究发现,不同土壤质地(沙含量、黏土含量)的电阻率与土壤孔隙度之间均呈线性相关[9]；Fikri等分析了电阻率与不同土壤容重、含水量及干密度之间的相关性[10]；Poudel等发现土壤电阻率与土壤含水量以及土壤中氯离子和硫酸盐离子之和具有较强的相关性[11]；李博伦等发现土壤中可溶盐总量、阳离子交换量对土壤电阻率也会产生影响[12]。综上可知，土壤电阻率虽受土壤含水量、孔隙度、质地、温度及部分理化性质的影响，但土壤水分的导电性对土壤电阻率具有重要影响[13]。因此，土壤含水量实时、大范围、高精度的估算是土壤电阻率估算的关键[14-16]，也是目前研究的热点与难点[17]。

传统的土壤含水量监测方法耗时耗力，难以获取大范围内土壤含水量的空间分布[18]。基于遥感技术的土壤含水量监测主要包括基于光学遥感的监测和基于微波遥感的监测[19]，其中微波遥感具有穿透性强、受天气影响小且对土壤水分敏感等优点，更适合土壤含水量的监测。常见的雷达土壤含水量反演方法包括变化检测法、神经网络法、查找表法及迭代优化法等[20]，其中变化检测法不需要地表粗糙度等辅助信息，在实际反演工作中广为关注，但目前鲜有结合雷达影像对土壤电阻率估算的研究。因此，本文面向直流输电接地极设计中面临的实际问题，通过Sentinel-1、Sentinel-2和SMAP卫星数据联合反演研究区范围内的土壤体积含水量，结合土壤容重获取土壤重量含水量，进而建立其与土壤电阻率之间的估算模型，最终获取土壤电阻率的时空分布特征，为直流输电接地极设计与优化提供参考。

1 研究区概况

陕西省煤电资源丰富，在满足省内需求的同时，可大规模向外输出；湖北省由于发展需求，需要从外省大量引入电能。陕西—湖北±800 kV特高压直流输电工程可解决两省电力需求与能源分布不均衡问题，提高煤电基地的电能输送能力和资源开发利用效率，节约电源建设和运行成本，同时带动相关产业发展，促进陕北地区经济的健康可持续发展。

在建设高压直流技术工程时需要构建相应接地系统，其中接地极极址的选择关系到接地系统的安全运行，陕西—湖北±800 kV特高压直流输电接地极的拟选极址位于陕西省榆林市境内。本文研究区主要包括山西省的河曲县和兴县的部分地区以及陕西省的府谷县。研究区属中温带半干旱大陆性季风气候，冬季平均气温低于0 ℃，降水量极少，地表干燥，水土流失严重，地形支离破碎，地貌以沟壑或山坡为主，形成特有的半干旱黄土—风沙地貌，区内覆盖大面积第四系黄土；东南部地势较高，植被覆盖率高，西北部地表裸露，植被覆盖率低。研究区域属华北地层区鄂尔多斯地层小区，出露地层有三叠系、侏罗系、第三系、第四系、全新统，不同地层的岩性组成、地层厚度和分布区域不同。为全面分析拟选极址区域土壤电阻率的特点，在布设测量点位时需考虑研究区域的地质概况，分多层对实地的土壤电阻率进行测量。图1为本文研究区域的地理位置以及电阻率实测点的具体位置。

图1 研究区域概况Fig.1 Overview of the study area

2 实验数据

2.1 Sentinel-1数据

Sentinel-1由同轨的两颗卫星组成(Sentinel-1A/B)，其搭载C波段合成孔径雷达，具有较高的时间重返能力和空间分辨率。Sentinel-1单颗卫星重访时间为12 d，A/B两颗卫星的重访周期为6 d。由于Sentinel-1卫星覆盖不均匀[21]，本文从欧空局(ESA)网站获取了覆盖全研究区的30景Sentinel-1A(Level-1 Ground Range Detected)影像数据，时间跨度为2019年1月8日-12月22日，空间分辨率为10 m。利用Sentinel-1工具箱对影像进行热噪声去除、辐射定标和地形校正等预处理。为减小地表粗糙度及雷达穿透深度对雷达回波信号的干扰，本文对预处理后的影像进行投影转换和重采样(有效消除雷达数据中散斑效应[22])，获得WGS 84坐标下100 m分辨率的影像数据。

为消除不同入射角的影响，需对后向散射系数进行角度校正，基于兰伯特定律的校正方法[23]简便、应用广泛，其校正公式为：

σref=σθicosn(θref)/cosn(θi)

(1)

式中：σref为归一化后的后向散射系数；θref为校正的目标角度；θi为待校正的本地入射角；σθi为待归一化的后向散射系数；指数n一般取2。

2.2 Sentinel-2数据

Sentinel-2包含4个10 m、6个20 m和3个60 m空间分辨率波段，此外，Sentinel-2还包括一个“QA60”波段，具有云掩码的相关信息，用于移除被云覆盖的区域。本文获取覆盖研究区域的61景Sentinel-2A/B Level 1C影像数据，时间跨度为2019年1月4日-12月30日，重访周期为5 d。利用Sen2Cor对其进行辐射定标和大气校正，然后从影像数据中选取空间分辨率为10 m的波段4(Red)和波段8(NIR)计算归一化植被指数(NDVI)。为抑制云量对NDVI的影响，本文采用三阶Savitzky-Golay(S-G)滤波对原始NDVI数据进行平滑处理，然后利用三次样条插值技术对平滑后的NDVI进行插值，用于计算Sentinel-1影像采集日期的NDVI[24]，最后将NDVI结果重采样到100 m。

2.3 实测土壤电阻率数据

为分析土壤电阻率与土壤含水量的关系，本研究在榆林市境内800 kV特高压直流接地极工程拟选的府谷县白云乡枣林峁村极址区域，对表层深度(0～200 m)、浅层深度(0～2 000 m)以及深层深度(0～30 000 m)的视电阻率及分层电阻率进行测量，测量时间为2017年4月25日-5月30日。对于表层深度，布设16个测点和1个质量检查点，采用对称四极电测深法的温纳尔装置进行测量；对于浅层深度，布设9个测点和1个质量检查点，采用交流电勘探类的瞬变电磁测深法进行测量；对于深层深度，布设5个测点和1个质量检查点，采用交流电勘探类的大地电磁测深法进行测量。

受地表干燥程度、黄土覆盖厚度、地下水位高度和基岩埋深厚度等多种因素影响，该极址土壤电阻率在横向和纵向均有较大变化，实测表层土壤视电阻率变化范围为42.1～519.4 Ω·m。

2.4 其他辅助数据

本文所用的9 km空间分辨率的SMAP土壤水分产品是由SMAP降尺度反演生产[25]，最终获取2019年共332景SMAP增强L3辐射计日尺度土壤水分产品数据(SMAP enhanced L3 radiometer global daily 9 km EASE-Grid soil moisture)。6月20日-7月24日的数据缺失，但本文仅用其估算研究区域全年土壤含水量的最大值和最小值，因此可以忽略缺失数据对实验的影响。

土壤容重数据源于中山大学陆—气相互作用研究组开发的全球土壤数据集，该数据集基于各国家及区域土壤数据库建立，然后采用标准化数据结构和数据处理程序以保证数据的一致性[26]。土壤容重数据的空间分辨率为1 km，为与其他数据匹配，将其重采样到100 m。

3 研究方法

本文研究流程(图2)包括：1)获取研究区Sentinel-1和Sentinel-2影像数据，并进行预处理，将预处理后的影像数据输入变化检测模型中，然后利用SMAP 9 km的土壤水分产品数据估算土壤含水量最大值和最小值，进而反演土壤体积含水量；2)将反演得到的土壤体积含水量结合土壤容重数据计算得到土壤重量含水量；3)利用土壤电阻率的实地测量值与遥感反演的土壤重量含水量建立土壤电阻率的估算模型；4)由估算模型对整个研究区域的土壤电阻率进行反演，通过阈值筛选获得适合特高压直流输电接地极安放区域。

3.1 土壤体积含水量计算

当使用多时相雷达数据反演土壤含水量时，变化检测法不需要地表粗糙度等先验知识辅助，本文采用该方法获取研究区域内的土壤含水量。该方法假设后向散射系数的变化主要由土壤含水量的变化主导，受植被和土壤表面粗糙度的变化影响较小[27]。因此，可以通过计算不同日期雷达后向散射系数的差值消除土壤表面粗糙度和植被的影响，建立后向散射系数差值与土壤含水量差值的关系[28]。比较同一像素点所有日期Sentinel-1影像的后向散射系数，后向散射系数的最小值对应土壤最干燥的日期。将给定日期的后向散射系数σ°(d)减去该像素后向散射系数的最小值σ°(dry)，获取后向散射的变化值Δσ°[22](式(2))。

Δσ°=σ°(d)-σ°(dry)

(2)

图3 后向散射系数差值与NDVI间关系示意[22]Fig.3 Schematic diagram of the relationship between backscattering coefficient difference and NDVI

通过当日后向散射系数差值与同一NDVI下最大后向散射差值之比，计算当日土壤含水量差值与最大土壤含水量差值之比；然后，利用辅助信息获取最大土壤含水量差值，进而得到当日土壤含水量差值；最后，结合当日土壤含水量差值与最小土壤含水量，计算当日土壤含水量，计算公式为：

(3)

Δmv=mvd-mvmin

(4)

Δmvmax=mvmax-mvmin

(5)

式中：Δσ°d为第d天与全年最干燥日期的后向散射系数差值；Δmv为第d天与最干燥日期的土壤含水量差值；Δmvmax为全年土壤含水量差值的最大值；mvd为第d天的土壤含水量；mvmax、mvmin分别为全年土壤含水量的最大值和最小值。通过多年水文观测或被动微波土壤水分产品，可估算得到Δmvmax和mvmin，然后通过上述公式可求得mvd。

3.2 土壤重量含水量及电阻率计算

土壤重量含水量包含土壤体积含水量和土壤容重信息，常用于土壤电阻率的拟合[30]。本文利用从中山大学陆—气相互作用研究组获取的土壤容重数据bd以及反演的土壤体积含水量mv计算土壤重量含水量Mmv(式(6))，然后利用土壤重量含水量与实测土壤电阻率进行建模，最终获取研究区域土壤电阻率的空间分布。

Mmv=mv/bd

(6)

4 结果与分析

4.1 土壤体积含水量和重量含水量空间分布

由原始NDVI与预处理后NDVI数据的时间序列(图4)可以看出，原始NDVI在1-2月离散程度较高，可能是由于1-2月植被覆盖度较低，卫星接收到的信号受地面噪声影响较大。经过滤波、插值预处理后的NDVI值既能保留原有NDVI的波动信息，又能实现对NDVI空值日期数据的补充。因此，预处理后的NDVI值可满足构建土壤含水量反演模型的需求。

图4 Sentinel-2计算的原始NDVI和NDVI插值数据的时间序列Fig.4 Time series of the original NDVI and interpolated NDVI data based on Sentinel-2

基于插值后的NDVI数据与Sentinel-1数据绘制后向散射差值与NDVI之间的散点图(图5)，以获取f(NDVI)。图5具有明显的分段特征，当NDVI<0.12时，曲线拟合结果为f(NDVI)= 580.41×NDVI2-128.29×NDVI+18.07,当NDVI≥0.12时，f(NDVI)= 1.3×NDVI3-19.8×NDVI2+15.8×NDVI+9.41。同时，利用SMAP 9 km分辨率的土壤水分产品估算得到Δmvmax为0.2 cm3/cm3，mvmin为0.02 cm3/cm3。根据式(2)-式(5)可反演得到研究区域的土壤体积含水量，对获取的土壤体积含水量按月平均，得到研究区每月以及全年的土壤体积含水量空间分布(图6)。

图5 研究区后向散射系数差值与NDVI之间关系Fig.5 Relationship between backscattering coefficient difference and NDVI in the study area

从图6可以看出，研究区域内土壤体积含水量呈现东南高、西北低的趋势，原因可能是东南地区植被覆盖度较高，土壤蒸发少[31]。研究区域土壤体积含水量整体偏低，大部分地区土壤体积含水量在0.04～0.14 cm3/cm3之间，5月土壤体积含水量高于全年均值。

图6 研究区土壤体积含水量分布Fig.6 Distribution of soil volumetric water content in the study area

为验证土壤含水量的反演精度，本文将Sentinel-1数据反演的土壤含水量产品(分辨率为100 m)重采样到9 km，然后计算每月的平均土壤含水量，与SMAP土壤含水量产品的月均值进行比较。如图7所示，Sentinel-1与SMAP土壤含水量产品均小于0.2 cm3/cm3，在土壤含水量较高的情况下，Sentinel-1反演的土壤含水量略低于SMAP土壤含水量，二者总体上具有较好的一致性(R=0.731,ubRMSD=0.025 cm3/cm3)。

注：虚线指±0.1 cm3/cm3的边距。

为分析研究区全年土壤含水量的变化趋势，本文利用16个采样点对应位置土壤含水量的平均值绘制折线图(图8)，可以看出，采样点全年的土壤体积含水量变化不大，土壤含水量处于较低状态。

图8 待选极址区域土壤含水量时间变化Fig.8 Time change of soil water content in the electrode site area to be selected

土壤电阻率与土壤重量含水量具有较好的相关性[14,32]，利用式(6)计算得到研究区域的土壤重量含水量空间分布(图9)。对比图6与图9可以看出，土壤重量含水量与土壤体积含水量空间分布大体相似，在整个研究区内也呈现出东南高、西北低的趋势。由于土壤重量含水量包含土壤容重的部分信息，在空间分布上包含的内容更丰富。研究区东南出现部分像素土壤重量含水量明显高于周围地区的情况，主要是受土壤容重数据的影响，异常区域土壤容重明显小于周围地区。由图8可知，土壤重量含水量较低且变化不大，其变化趋势与土壤体积含水量相近。

图9 研究区土壤重量含水量分布Fig.9 Distribution of soil gravimetric water content in the study area

4.2 土壤电阻率空间分布

本文土壤电阻率实测时间为2017年4-5月，由于2017年光学影像受云覆盖影响，对土壤含水量反演影响较大，最终获取的土壤重量含水量为2019年4-5月的平均值。虽然土壤电阻率与土壤重量含水量之间存在时间差异，但由于研究区域土壤重量含水量整体较低，2017年与2019年土壤重量含水量变化不明显。同时，本文利用土壤重量含水量的月均值，进一步减小了不同年份之间土壤重量含水量的变化差异。通过分析2017年与2019年4-5月SMAP土壤含水量之间的相关性(图10)可知，2017年和2019年4-5月的SMAP土壤含水量均值具有较好的相关性，R为0.825，ubRMSD为0.008 cm3/cm3。由于土壤电阻率主要受土壤类型、土壤含水量以及土壤成分等因素影响，短期内土壤的类型、成分不会发生剧烈变化，因此，在土壤含水量变化不明显的基础上，可利用2019年的Sentinel-1与Sentinel-2获取土壤重量含水量，进而与土壤实测电阻率进行建模。此外，月均值可减小遥感的瞬时特性所引起的代表性较差问题。

图10 SMAP 2017年和 2019年4-5月平均土壤含水量的比较Fig.10 Comparison of average soil water content of SMAP in April and May between 2017 and 2019

目前鲜有开展较低土壤重量含水量(<0.1 g/g)与土壤电阻率关系的研究，鉴于此，本文对0.03～0.07 g/g范围内的土壤重量含水量与土壤电阻率之间的关系进行建模(图11)，进而获取研究区内土壤电阻率的空间分布(图12)。从图11可以看出，当土壤重量含水量较低时，土壤电阻率随土壤含水量变化迅速，当土壤重量含水量偏高时，土壤电阻率随土壤含水量变化减缓，但整体上表现出土壤重量含水量越低土壤电阻率越高的趋势。

图11 土壤电阻率与重量含水量拟合关系Fig.11 Fitting relationship between soil resistivity and gravimetric water content

图12 土壤电阻率分布Fig.12 Distribution of soil resistivity

从图12可以看出，土壤电阻率分布与土壤重量含水量分布具有较好的一致性。特高压直流输电工程极址一般选择在电阻率小于100 Ω·m的区域[33]，但考虑到研究区整体电阻率偏高，所以本文以200 Ω·m为阈值对研究区域内的电阻率进行划分(图13)。陕西—湖北±800 kV特高压直流输电工程极址要从陕西省选择，山西省不在考虑范围内。从图13可以看出，拟选的枣林峁村极址区域位于土壤电阻率较低的区域，极址选择合理。所以，通过遥感的方法可以为特高压直流输电接地极极址选择提供数据支持，具体极址的选择还需要综合考虑地形地貌、现有交通、电力线路、油气管道等因素。

图13 接地极极址待选区域Fig.13 Areas to be selected for grounding electrode site

5 结论

土壤电阻率是影响接地系统设计的重要因素，而土壤含水量是影响土壤电阻率的关键因子。本文通过Sentinel-1和Sentinel-2卫星相结合的方式，使用基于时序遥感观测的方法稳定获取土壤体积含水量信息，然后结合土壤容重等辅助信息获取土壤重量含水量，为及时、准确地估计土壤电阻率提供了一种新的方法和途径。研究获取的高分辨率土壤水分信息与SMAP粗分辨率土壤水分产品具有较好的一致性，二者的相关系数R为0.731,ubRMSD为0.025 cm3/cm3。

本文方法结合雷达与光学影像，在获取土壤电阻率的空间分布方面具有独特优势，但由于利用遥感数据通常仅能获取土壤的体积含水量，在进一步获取土壤电阻率时仍然依赖于经验性的统计关系。这种统计关系后续需要在不同植被类型、土壤质地区域进行检验或者修正，以完善土壤电阻率与土壤重量含水量之间的转换模型，提升本文方法在不同地理区域的适用性和鲁棒性。