伊天宇
(中國(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)下一層土壤含水率。
在時(shí)間域中,電磁波在介質(zhì)中的傳播速度v由因子k決定[1]:
實(shí)驗(yàn)表明介電常數(shù)還隨頻率變化,用Debye公式來(lái)描述變化規(guī)律:
式中:εS、ε∞——直流和極高頻下的介質(zhì)介電常數(shù);
frel——弛豫頻率。
Debye模型表明:在低頻和高頻狀態(tài)下,介質(zhì)的介電常數(shù)近似為實(shí)數(shù),頻散不明顯,能量耗散較小。
本文主要基于Topp[4]通過實(shí)驗(yàn)給出的土壤介電常數(shù)與土壤含水率之間的關(guān)系來(lái)建立模型:
式中:εb——土壤整體介電常數(shù);
θ——土壤含水率。
含水土壤中,水的存在將會(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)建含水土壤體介電模型。
一般情況下隨機(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)所示。
探測(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.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ù)
由前面的表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)確的得到表層土壤的含水率。
分別將含水率設(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)系逐層向下求取含水率差值,與表層含水率依次相加便可得到不同深度的含水率值。
介質(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ì)。