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

(1.北京洛斯達(dá)科技發(fā)展有限公司，北京 100120；2.國家基礎(chǔ)地理信息中心，北京 100830；3.中國科學(xué)院空天信息創(chuàng)新研究院，遙感科學(xué)國家重點實驗室，北京 100101)

0 引言

特高壓直流輸電接地極的選擇關(guān)系著整個接地系統(tǒng)的安全運行。由于接地極入地電流大，為滿足溫升、跨步電壓等需求，要求接地極極址的可用面積大、土壤導(dǎo)電(熱)性能好、熱容率高、表層土壤厚等[1]，因此，接地極極址土壤的物理和化學(xué)性質(zhì)均影響直流輸電接地極的接地性能[2-4]。其中，土壤電阻率是影響土壤導(dǎo)電性能的主要因素，進(jìn)而影響整個接地系統(tǒng)的安全[5]，如何高效、準(zhǔn)確地獲取土壤電阻率的空間分布對于特高壓直流輸電接地極極址的選擇至關(guān)重要。

目前，部分學(xué)者在土壤電阻率和導(dǎo)電性方面進(jìn)行了一系列研究。例如：孫宇瑞發(fā)現(xiàn)在不同土壤含水量條件下，土壤電導(dǎo)率與土壤含水量相關(guān)關(guān)系不同，含水量在15%～20%之間時，兩者呈顯著的線性關(guān)系[6]；劉春泉等研究表明，土壤含水量和氣溫是影響土壤電阻率的主要因素[7]；段旭等發(fā)現(xiàn)土壤電阻率與土壤的總孔隙度、體積含水量均具有良好的相關(guān)性[8]；Lech等研究發(fā)現(xiàn),不同土壤質(zhì)地(沙含量、黏土含量)的電阻率與土壤孔隙度之間均呈線性相關(guān)[9]；Fikri等分析了電阻率與不同土壤容重、含水量及干密度之間的相關(guān)性[10]；Poudel等發(fā)現(xiàn)土壤電阻率與土壤含水量以及土壤中氯離子和硫酸鹽離子之和具有較強的相關(guān)性[11]；李博倫等發(fā)現(xiàn)土壤中可溶鹽總量、陽離子交換量對土壤電阻率也會產(chǎn)生影響[12]。綜上可知，土壤電阻率雖受土壤含水量、孔隙度、質(zhì)地、溫度及部分理化性質(zhì)的影響，但土壤水分的導(dǎo)電性對土壤電阻率具有重要影響[13]。因此，土壤含水量實時、大范圍、高精度的估算是土壤電阻率估算的關(guān)鍵[14-16]，也是目前研究的熱點與難點[17]。

傳統(tǒng)的土壤含水量監(jiān)測方法耗時耗力，難以獲取大范圍內(nèi)土壤含水量的空間分布[18]?；谶b感技術(shù)的土壤含水量監(jiān)測主要包括基于光學(xué)遙感的監(jiān)測和基于微波遙感的監(jiān)測[19]，其中微波遙感具有穿透性強、受天氣影響小且對土壤水分敏感等優(yōu)點，更適合土壤含水量的監(jiān)測。常見的雷達(dá)土壤含水量反演方法包括變化檢測法、神經(jīng)網(wǎng)絡(luò)法、查找表法及迭代優(yōu)化法等[20]，其中變化檢測法不需要地表粗糙度等輔助信息，在實際反演工作中廣為關(guān)注，但目前鮮有結(jié)合雷達(dá)影像對土壤電阻率估算的研究。因此，本文面向直流輸電接地極設(shè)計中面臨的實際問題，通過Sentinel-1、Sentinel-2和SMAP衛(wèi)星數(shù)據(jù)聯(lián)合反演研究區(qū)范圍內(nèi)的土壤體積含水量，結(jié)合土壤容重獲取土壤重量含水量，進(jìn)而建立其與土壤電阻率之間的估算模型，最終獲取土壤電阻率的時空分布特征，為直流輸電接地極設(shè)計與優(yōu)化提供參考。

1 研究區(qū)概況

