土壤含水率與GPR信號(hào)屬性關(guān)系數(shù)值模擬研究

伊天宇

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

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

在物理模擬方面，目前已有利用GPR探測(cè)土壤含水率的實(shí)例。馬福建等[2]利用GPR探測(cè)人工模型得到了反射系數(shù)與含水率之間的定量表達(dá)式。吳信民[3]提出相位移時(shí)間的概念，并建立了土體含水率與相位移時(shí)間關(guān)系的理論模型。而在數(shù)值模擬方面，方慧[1]定性的研究了介質(zhì)含水率與雷達(dá)信號(hào)的關(guān)系，其中用到了雷達(dá)振幅屬性?？偟膩?lái)說，運(yùn)用GPR探測(cè)土壤含水率定性研究居多，而定量研究相對(duì)較少。

基于此，本文采用隨機(jī)過程理論建立土壤含水平穩(wěn)隨機(jī)介質(zhì)模型，并根據(jù)雷達(dá)振幅和含水率差值來(lái)討論雷達(dá)信號(hào)屬性與土壤含水率之間的關(guān)系。振幅是雷達(dá)信號(hào)中絕對(duì)意義的量，含水率差值是相對(duì)意義上的量。只要能夠得到淺層土壤含水率，即可通過關(guān)系式推導(dǎo)下一層土壤含水率。

1 基本理論

1.1 水的電性特征

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

實(shí)驗(yàn)表明介電常數(shù)還隨頻率變化，用Debye公式來(lái)描述變化規(guī)律：

式中：εS、ε∞——直流和極高頻下的介質(zhì)介電常數(shù)；

frel——弛豫頻率。

Debye模型表明：在低頻和高頻狀態(tài)下，介質(zhì)的介電常數(shù)近似為實(shí)數(shù)，頻散不明顯，能量耗散較小。

1.2 Topp公式

本文主要基于Topp[4]通過實(shí)驗(yàn)給出的土壤介電常數(shù)與土壤含水率之間的關(guān)系來(lái)建立模型：

式中：εb——土壤整體介電常數(shù)；

θ——土壤含水率。

2 方法

2.1 建模思路

含水土壤中，水的存在將會(huì)影響土壤的電學(xué)特性，即會(huì)影響土壤的介電常數(shù)和電導(dǎo)率。由于水隨機(jī)分布在土壤中，該固液-雙相介質(zhì)屬于隨機(jī)介質(zhì)，所以根據(jù)式及隨機(jī)介質(zhì)理論建立含水土壤平穩(wěn)隨機(jī)介質(zhì)模型。具體實(shí)現(xiàn)如下：

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

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

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

表1 土壤含水率與相對(duì)介電常數(shù)對(duì)應(yīng)數(shù)據(jù)

（4）實(shí)現(xiàn)層內(nèi)含水率隨機(jī)變化[5]，逐層求取y方向上每層網(wǎng)格內(nèi)的含水率：

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

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

δf(x,z)——隨機(jī)擾動(dòng)的含水率；

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

（5）產(chǎn)生隨機(jī)變化的介電常數(shù)體，將模型每個(gè)位置的含水率代入式。即可求得該位置的介電常數(shù)，從而構(gòu)建含水土壤體介電模型。

2.2 含水土壤體的算法實(shí)現(xiàn)

一般情況下隨機(jī)過程不存在傅里葉變換，但可從自相關(guān)函數(shù)的功率譜來(lái)描述，其算法原理如下：

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

（2）對(duì)自相關(guān)函數(shù)做二維傅氏變換，得到φ(kx+kz)；

（3）用隨機(jī)數(shù)發(fā)生器產(chǎn)生在[0,2π]服從均勻分布的二維隨機(jī)場(chǎng)Θ(kx,kz)；

（4）計(jì)算隨機(jī)功率譜函數(shù)：

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

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

