土壤含水率與GPR信號屬性關系數值模擬研究

伊天宇

(中國海洋大學海洋地球科學學院，山東青島266100)

實際生活中，應用時域反射法檢測物質含水量已較為普遍。目前，探地雷達正演模擬通常假設目標體及其周圍介質是均勻的，然而物質天然存在著不均勻性，雷達波在不均勻介質中的傳播特性會發生改變。特別是，介質中的水會影響整體介質的介電常數。因此，應用探地雷達技術探測物質含水情況具有良好的地球物理前提[1]。

在物理模擬方面，目前已有利用GPR探測土壤含水率的實例。馬福建等[2]利用GPR探測人工模型得到了反射系數與含水率之間的定量表達式。吳信民[3]提出相位移時間的概念，并建立了土體含水率與相位移時間關系的理論模型。而在數值模擬方面，方慧[1]定性的研究了介質含水率與雷達信號的關系，其中用到了雷達振幅屬性?？偟膩碚f，運用GPR探測土壤含水率定性研究居多，而定量研究相對較少。

基于此，本文采用隨機過程理論建立土壤含水平穩隨機介質模型，并根據雷達振幅和含水率差值來討論雷達信號屬性與土壤含水率之間的關系。振幅是雷達信號中絕對意義的量，含水率差值是相對意義上的量。只要能夠得到淺層土壤含水率，即可通過關系式推導下一層土壤含水率。

1 基本理論

1.1 水的電性特征

在時間域中，電磁波在介質中的傳播速度v由因子k決定[1]：

實驗表明介電常數還隨頻率變化，用Debye公式來描述變化規律：

式中：εS、ε∞——直流和極高頻下的介質介電常數；

frel——弛豫頻率。

Debye模型表明：在低頻和高頻狀態下，介質的介電常數近似為實數，頻散不明顯，能量耗散較小。

1.2 Topp公式

本文主要基于Topp[4]通過實驗給出的土壤介電常數與土壤含水率之間的關系來建立模型：

式中：εb——土壤整體介電常數；

θ——土壤含水率。

2 方法

2.1 建模思路

含水土壤中，水的存在將會影響土壤的電學特性，即會影響土壤的介電常數和電導率。由于水隨機分布在土壤中，該固液-雙相介質屬于隨機介質，所以根據式及隨機介質理論建立含水土壤平穩隨機介質模型。具體實現如下：

（1）含水土壤體規模：0.5m×1m×0.5m，如圖1所示；

（2）模型分層如圖1所示；

（3）每層賦予背景含水率，共設計8種含水率并代入式。計算，如表1所示；

表1 土壤含水率與相對介電常數對應數據

（4）實現層內含水率隨機變化[5]，逐層求取y方向上每層網格內的含水率：

式中：f(x,z)——含水率；

f0——每層的背景含水率；

δf(x,z)——隨機擾動的含水率；

σ2(x,z)——具有一定均值和方差的二階平穩的空間隨機過程。

（5）產生隨機變化的介電常數體，將模型每個位置的含水率代入式。即可求得該位置的介電常數，從而構建含水土壤體介電模型。

2.2 含水土壤體的算法實現

一般情況下隨機過程不存在傅里葉變換，但可從自相關函數的功率譜來描述，其算法原理如下：

（1）選擇自相關函數φ(x,z)，這里選用高斯型橢圓自相關函數[6]，其表達式為：

（2）對自相關函數做二維傅氏變換，得到φ(kx+kz)；

（3）用隨機數發生器產生在[0,2π]服從均勻分布的二維隨機場Θ(kx,kz)；

（4）計算隨機功率譜函數：

（5）計算R(kx,kz)的二維反傅氏變換得到=(kx,kz)；

（6） 計 算 均 值μ=E[,z)]和 方 差d2=E[((x,z)-μ)2]并規范化，即使μ=0,d2≈0.03；

（7）再由式得到模型。y方向單層網格點隨機擾動建立結果如圖2（a）所示，含水率結果如圖2（b）所示（第一層含水率15.58%，第二層含水率20.23%），相對介電常數結果如圖2（c）所示。

2.3 正演模擬參數設置

探測方式為剖面法，天線間距為20cm，以5cm為探測間距沿測線移動，測線選取為模型中線。采用三維時域有限差分法對含水土體隨機介質進行正演模擬[7]。

為滿足穩定性條件，通常剖分單元每一個邊的長度不得大于發射天線主頻波長的10，并且要求時間步長的設計滿足下式：

式中：Δx、Δy和Δz——剖分單元的長度；

c——光速。

考慮天線主頻，選取空間網格dx=0.001m，時間步長dt由式計算來得到，以滿足穩定性關系。

震源采用500MHz雷克子波。

2.4 數據處理