陜西省煤電資源豐富，在滿足省內(nèi)需求的同時，可大規(guī)模向外輸出；湖北省由于發(fā)展需求，需要從外省大量引入電能。陜西—湖北±800 kV特高壓直流輸電工程可解決兩省電力需求與能源分布不均衡問題，提高煤電基地的電能輸送能力和資源開發(fā)利用效率，節(jié)約電源建設(shè)和運行成本，同時帶動相關(guān)產(chǎn)業(yè)發(fā)展，促進(jìn)陜北地區(qū)經(jīng)濟(jì)的健康可持續(xù)發(fā)展。

在建設(shè)高壓直流技術(shù)工程時需要構(gòu)建相應(yīng)接地系統(tǒng)，其中接地極極址的選擇關(guān)系到接地系統(tǒng)的安全運行，陜西—湖北±800 kV特高壓直流輸電接地極的擬選極址位于陜西省榆林市境內(nèi)。本文研究區(qū)主要包括山西省的河曲縣和興縣的部分地區(qū)以及陜西省的府谷縣。研究區(qū)屬中溫帶半干旱大陸性季風(fēng)氣候，冬季平均氣溫低于0 ℃，降水量極少，地表干燥，水土流失嚴(yán)重，地形支離破碎，地貌以溝壑或山坡為主，形成特有的半干旱黃土—風(fēng)沙地貌，區(qū)內(nèi)覆蓋大面積第四系黃土；東南部地勢較高，植被覆蓋率高，西北部地表裸露，植被覆蓋率低。研究區(qū)域?qū)偃A北地層區(qū)鄂爾多斯地層小區(qū)，出露地層有三疊系、侏羅系、第三系、第四系、全新統(tǒng)，不同地層的巖性組成、地層厚度和分布區(qū)域不同。為全面分析擬選極址區(qū)域土壤電阻率的特點，在布設(shè)測量點位時需考慮研究區(qū)域的地質(zhì)概況，分多層對實地的土壤電阻率進(jìn)行測量。圖1為本文研究區(qū)域的地理位置以及電阻率實測點的具體位置。

圖1 研究區(qū)域概況Fig.1 Overview of the study area

2 實驗數(shù)據(jù)

2.1 Sentinel-1數(shù)據(jù)

Sentinel-1由同軌的兩顆衛(wèi)星組成(Sentinel-1A/B)，其搭載C波段合成孔徑雷達(dá)，具有較高的時間重返能力和空間分辨率。Sentinel-1單顆衛(wèi)星重訪時間為12 d，A/B兩顆衛(wèi)星的重訪周期為6 d。由于Sentinel-1衛(wèi)星覆蓋不均勻[21]，本文從歐空局(ESA)網(wǎng)站獲取了覆蓋全研究區(qū)的30景Sentinel-1A(Level-1 Ground Range Detected)影像數(shù)據(jù)，時間跨度為2019年1月8日-12月22日，空間分辨率為10 m。利用Sentinel-1工具箱對影像進(jìn)行熱噪聲去除、輻射定標(biāo)和地形校正等預(yù)處理。為減小地表粗糙度及雷達(dá)穿透深度對雷達(dá)回波信號的干擾，本文對預(yù)處理后的影像進(jìn)行投影轉(zhuǎn)換和重采樣(有效消除雷達(dá)數(shù)據(jù)中散斑效應(yīng)[22])，獲得WGS 84坐標(biāo)下100 m分辨率的影像數(shù)據(jù)。

為消除不同入射角的影響，需對后向散射系數(shù)進(jìn)行角度校正，基于蘭伯特定律的校正方法[23]簡便、應(yīng)用廣泛，其校正公式為：

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

(1)

式中：σref為歸一化后的后向散射系數(shù)；θref為校正的目標(biāo)角度；θi為待校正的本地入射角；σθi為待歸一化的后向散射系數(shù)；指數(shù)n一般取2。

2.2 Sentinel-2數(shù)據(jù)