（7）再由式得到模型。y方向單層網(wǎng)格點(diǎn)隨機(jī)擾動(dòng)建立結(jié)果如圖2（a）所示，含水率結(jié)果如圖2（b）所示（第一層含水率15.58%，第二層含水率20.23%），相對(duì)介電常數(shù)結(jié)果如圖2（c）所示。

2.3 正演模擬參數(shù)設(shè)置

探測(cè)方式為剖面法，天線間距為20cm，以5cm為探測(cè)間距沿測(cè)線移動(dòng)，測(cè)線選取為模型中線。采用三維時(shí)域有限差分法對(duì)含水土體隨機(jī)介質(zhì)進(jìn)行正演模擬[7]。

為滿足穩(wěn)定性條件，通常剖分單元每一個(gè)邊的長(zhǎng)度不得大于發(fā)射天線主頻波長(zhǎng)的10，并且要求時(shí)間步長(zhǎng)的設(shè)計(jì)滿足下式：

式中：Δx、Δy和Δz——剖分單元的長(zhǎng)度；

c——光速。

考慮天線主頻，選取空間網(wǎng)格dx=0.001m，時(shí)間步長(zhǎng)dt由式計(jì)算來(lái)得到，以滿足穩(wěn)定性關(guān)系。

震源采用500MHz雷克子波。

2.4 數(shù)據(jù)處理

2.4.1 最優(yōu)頻率選取

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

2.4.2 層速度求取

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

因而有對(duì)應(yīng)關(guān)系見表2。

2.4.3 振幅求取

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

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

其中H表示Hilbert變換，?表示卷積。x(t)與其Hilbert變換是正交的。Hilbert變換相當(dāng)于通過一個(gè)沖擊響應(yīng)為的線性網(wǎng)絡(luò)，即Hilbert濾波器。

對(duì)于x(t)，若其Hilbert變換為R(t)，則有瞬時(shí)包絡(luò)：

圖4為其道Ex分量包絡(luò)圖。

表2 土壤含水率與層速度對(duì)應(yīng)數(shù)據(jù)

3 結(jié)果與分析

3.1 土壤含水率與層速度

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

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

式中：θ——含水率；

v——層速度，ms。

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

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

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

表3 含水率差值與Ex振幅對(duì)應(yīng)數(shù)據(jù)

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

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

A——振幅。

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

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

4 結(jié)論

介質(zhì)含水率變化對(duì)雷達(dá)信號(hào)的影響在多數(shù)情況下被認(rèn)為是干擾，由于水具有特殊性質(zhì)，使電磁波在介質(zhì)的傳播過程更加復(fù)雜。通過GPR信號(hào)屬性研究土壤含水率，重點(diǎn)探討了界面反射振幅與土壤含水率差值可能存在的關(guān)系。界面反射振幅與土壤含水率差值都與土壤含水相鄰兩層相關(guān)，都屬于相對(duì)意義上的量，在實(shí)際中有更好的應(yīng)用前景。

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

但在實(shí)際中，要想應(yīng)用GPR信號(hào)屬性探測(cè)土壤含水率，土壤中必須有一個(gè)明顯的標(biāo)志反射層，如果土壤呈近乎均勻分布，探測(cè)土壤含水率的工作也難以為繼。正常情況下，土壤含水率越往深越大。但在正演模擬中，若下一層土壤含水率比上一層土壤含水率低，雷達(dá)信號(hào)便會(huì)出現(xiàn)極性反轉(zhuǎn)的情況，這表明GPR的其他信號(hào)屬性也有可能較好的描述土壤含水率。

總的來(lái)說，應(yīng)用雷達(dá)信號(hào)屬性測(cè)量土壤含水率是可行的，有著廣闊的應(yīng)用前景。其無(wú)損、非接觸、高效的特性表現(xiàn)出其它檢測(cè)介質(zhì)含水率的傳統(tǒng)技術(shù)所不具備的優(yōu)勢(shì)。