不同充填基質復耕土壤含水率反演與特征分析

王 陽，周 濤，程 琦，宋桂芳,張世文*

(1.安徽理工大學 地球與環境學院，安徽 淮南 232001；2.安徽理工大學 空間信息與測繪工程學院，安徽 淮南 232001；3.安徽理工大學 安徽省高潛水位礦區水土資源綜合利用與生態保護工程實驗室，安徽 淮南 232001)

隨著我國城鎮化步伐的加快，耕地后備資源不斷減少，人地矛盾愈演愈烈，加之人類不合理的耕地利用方式，耕地質量日趨下降，耕地供需平衡壓力日益加大，土地復耕已成為落實耕地占補平衡的一項有力舉措。其中，土壤含水率是判斷復耕、復種效果指標體系中的重要一環[1]。快速、精準、無損探測復耕土壤含水率對復耕農業水資源管理、復耕效果監測以及復耕農田水文研究具有重要意義。目前，含水率測量主要使用烘干法，雖然取得的結果準確性較高，但烘干法效率低下，大范圍探測較為困難。而中子射線法和TDR等方法操作復雜且需要進行事先標定[2]。衛星遙感等空間信息處理技術測量含水率的范圍較大(通常是對大于1公頃的地塊)，對于構造較為復雜、充填基質多樣的新增復耕耕地，無法滿足復耕質量檢測精度的要求。探地雷達作為一種新型的近地微波物探技術，逐漸運用于資源環境監測，對于含水率的探測較為快速和準確，既彌補了烘干法在效率方面的不足，也能夠對充填多樣、層次復雜的土壤含水率產生響應。

近年來，探地雷達在土壤含水率探測方面的研究進展較快，同時分析方法已逐漸從波速、距離估計法預測向電磁波屬性特征分析轉變。Wang等[3]利用探地雷達對青藏高原季節性凍土區進行含水率探測，結合振幅分析、波形分析以及電磁波衰減屬性等分析方法，研究得出振幅屬性以及瞬時頻率等屬性與含水率間具有一定的相關性。Shen等[4]通過分析電磁波不同特征屬性與有機碳、含水率的響應關系，研究得出所有特征屬性都與土壤含水率呈現顯著相關，其中振幅分析中能量屬性與含水率相關性最高。Fabio等[5]通過室內實驗分析探地雷達譜峰頻率分量與土壤中黏土含量的響應關系，實驗結果表明電磁波譜峰頻率分量與土壤中黏土含量具有較好的響應關系，隨著黏土含量的增加，頻率峰值逐漸減小。Ferrara等[6]在自然條件(野外環境)下通過研究電磁波信號的早期振幅包絡平均值與不同深度土壤含水率的響應關系，結果表明早期振幅包絡平均值能夠準確反演出0～50 cm范圍內的土壤含水率。羅古拜[7]利用GPR反射波法探測結果以及烘干法結果進行建模，建立了介電常數和含水率的模型。喬新濤等[8]采用包絡振幅平均值對黃土高原復墾礦區土壤進行探測，對黃土區土壤水分反演精度較高。程琦等[9]使用AEA法對以粉煤灰為充填基質的復墾土壤進行探測，發現電磁波特征波段與耕層含水率存在較強響應，反演擬合效果較好。目前研究大多基于自然土壤下的含水率探測，而對于不同充填基質下復耕重構土壤的探測鮮少。復耕重構土壤因受到人為干預導致其土壤環境較正常土壤更為復雜，內部充填基質不再均一，雷達信號受影響因素頗多。而新構土壤含水量的研究對于掌握礦區土壤水土保持生態修復十分關鍵。郭婷婷[10]以礦區重構土壤含水率為研究對象，結果表明不同粒徑的粉煤灰、煤矸石呈現不同的土壤水分下滲規律，從而影響土壤田間持水量。胡振琪等[11]在自然條件下開展不同夾心土層位置的黃河泥沙夾層式充填復墾，結果表明土壤水分累計入滲量隨夾層位置的變化而變化，HYDRUSHY-1D結合含水量變化反演土壤水力參數精度較高。不同的土壤重構、充填方式對于土壤含水率保持至關重要。