Sentinel-2包含4個10 m、6個20 m和3個60 m空間分辨率波段，此外，Sentinel-2還包括一個“QA60”波段，具有云掩碼的相關(guān)信息，用于移除被云覆蓋的區(qū)域。本文獲取覆蓋研究區(qū)域的61景Sentinel-2A/B Level 1C影像數(shù)據(jù)，時間跨度為2019年1月4日-12月30日，重訪周期為5 d。利用Sen2Cor對其進(jìn)行輻射定標(biāo)和大氣校正，然后從影像數(shù)據(jù)中選取空間分辨率為10 m的波段4(Red)和波段8(NIR)計算歸一化植被指數(shù)(NDVI)。為抑制云量對NDVI的影響，本文采用三階Savitzky-Golay(S-G)濾波對原始NDVI數(shù)據(jù)進(jìn)行平滑處理，然后利用三次樣條插值技術(shù)對平滑后的NDVI進(jìn)行插值，用于計算Sentinel-1影像采集日期的NDVI[24]，最后將NDVI結(jié)果重采樣到100 m。

2.3 實測土壤電阻率數(shù)據(jù)

為分析土壤電阻率與土壤含水量的關(guān)系，本研究在榆林市境內(nèi)800 kV特高壓直流接地極工程擬選的府谷縣白云鄉(xiāng)棗林峁村極址區(qū)域，對表層深度(0～200 m)、淺層深度(0～2 000 m)以及深層深度(0～30 000 m)的視電阻率及分層電阻率進(jìn)行測量，測量時間為2017年4月25日-5月30日。對于表層深度，布設(shè)16個測點和1個質(zhì)量檢查點，采用對稱四極電測深法的溫納爾裝置進(jìn)行測量；對于淺層深度，布設(shè)9個測點和1個質(zhì)量檢查點，采用交流電勘探類的瞬變電磁測深法進(jìn)行測量；對于深層深度，布設(shè)5個測點和1個質(zhì)量檢查點，采用交流電勘探類的大地電磁測深法進(jìn)行測量。

受地表干燥程度、黃土覆蓋厚度、地下水位高度和基巖埋深厚度等多種因素影響，該極址土壤電阻率在橫向和縱向均有較大變化，實測表層土壤視電阻率變化范圍為42.1～519.4 Ω·m。

2.4 其他輔助數(shù)據(jù)

本文所用的9 km空間分辨率的SMAP土壤水分產(chǎn)品是由SMAP降尺度反演生產(chǎn)[25]，最終獲取2019年共332景SMAP增強L3輻射計日尺度土壤水分產(chǎn)品數(shù)據(jù)(SMAP enhanced L3 radiometer global daily 9 km EASE-Grid soil moisture)。6月20日-7月24日的數(shù)據(jù)缺失，但本文僅用其估算研究區(qū)域全年土壤含水量的最大值和最小值，因此可以忽略缺失數(shù)據(jù)對實驗的影響。

土壤容重數(shù)據(jù)源于中山大學(xué)陸—氣相互作用研究組開發(fā)的全球土壤數(shù)據(jù)集，該數(shù)據(jù)集基于各國家及區(qū)域土壤數(shù)據(jù)庫建立，然后采用標(biāo)準(zhǔn)化數(shù)據(jù)結(jié)構(gòu)和數(shù)據(jù)處理程序以保證數(shù)據(jù)的一致性[26]。土壤容重數(shù)據(jù)的空間分辨率為1 km，為與其他數(shù)據(jù)匹配，將其重采樣到100 m。

3 研究方法

本文研究流程(圖2)包括：1)獲取研究區(qū)Sentinel-1和Sentinel-2影像數(shù)據(jù)，并進(jìn)行預(yù)處理，將預(yù)處理后的影像數(shù)據(jù)輸入變化檢測模型中，然后利用SMAP 9 km的土壤水分產(chǎn)品數(shù)據(jù)估算土壤含水量最大值和最小值，進(jìn)而反演土壤體積含水量；2)將反演得到的土壤體積含水量結(jié)合土壤容重數(shù)據(jù)計算得到土壤重量含水量；3)利用土壤電阻率的實地測量值與遙感反演的土壤重量含水量建立土壤電阻率的估算模型；4)由估算模型對整個研究區(qū)域的土壤電阻率進(jìn)行反演，通過閾值篩選獲得適合特高壓直流輸電接地極安放區(qū)域。