2.4.1 最優頻率選取

Topp公式在高頻域(500～1000MHz)輕質土壤含水量與土壤介電常數的關系擬合中效果最好[3]，故本文選取天線頻率為500MHz，如圖3（截取部分圖像）所示。

2.4.2 層速度求取

依次將各含水率（共8種）設置在第一層，第二層設置為與第一層含水率不同。根據第一層與第二層界面反射波時至與震源時間差值t（即雙程走時）及圖像上每層深度h與雷達天線距d來求取層速度：

因而有對應關系見表2。

2.4.3 振幅求取

本文研究的振幅處理方式為Hilbert變換[8]，該變換是研究瞬時包絡、瞬時相位及瞬時頻率的重要手段。

對于任意連續實函數x(t)的Hilbert變換定義為：

其中H表示Hilbert變換，?表示卷積。x(t)與其Hilbert變換是正交的。Hilbert變換相當于通過一個沖擊響應為的線性網絡，即Hilbert濾波器。

對于x(t)，若其Hilbert變換為R(t)，則有瞬時包絡：

圖4為其道Ex分量包絡圖。

表2 土壤含水率與層速度對應數據

3 結果與分析

3.1 土壤含水率與層速度

由前面的表2可知，層速度對土壤含水率非常敏感，即介電常數依舊是對土壤含水率的指示性指標，利用MATLAB將表2數據進行回歸分析可得層速度與土壤含水率的關系模型，可以看出呈現明顯的負相關線性關系，擬合曲線見圖5。

回歸分析選用線性模型，R2=0.9927，擬合非常好，表達式為：

式中：θ——含水率；

v——層速度，ms。

由圖5可以看出土壤含水率與層速度存在線性模型關系，含水率越大，層速度越小。理論上，如果地下存在明顯反射層，應用寬角法等測取表層土壤層速度，便可很方便準確的得到表層土壤的含水率。

3.2 土壤含水率差值與初至反射波振幅

分別將含水率設置在第一二層，保證下一層比上一層的含水率要大，類似于圖3，可得到一二層分界面的反射波，再進行Hilbert變換求取振幅信息，這里選用Ex分量進行處理，處理得到的含水率差值與反射波振幅的數據，見表3。利用MATLAB軟件進行回歸分析，擬合曲線見圖6。

表3 含水率差值與Ex振幅對應數據

回歸分析選用線性模型，R2=0.9939，擬合很好，并由擬合參數得出表達式為：

式中：θ1、θ2——相鄰2個土壤層的含水率，%；

A——振幅。

由圖6可以看出含水率差值與反射振幅存在明顯的正相關線性關系。隨著含水率差值的增加，界面反射振幅也會跟著增加。反射振幅具有相對意義。由表3可知，反射振幅對含水率差值非常敏感。

在實際野外土壤中，水分是連續分布的，只要地下有明顯的分界面，可以先利用寬角法等[9]對分界面上覆土壤進行速度測量，再代入式便可得到分界面上覆土壤的整體含水率。要繼續遞推往下的土壤含水率，需要從GPR雷達記錄中提取反射序列，明確反射界面。由于重力的因素，含水率往往是越往深越大，故可依據反射振幅與含水率差值的關系逐層向下求取含水率差值，與表層含水率依次相加便可得到不同深度的含水率值。

4 結論

介質含水率變化對雷達信號的影響在多數情況下被認為是干擾，由于水具有特殊性質，使電磁波在介質的傳播過程更加復雜。通過GPR信號屬性研究土壤含水率，重點探討了界面反射振幅與土壤含水率差值可能存在的關系。界面反射振幅與土壤含水率差值都與土壤含水相鄰兩層相關，都屬于相對意義上的量，在實際中有更好的應用前景。

應用界面反射振幅與土壤含水率的關系式，只要得到表層的土壤含水率，依照反射系數分辨清楚反射振幅，并應用合適的振幅補償，就能將層深引起的振幅衰減降低，從而突出含水率差值對反射振幅的影響，進而逐層換算含水率差值，由淺到深，逐層計算，進而得到深層土壤的含水層。

但在實際中，要想應用GPR信號屬性探測土壤含水率，土壤中必須有一個明顯的標志反射層，如果土壤呈近乎均勻分布，探測土壤含水率的工作也難以為繼。正常情況下，土壤含水率越往深越大。但在正演模擬中，若下一層土壤含水率比上一層土壤含水率低，雷達信號便會出現極性反轉的情況，這表明GPR的其他信號屬性也有可能較好的描述土壤含水率。

總的來說，應用雷達信號屬性測量土壤含水率是可行的，有著廣闊的應用前景。其無損、非接觸、高效的特性表現出其它檢測介質含水率的傳統技術所不具備的優勢。