因此，為探究不同充填基質土壤環境下探地雷達含水率探測問題，本文提出了一種利用零點調節、一維濾波、道間均衡、背景消除、小波去噪、FIR濾波、Hilbert變換等信號處理流程。通過構建復耕耕地不同充填基質下土壤含水率與振幅屬性的關系模型，并開展野外試驗探究模型適用深度以及不同充填基質下土壤含水率差異規律。然后選取最優地統計插值方法對復耕區水分空間變異性進行了分析。研究結果可在復耕農業水資源管理、復耕效果監測以及復耕農田水文研究等方面提供理論與技術支撐。

1 實驗設計與方法原理

1.1 實驗設計

實驗開展在安徽省淮南市安徽理工大學試驗區，屬亞熱帶季風性氣候，土壤質地為黏土，土壤容重偏大，介于1.3～1.6 g/cm3。實驗時間為6～7月，實驗選擇在最后一次降雨的3 d后開展。設置2個1 m×1 m的實驗樣方，為還原真實復耕耕地地下充填物的復雜與多樣性，填埋物采用粉煤灰、建筑用紅磚、瓷磚等。A樣方以瓷磚和紅磚分別充填；B樣方以粉煤灰和紅磚碎料分別充填，覆土厚度均為40 cm。具體充填方式、方位、雷達測線及取樣點如圖1所示。雷達天線中心頻率為500 MHz-Shield，時窗設置為50ns。由于充填層厚度要大于最小垂直分辨率(500 MHz中心頻率換算后垂直分辨率約為5 cm)，故充填基質厚度統一設置為20 cm厚。充填完畢后，等待土層的自然沉降，于一個月后對試驗區進行雷達探測。雷達探測后重新開挖坡面，在每種充填模式下延測線方向探測4個層次的土壤容重和含水率，分別為0～10 cm、10～20 cm、20～30 cm和30～40 cm。為保證數據的準確性，當天進行實驗室烘干和測量記錄。

1.2 原理與方法

(1)雷達探測原理。探地雷達(Ground Penetrating，GPR)是對地下物體發射高頻電磁波(1 MHz～10 GHz)的一種物探方法。由于電磁波在傳播過程中遇到不同介質時會發生波的折射和反射，當以耕層土壤為研究對象時，會因為不同的充填基質土壤介質產生不同的反射信號。

電磁波速度與介電常數的關系式為

(1)

本實驗電磁波垂直入射，表達式為

(2)

(2)數據預處理。在得到每個剖面的探測圖像后，先用Reflexw軟件對每個圖像都進行零點調節、一維濾波、道間均衡、背景消除、小波去噪等方法處理。零點調節是對雷達剖面圖像進行上下移動，移動到大氣與地面的分界面，切除地上部分，保留地下部分。一維濾波的作用主要是壓制干擾信號提高信噪比，提取地下介質的響應特征信號。道間均衡通過道間的相關性進行信號的加強，對于道間不相關的信號，處理后信號減弱，此方法對雷達剖面圖像中相關性較好的弱信號有極大的增強作用。背景消除是為了凸顯異常信號，消除背景噪聲和由于阻抗不匹配產生的駐波干擾，達到壓縮水平干擾信號，提取有用信號的目的。小波去噪也是用于壓制噪聲信號，提高有用信號。水分反演部分預處理主要是利用Matlab軟件，發現第一個正半周期是信號與土壤介電特性相關性最高的部分，能夠使電磁波信號的信噪比最大化，并使淺層界面反射的干擾最小化。因此只需要進行去直流漂移處理與時間零點校正，其中去直流漂移目的是使有效信號不受漂移現象的影響。

(3)FIR(Finite Impulse Response)濾波。FIR濾波系統函數為H(z)：

(3)

(4)

h(n)=ω(n)×hd(n)，

(5)

根據時域相乘關系與頻域卷積關系得

(6)