3.1 土壤體積含水量計算

當(dāng)使用多時相雷達(dá)數(shù)據(jù)反演土壤含水量時，變化檢測法不需要地表粗糙度等先驗知識輔助，本文采用該方法獲取研究區(qū)域內(nèi)的土壤含水量。該方法假設(shè)后向散射系數(shù)的變化主要由土壤含水量的變化主導(dǎo)，受植被和土壤表面粗糙度的變化影響較小[27]。因此，可以通過計算不同日期雷達(dá)后向散射系數(shù)的差值消除土壤表面粗糙度和植被的影響，建立后向散射系數(shù)差值與土壤含水量差值的關(guān)系[28]。比較同一像素點所有日期Sentinel-1影像的后向散射系數(shù)，后向散射系數(shù)的最小值對應(yīng)土壤最干燥的日期。將給定日期的后向散射系數(shù)σ°(d)減去該像素后向散射系數(shù)的最小值σ°(dry)，獲取后向散射的變化值Δσ°[22](式(2))。

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

(2)

圖3 后向散射系數(shù)差值與NDVI間關(guān)系示意[22]Fig.3 Schematic diagram of the relationship between backscattering coefficient difference and NDVI

通過當(dāng)日后向散射系數(shù)差值與同一NDVI下最大后向散射差值之比，計算當(dāng)日土壤含水量差值與最大土壤含水量差值之比；然后，利用輔助信息獲取最大土壤含水量差值，進(jìn)而得到當(dāng)日土壤含水量差值；最后，結(jié)合當(dāng)日土壤含水量差值與最小土壤含水量，計算當(dāng)日土壤含水量，計算公式為：

(3)

Δmv=mvd-mvmin

(4)

Δmvmax=mvmax-mvmin

(5)

式中：Δσ°d為第d天與全年最干燥日期的后向散射系數(shù)差值；Δmv為第d天與最干燥日期的土壤含水量差值；Δmvmax為全年土壤含水量差值的最大值；mvd為第d天的土壤含水量；mvmax、mvmin分別為全年土壤含水量的最大值和最小值。通過多年水文觀測或被動微波土壤水分產(chǎn)品，可估算得到Δmvmax和mvmin，然后通過上述公式可求得mvd。

3.2 土壤重量含水量及電阻率計算

土壤重量含水量包含土壤體積含水量和土壤容重信息，常用于土壤電阻率的擬合[30]。本文利用從中山大學(xué)陸—氣相互作用研究組獲取的土壤容重數(shù)據(jù)bd以及反演的土壤體積含水量mv計算土壤重量含水量Mmv(式(6))，然后利用土壤重量含水量與實測土壤電阻率進(jìn)行建模，最終獲取研究區(qū)域土壤電阻率的空間分布。

Mmv=mv/bd

(6)

4 結(jié)果與分析

4.1 土壤體積含水量和重量含水量空間分布

由原始NDVI與預(yù)處理后NDVI數(shù)據(jù)的時間序列(圖4)可以看出，原始NDVI在1-2月離散程度較高，可能是由于1-2月植被覆蓋度較低，衛(wèi)星接收到的信號受地面噪聲影響較大。經(jīng)過濾波、插值預(yù)處理后的NDVI值既能保留原有NDVI的波動信息，又能實現(xiàn)對NDVI空值日期數(shù)據(jù)的補充。因此，預(yù)處理后的NDVI值可滿足構(gòu)建土壤含水量反演模型的需求。

圖4 Sentinel-2計算的原始NDVI和NDVI插值數(shù)據(jù)的時間序列Fig.4 Time series of the original NDVI and interpolated NDVI data based on Sentinel-2

