面向高壓直流輸電接地極設計的土壤電阻率遙感估算

劉 斌，王 振，譚 榮 榮，李 志 斌，趙 天 杰

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

本文方法結合雷達與光學影像，在獲取土壤電阻率的空間分布方面具有獨特優勢，但由于利用遙感數據通常僅能獲取土壤的體積含水量，在進一步獲取土壤電阻率時仍然依賴于經驗性的統計關系。這種統計關系后續需要在不同植被類型、土壤質地區域進行檢驗或者修正，以完善土壤電阻率與土壤重量含水量之間的轉換模型，提升本文方法在不同地理區域的適用性和魯棒性。