式中,H(ejω)為FIR濾波器頻率響應；WR[ej(ω-θ)]為窗譜。

(4)提取GPR信號振幅包絡值。探地雷達通過500 MHz-Shield天線對耕層發射電磁波。土壤中含水率的變化是影響雷達反射波振幅的主要原因。利用Hilbert變換后解析信號得到信號振幅包絡值。公式如下：

(7)

解析信號為

R(t)=x(t)+iH(t),

(8)

則其振幅包絡為

(9)

(5)土壤干容重與體積含水率測定。為驗證AEA法反演含水率的效果，在每個樣方同一深度，利用環刀均勻提取3份樣品，裝入密封袋中。提取當天在實驗室將樣品放置在烘干機內，在105 ℃環境下烘干24 h，烘干后得到土壤干重，取平均值。然后利用以下公式計算出干容重和體積含水率：

(10)

(11)

θv=θm×ρ,

(12)

式中,θv、θm分別為土壤樣品的體積含水率和質量含水率；mwet和mdry分別是土壤濕重和土壤干重；v是環刀體積，為100 cm3；ρ為土壤樣品的容重，單位為 g/cm3。

2 結果與分析

2.1 基于AEA法的含水率反演模型優選

為了探究雷達電磁波與含水率的響應關系，提取不同測線單道波形，截取特定時窗下雷達信號波，經預處理后，求得振幅包絡值，并將測線處附近各測點的振幅包絡值取平均值。隨機選取測線上各個層次20個樣點為建模集，10個樣點為驗證集，建立多種回歸模型如表1所示。從建模集的R2來看，不同擬合模型的精度差異較小，對數擬合模型精度最高，R2為0.939；線性擬合模型精度最低，為0.889。多項式擬合模型RMSE最低，為1.987。從驗證集的R2來看，擬合效果均在0.89以上，說明基于0～12 ns的波段提取的信號振幅包絡平均值與土壤含水率有較強的響應關系，通過AEA法對含水率進行反演是較為精準的。從驗證集RMSE來看，對數擬合模型RMSE最低，為1.707。土壤含水率反演模型優選如圖3所示。圖3a是振幅包絡平均值與實測含水率的4種回歸擬合，圖3b是含水率實測值與模型預測值對比圖，其中橫軸表示10個驗證集點，縱軸表示含水率。其中，對數擬合模型的預測效果最佳。

表1 AEA-1含水率建模精度

2.2 基于AEA法的含水率反演

為探究AEA法對不同層次土壤含水率的識別能力，使用以上確定的對數模型分別對試驗區土壤剖面0～10 cm、10～20 cm、20～30 cm和30～40 cm各層的含水率進行反演，結果如圖4所示。由圖4可知，從R2來看，0～10 cm土壤含水率反演精度最高，為0.89；30～40 cm反演精度最低，為0.65。以RMSE作為反演精度標準，0～10 cm土壤含水率反演精度最高，為1.82；20～30 cm土壤含水率反演精度最低，為3.50。在20～30 cm和30～40 cm反演中，均存在部分差異過大點位。基于0～4 ns和4.2～8.1 ns的波段提取的信號振幅包絡平均值與0～10 cm和10～20 cm土壤含水率有較強的響應關系，通過AEA法對含水率進行反演是可行且較為精準的。8.1～11.6 ns和11.6～15.4 ns波段提取的信號振幅包絡平均值與20～30 cm和30～40 cm土壤含水率響應關系較差。

為更直觀體現不同層次含水率的差異及反演模型的驗證精度，以A樣方為例，以對數回歸模型作為矩形區域雷達數據反演模型，分別進行各個層次土壤實測值與AEA預測值的插值，其中圖示顏色越深表明含水率越高。應用Sufer 14軟件進行克里金網格插值得到該實驗樣方各個層次的含水率平面分布圖，如圖5所示。圖5中分別為烘干法實測含水率分布圖和雷達信號反演含水率分布圖。根據插值結果，兩者含水率分布較為接近，反演效果較好，但是存在部分差異較大現象，其中20～30 cm和30～40 cm兩者之間差異性相較于上兩層較大。