基于插值后的NDVI數(shù)據(jù)與Sentinel-1數(shù)據(jù)繪制后向散射差值與NDVI之間的散點圖(圖5)，以獲取f(NDVI)。圖5具有明顯的分段特征，當(dāng)NDVI<0.12時，曲線擬合結(jié)果為f(NDVI)= 580.41×NDVI2-128.29×NDVI+18.07,當(dāng)NDVI≥0.12時，f(NDVI)= 1.3×NDVI3-19.8×NDVI2+15.8×NDVI+9.41。同時，利用SMAP 9 km分辨率的土壤水分產(chǎn)品估算得到Δmvmax為0.2 cm3/cm3，mvmin為0.02 cm3/cm3。根據(jù)式(2)-式(5)可反演得到研究區(qū)域的土壤體積含水量，對獲取的土壤體積含水量按月平均，得到研究區(qū)每月以及全年的土壤體積含水量空間分布(圖6)。

圖5 研究區(qū)后向散射系數(shù)差值與NDVI之間關(guān)系Fig.5 Relationship between backscattering coefficient difference and NDVI in the study area

從圖6可以看出，研究區(qū)域內(nèi)土壤體積含水量呈現(xiàn)東南高、西北低的趨勢，原因可能是東南地區(qū)植被覆蓋度較高，土壤蒸發(fā)少[31]。研究區(qū)域土壤體積含水量整體偏低，大部分地區(qū)土壤體積含水量在0.04～0.14 cm3/cm3之間，5月土壤體積含水量高于全年均值。

圖6 研究區(qū)土壤體積含水量分布Fig.6 Distribution of soil volumetric water content in the study area

為驗證土壤含水量的反演精度，本文將Sentinel-1數(shù)據(jù)反演的土壤含水量產(chǎn)品(分辨率為100 m)重采樣到9 km，然后計算每月的平均土壤含水量，與SMAP土壤含水量產(chǎn)品的月均值進(jìn)行比較。如圖7所示，Sentinel-1與SMAP土壤含水量產(chǎn)品均小于0.2 cm3/cm3，在土壤含水量較高的情況下，Sentinel-1反演的土壤含水量略低于SMAP土壤含水量，二者總體上具有較好的一致性(R=0.731,ubRMSD=0.025 cm3/cm3)。

注：虛線指±0.1 cm3/cm3的邊距。

為分析研究區(qū)全年土壤含水量的變化趨勢，本文利用16個采樣點對應(yīng)位置土壤含水量的平均值繪制折線圖(圖8)，可以看出，采樣點全年的土壤體積含水量變化不大，土壤含水量處于較低狀態(tài)。

圖8 待選極址區(qū)域土壤含水量時間變化Fig.8 Time change of soil water content in the electrode site area to be selected

土壤電阻率與土壤重量含水量具有較好的相關(guān)性[14,32]，利用式(6)計算得到研究區(qū)域的土壤重量含水量空間分布(圖9)。對比圖6與圖9可以看出，土壤重量含水量與土壤體積含水量空間分布大體相似，在整個研究區(qū)內(nèi)也呈現(xiàn)出東南高、西北低的趨勢。由于土壤重量含水量包含土壤容重的部分信息，在空間分布上包含的內(nèi)容更豐富。研究區(qū)東南出現(xiàn)部分像素土壤重量含水量明顯高于周圍地區(qū)的情況，主要是受土壤容重數(shù)據(jù)的影響，異常區(qū)域土壤容重明顯小于周圍地區(qū)。由圖8可知，土壤重量含水量較低且變化不大，其變化趨勢與土壤體積含水量相近。

圖9 研究區(qū)土壤重量含水量分布Fig.9 Distribution of soil gravimetric water content in the study area

4.2 土壤電阻率空間分布