2.3 不同充填基質下水平、垂直向含水率規律分析

(1)不同充填基質垂直向土壤含水率規律。為探究基于AEA法的不同充填基質下垂直向水分演變規律，于兩個樣方各個充填基質上布設測線，以確定的反演模型進行反演并揭示規律。運用Origin軟件繪制不同基質測線中3道雷達電磁波所反演的含水率如圖6所示。從圖6a中可以看出，在沉降穩定后，不同充填基質下土壤各層次含水率均遵循隨著深度的增加，含水率上升這一趨勢。從圖6b正態曲線來看，0～10 cm的隨機變量正態分布較為密集，σ較小；30～40 cm的隨機變量正態分布較為疏松，σ較大。各種充填模式下的土壤含水率遵循深層的離散程度比淺層離散程度大。從A樣方與B樣方相同的充填基質紅磚可以看出，由于沉降的作用，B樣方的含水率要高于A樣方的含水率。通過對比不同充填基質下不同深度土壤含水率，以粉煤灰為充填基質的土壤含水率在各層均大于其他充填基質，其他不同充填基質下的含水率差異較小。

(2)不同充填基質水平向土壤含水率規律。為探究基于AEA法的不同充填基質下水平向水分演變規律，于兩個樣方各個充填基質上布設成散射形由內向外測線。提取測線上雷達信號電磁波段，以確定的反演模型進行反演并揭示規律。并將反演的結果以插值的形式表現，不同充填基質下表層土壤含水率分布圖如圖7所示，顏色越深表示含水率越高。從圖7中可以看出，不同充填模式下，樣方周圍的土壤含水率大體上遵循距離樣方越遠，含水率越低的規律。其中，從A樣方的瓷磚與紅磚充填可以看出，以紅磚為充填基質的土壤含水率要高于瓷磚；從B樣方的紅磚充填與粉煤灰充填可以看出，以粉煤灰為充填基質的土壤含水率高于紅磚。從測線向土壤含水率變化規律來看，以瓷磚作為充填基質的土壤含水率除充填位置中心外，與周邊土壤含水率無明顯差異。以紅磚為充填基質的土壤含水率比周邊土壤含水率高，延測線方向，呈現先減少后增加，最后與周邊土壤含水率無差異。以粉煤灰為充填基質的土壤含水率比周邊土壤含水率高，越往外含水率越少。通過含水率等值線分布可以看出，A樣方以紅磚為充填基質的土壤含水率等值線范圍大約為40 cm；B樣方以粉煤灰為充填基質的土壤含水率等值線范圍大約為50 cm，以紅磚為充填基質的土壤含水率等值線范圍與A樣方紅磚相近，約為30 cm。

3 討論

探地雷達通過主機及天線收發電磁波，通過截取時窗、時間、頻率對土壤進行無損探測。在本次研究中發現500 MHz雷達電磁波信號特征時窗波段與土壤體積含水率具有較強的相關性，通過振幅包絡值與含水率的相關性擬合可以建立反演模型。通過隨機選取的建模集與驗證集，分別建立多種回歸模型，對數擬合模型精度最高，R2為0.939。運用以上確定的模型對不同深度土壤含水率進行反演，基于0～4 ns和4.2～8.1 ns的波段提取的信號振幅包絡平均值與0～10 cm和10～20 cm土壤含水率有較強的響應關系。從正態曲線疏松程度σ來看，30～40 cm的隨機變量正態分布較0～10 cm大，0～10 cm反演精度較好，通過AEA法對含水率進行反演是可行且較為精準的，R2分別為0.89和0.80；8.1～11.6 ns和11.6～15.4 ns波段提取的信號振幅包絡平均值與20～30 cm和30～40 cm土壤含水率響應關系較差，RMSE均在3.20以上，與張金珠等[12]研究結果一致。