本文土壤電阻率實測時間為2017年4-5月，由于2017年光學(xué)影像受云覆蓋影響，對土壤含水量反演影響較大，最終獲取的土壤重量含水量為2019年4-5月的平均值。雖然土壤電阻率與土壤重量含水量之間存在時間差異，但由于研究區(qū)域土壤重量含水量整體較低，2017年與2019年土壤重量含水量變化不明顯。同時，本文利用土壤重量含水量的月均值，進(jìn)一步減小了不同年份之間土壤重量含水量的變化差異。通過分析2017年與2019年4-5月SMAP土壤含水量之間的相關(guān)性(圖10)可知，2017年和2019年4-5月的SMAP土壤含水量均值具有較好的相關(guān)性，R為0.825，ubRMSD為0.008 cm3/cm3。由于土壤電阻率主要受土壤類型、土壤含水量以及土壤成分等因素影響，短期內(nèi)土壤的類型、成分不會發(fā)生劇烈變化，因此，在土壤含水量變化不明顯的基礎(chǔ)上，可利用2019年的Sentinel-1與Sentinel-2獲取土壤重量含水量，進(jìn)而與土壤實測電阻率進(jìn)行建模。此外，月均值可減小遙感的瞬時特性所引起的代表性較差問題。

圖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)與土壤電阻率關(guān)系的研究，鑒于此，本文對0.03～0.07 g/g范圍內(nèi)的土壤重量含水量與土壤電阻率之間的關(guān)系進(jìn)行建模(圖11)，進(jìn)而獲取研究區(qū)內(nèi)土壤電阻率的空間分布(圖12)。從圖11可以看出，當(dāng)土壤重量含水量較低時，土壤電阻率隨土壤含水量變化迅速，當(dāng)土壤重量含水量偏高時，土壤電阻率隨土壤含水量變化減緩，但整體上表現(xiàn)出土壤重量含水量越低土壤電阻率越高的趨勢。

圖11 土壤電阻率與重量含水量擬合關(guān)系Fig.11 Fitting relationship between soil resistivity and gravimetric water content

圖12 土壤電阻率分布Fig.12 Distribution of soil resistivity

從圖12可以看出，土壤電阻率分布與土壤重量含水量分布具有較好的一致性。特高壓直流輸電工程極址一般選擇在電阻率小于100 Ω·m的區(qū)域[33]，但考慮到研究區(qū)整體電阻率偏高，所以本文以200 Ω·m為閾值對研究區(qū)域內(nèi)的電阻率進(jìn)行劃分(圖13)。陜西—湖北±800 kV特高壓直流輸電工程極址要從陜西省選擇，山西省不在考慮范圍內(nèi)。從圖13可以看出，擬選的棗林峁村極址區(qū)域位于土壤電阻率較低的區(qū)域，極址選擇合理。所以，通過遙感的方法可以為特高壓直流輸電接地極極址選擇提供數(shù)據(jù)支持，具體極址的選擇還需要綜合考慮地形地貌、現(xiàn)有交通、電力線路、油氣管道等因素。

圖13 接地極極址待選區(qū)域Fig.13 Areas to be selected for grounding electrode site

5 結(jié)論

土壤電阻率是影響接地系統(tǒng)設(shè)計的重要因素，而土壤含水量是影響土壤電阻率的關(guān)鍵因子。本文通過Sentinel-1和Sentinel-2衛(wèi)星相結(jié)合的方式，使用基于時序遙感觀測的方法穩(wěn)定獲取土壤體積含水量信息，然后結(jié)合土壤容重等輔助信息獲取土壤重量含水量，為及時、準(zhǔn)確地估計土壤電阻率提供了一種新的方法和途徑。研究獲取的高分辨率土壤水分信息與SMAP粗分辨率土壤水分產(chǎn)品具有較好的一致性，二者的相關(guān)系數(shù)R為0.731,ubRMSD為0.025 cm3/cm3。

本文方法結(jié)合雷達(dá)與光學(xué)影像，在獲取土壤電阻率的空間分布方面具有獨特優(yōu)勢，但由于利用遙感數(shù)據(jù)通常僅能獲取土壤的體積含水量，在進(jìn)一步獲取土壤電阻率時仍然依賴于經(jīng)驗性的統(tǒng)計關(guān)系。這種統(tǒng)計關(guān)系后續(xù)需要在不同植被類型、土壤質(zhì)地區(qū)域進(jìn)行檢驗或者修正，以完善土壤電阻率與土壤重量含水量之間的轉(zhuǎn)換模型，提升本文方法在不同地理區(qū)域的適用性和魯棒性。