在更深層的土壤含水率反演中，AEA法的適用性依舊需要進行進一步的研究。隨著土壤深度的加深，該法對含水率的反演精度逐漸減弱。在先前的研究中發現，影響AEA方法探測深度的深淺主要有兩方面因素，一方面是與探地雷達的主頻信號有關，而這種情況可以通過頻譜補償來提高最終的反演精度。比如程琦等[13]利用CZT法進行頻譜細化處理，根據增強后的頻譜信息進行含水率的反演，結果表明，反演精度有了大幅度的提升；而另一方面，時窗的選取劃分也會對含水率的反演精度造成影響。由于在不同深度的反演過程中，在時窗劃分中存在一定的誤差，無法根據實測厚度進行準確的時窗拾取，而隨著誤差的傳遞，會引起時窗選取的偏差不斷加大，使得最終的精度大幅度下降，所以在實際探測中，若需要獲取較為精確的反演效果，可以在探測時取一些點位進行標定處理；或者可以采用不同方法對含水率進行聯合反演，比如崔凡等[14]根據不同深度的土壤分別采用AEA法和AMRA法進行含水率的聯合反演，從而獲取較高的精度。

本文利用AEA法對不同充填類型的土壤含水率進行了反演并揭示規律。實驗結果表明，不同充填基質下土壤各層次含水率均遵循隨著深度的增加，含水率上升這一趨勢，與崔凡等[15]和張世文等[16]的研究結果一致。樣方周圍的土壤含水率大體上遵循距離樣方越遠，含水率越低的規律，與Cui等[17]的研究結果一致。以粉煤灰為充填基質的土壤含水率比周邊土壤含水率高，越往外含水率逐漸減少，與Tripathi等[18-19]的研究成果一致。在本次實驗中，A樣方以紅磚為充填基質的土壤含水率等值線范圍大約為40 cm；B樣方以粉煤灰為充填基質的土壤含水率等值線范圍大約為50 cm，以紅磚為充填基質的土壤含水率等值線范圍與A樣方紅磚相近，約為30 cm。由此可以看出粉煤灰充填的土壤含水率等值線范圍大于紅磚和瓷磚，Algeo等[20]的研究成果可作為佐證。張敬凱等[21]對采煤塌陷地帶干旱期和豐水期的塌陷區和非塌陷區土壤含水率進行研究，與本實驗B樣方含水率呈現的規律較為一致，夏季雨水較為充沛，水分入滲比周邊無充填基質土壤多，故從雷達測線上呈現距離樣方越近含水率越高現象。該法對于不同充填基質的覆土土壤均有較好的適用性，雖含水率的分布有所不同，仍得到了較高的精度。

4 結論

本文通過GPR對不同充填基質各層土壤進行探測，得出以下結論：

(1)通過對雷達數據預處理，進行Hilbert變換，提取特征波段，求得振幅包絡平均值，建立與含水率的響應關系。分別建立兩個實驗樣方不同充填類型下的不同層次土壤與AEA-1的關系模型。通過烘干法的實測值與AEA-1預測值對比，RMSE均在2.01以下。結果表明AEA法適用于不同充填類型的復雜土壤環境，這為通過GPR快速無損監測新增耕地土壤含水率提供了一種可行途徑。

(2)對比不同充填基質各層次AEA法和實測法所得到的結果發現，無論何種充填基質，AEA法反演的0～10 cm、10～20 cm含水率與實測值R2均大于0.80，但在20～30 cm和30～40 cm含水率與實測值R2均小于0.80，對于較深層土壤含水率反演建模仍需進一步研究。

(3)不同充填基質下各層土壤含水率均遵循隨著深度增加，含水率上升這一趨勢；以粉煤灰為充填基質的土壤含水率大于紅磚和瓷磚。以樣方中心沉降處為起點向外延伸的測線反演結果顯示，以紅磚、粉煤灰為充填基質的土壤含水率比周邊無充填基質的土壤高。通過測線含水率反演分析等值線圖，紅磚為充填基質的土壤含水率等值線范圍大約為40 cm，B樣方以粉煤灰為充填基質的土壤含水率等值線范圍大約為50 cm。