探地雷達樁底橢圓周測線3D wiggle剖面生成方法研究

李 好

(1．中煤科工集團 重慶研究院有限公司，重慶 400039；2．瓦斯災(zāi)害監(jiān)控與應(yīng)急技術(shù)國家重點實驗室，重慶 400037)

0 引言

各類人工挖孔橢圓樁在鐵路、市政等工程建設(shè)中得到了廣泛地應(yīng)用，尤其是改革開放四十多年以來，大量普通鐵路、高速鐵路工程的設(shè)計、建造過程中，由于鐵路橋梁的樁基礎(chǔ)承受軌道上機車的垂直應(yīng)力往往遠遠大于公路橋梁的樁基礎(chǔ)，為了防止樁身的彎矩過大，常需要把鐵路橋梁的人工挖孔灌注樁的橫截面設(shè)計為橢圓形，同時加之人工挖孔灌注樁本身的諸多技術(shù)優(yōu)點[1-3]，所以人工挖孔橢圓灌注樁基礎(chǔ)在各類鐵路工程橋梁工程中了廣泛地應(yīng)用；另外，橢圓樁在高層建筑、市政工程其他工程領(lǐng)域也得了廣泛地應(yīng)用[4-5]。但是在建造過程中樁底開挖至設(shè)計標高時，需要對基巖完整性做進一步判別，尤其是基巖為灰?guī)r時，更需要準確判別樁底基巖完整性，探測查明是否具有足夠厚度的底板厚度，一定范圍內(nèi)是存在空洞或裂隙發(fā)育等情況。目前對于樁底基巖完整性檢測的方法有鉆探法、地震波反射法、探地雷達法、高密度電法等幾種主要方法[6-10]，實際工作中，探地雷達法作為必選方法得到了廣泛應(yīng)用，鉆探法作為探地雷達法、地震波反射法圈定的異常的驗證方法，通常是物探類方法和鉆探類方法結(jié)合使用，其中的探地雷達法測線布置方式一般為1條環(huán)測線、多條直測線的測網(wǎng)模式[11-12]；其樁底為圓周測線的情況，李好[13]已解決其可視化科學(xué)計算及正確顯示問題。但是隨著國內(nèi)、外大力發(fā)展高鐵等交通基礎(chǔ)設(shè)施，有數(shù)量驚人的挖孔橢圓樁的樁底在進行探地雷達周測線探測工作，故解決關(guān)于更復(fù)雜的橢圓周測線的3D wiggle剖面可視化科學(xué)計算及正確顯示問題，具有較強的工程實際意義，同時與傳統(tǒng)的2D wiggle剖面相比，探地雷達樁底橢圓周測線3D wiggle剖面本身屬于非直線型物探測線剖面，只有進行3D可視化顯示才能準確表述其測線、測道、各采樣點的空間分布位置，才便于以圖準確表述相應(yīng)的異常區(qū)及準確表述異常位于樁身哪一面或具體什么方位等重要信息，但是傳統(tǒng)的2D wiggle剖面只能準確到距測線多少米處，具體位于樁的哪一面或哪個方位無法直接獲得，另外3D wiggle剖面圖還可以通過任意方向的圖形旋轉(zhuǎn)、放縮、移動等功能進行不同精度、不同視角的可視化和精細化地質(zhì)解釋。

筆者基于當前相關(guān)工程領(lǐng)域?qū)υ擁椉夹g(shù)的剛性需求和3D科學(xué)可視化本身的優(yōu)勢，進而開展了橢圓周測線3D wiggle剖面繪制的算法研究工作，并編寫了相應(yīng)的matlab程序，最后通過與2D wiggle剖面對比，分析了本文研究成果的價值與優(yōu)勢。

1 橢圓樁底探地雷達反射法

探地雷達反射法是利用高頻電磁波(主頻10 MHz～1 000 MHz)以寬頻帶短脈沖形式，由地面通過發(fā)射天線將高頻電磁波傳入地下，經(jīng)地下地層或目標體反射回地面的電磁波被接收天線，當?shù)叵陆橘|(zhì)的波速已知時，可根據(jù)測得的脈沖波旅行時間，求出反射體的深度。電磁波在介質(zhì)中傳播時，其電磁波強度與wiggle將隨所通過介質(zhì)的電性及幾何形態(tài)而變化。所以根據(jù)接收到反射波的旅行時間、幅度及wiggle特征，便可推斷地下介質(zhì)的結(jié)構(gòu)特征。對于樁底下方基巖完整性情況的探測屬于精細化地質(zhì)探測，縱、橫向分辨率均要求較高，故關(guān)于樁底基巖探地雷達探測時，在保證探測深度情況下，盡量選擇頻率高的天線，通常選擇200 MHz～300 MHz范圍內(nèi)的天線便可滿足探測精度和深度的技術(shù)要求。在測線設(shè)計時，一般布置4條直測線呈“豐”字型、1條橢圓周測線，共5條測線(圖1)。當時間觸發(fā)連續(xù)采集的橢圓周測線布置時：橢圓周測線既以A點為起點又以其為終點，順時針旋轉(zhuǎn)布置測線，其半圓周測線半徑大小根據(jù)探地雷達天線體積大小合理設(shè)置；當按鍵觸發(fā)點采的橢圓周測線的測點布置時：布置方式與時間觸發(fā)連續(xù)采集的橢圓周測線布置方式一致，其離散測點還是布置在橢圓周測線上，其測點間距設(shè)置為10 cm～20 cm為宜。值得注意的是橢圓樁的橢圓并非數(shù)學(xué)上定義的橢圓，而是腰圓，其是指過圓心將一個圓平分成兩個半圓弧且相互反向平移，用二根等長平行線將兩個半圓弧的端點連接而形成的封閉圖形，其等長的平行線就是腰線(圖1)。

圖1 挖孔橢圓樁底探測雷達測線布置示意圖Fig.1 Arrangement diagram of GPR surveying lines for horizontal oval section digging pile bottom

為了保證時間觸發(fā)連續(xù)采集的距離歸一化精準定位和按鍵觸發(fā)點采的橢圓周測線的測點閉合，要求必須以天線上唯一的點作為參考點進行全過程打標，且天線上參考點的起始點和結(jié)束點一定要重疊閉合，其中按鍵觸發(fā)點采方式可不重疊閉合，但是第一個測點與最后一個測點的最小間距應(yīng)小于測點間距。

2 剖面計算

在利用探地雷達測線數(shù)據(jù)進行數(shù)據(jù)3D wiggle剖面可視化計算時，圖1所示的4條直測線，在三維直角坐標系下也是容易計算得到直測線各道的平衡點及其各采樣點坐標值，故筆者僅研究橢圓周測線的3D wiggle剖面生成方法。

橢圓樁底GPR 3D wiggle剖面的總體計算思路：先確定探地雷達橢圓周測線各道平衡點在x1o1y1直角坐標系中的坐標位置；然后計算各道的各采樣點在x1o1y1直角坐標系中的坐標位置；再將坐標變換到xoy直角坐標系中，最后通過時深轉(zhuǎn)換變換到三維直角坐標系o(x,y,z)中。

2.1 測道平衡點位置坐標計算

以左半圓的圓心為原點，測線1為x軸、測線4為y軸，建立新的x1o1y1直角坐標系(圖2)。不管是連續(xù)采集還是點采，一條已采集數(shù)據(jù)的橢圓樁的樁底周測線的道數(shù)是確定的。假設(shè)橢圓周測線道數(shù)為n，又由于橢圓樁的腰線BC或AD、左右半圓弧的半徑r均已知，于是每道所對應(yīng)的測線長度也是已知的。進一步假設(shè)：BC=AD=d，且橢圓周測線的起點和終點一般都設(shè)定在A點，故首先有橢圓周測線的第1道的平衡點落在A點上，為方便計算，假設(shè)正振幅方向為A點的右側(cè)、負振幅方向為A點的左側(cè)，其他測道與其保持一致。首先可得每一道對應(yīng)的橢圓周測線長度為2(πr+d)/n，于是基準振幅大小可定義為(πr+d)/n，然后計算確定各道平衡點在x1o1y1直角坐標系中坐標值。依據(jù)橢圓周測線幾何特征，該條測線先后有2種幾何線形，故可分4段分別求取其平衡點的相應(yīng)位置坐標值。

圖2 以樁底左半圓的圓心為原點的二維直角坐標系Fig.2 2D cartesian coordinate system with the center of the pile bottom as its left semicircle origin

2.1.1 落在左半圓弧上的平衡點位置坐標計算

首先，由于采集硬件等影響，探地雷達實際采集的道數(shù)據(jù)總是存在一定的漂移分量，即道平衡點兩側(cè)的最大正負振幅絕對值并不相等和正振幅代數(shù)和加上負振幅代數(shù)和并不等于零，所以在探地雷達數(shù)據(jù)正式處理之前，還必須消除直流漂移、歸一化、標準化預(yù)處理；此外定義當測道平衡點落在測線上時，則該測道落在該測線上。

其次，需要確定有多少道落在左半圓弧上，雖然針對某一具體的橢圓圓周測線已知有多少道、其第1道數(shù)據(jù)的平衡點設(shè)置剛好落在左半圓弧起點上，但是直接構(gòu)造一個通項數(shù)學(xué)公式來刻畫還是較困難。故利用前面假設(shè)的道數(shù)據(jù)的對稱性、落在測線上的定義與已有數(shù)學(xué)理論的內(nèi)在規(guī)律性。經(jīng)研究發(fā)現(xiàn)，橢圓周測線上多少道落在左半圓弧上的問題剛好滿足“四舍五入”數(shù)值計算的修約法則；另外，由于第1道落在左半圓弧上的A點上，依據(jù)上述定義，雖然算落在左半圓弧上，但是其只占用半道對應(yīng)的測線長度。故可構(gòu)造出計算落在左半圓弧上道數(shù)n1通項公式，見式(1)。

(1)

其中，round()表示“四舍五入”數(shù)值修約函數(shù)，進一步化簡為式(2)。

(2)

采用作者之前的類似幾何圖示分析法[13]，易得各道平衡點的在x1o1y1直角坐標系的位置坐標值公式，見式(3)～式(4)。

(3)

(4)

2.1.2 落在上腰線上的平衡點位置坐標計算

利用上述的“四舍五入”數(shù)值，計算修約法則與計算多少平衡點落在上腰線上的問題的等價性，基于n1已知，同理便可構(gòu)造出計算多少個平衡點落在上腰線上的通項公式(5)。

(5)

式(5)進一步化簡為式(6)。

(6)

計算上腰線上的各平衡點的位置，關(guān)鍵是找到其上的第1道的平衡點的位置，要直接確定上腰線上第1道的平衡點位置不可能的，只能從上一分段最后一道的平衡點位置坐標出發(fā)分析，從理論上分析可知，左半圓弧上最后道落在左半圓弧上的對應(yīng)最小弧長不超過(πr+d)/n，故其平衡點至B段的弧長計算公式為式(7)。

(7)

落在上腰線上的平衡點位置坐標計算公式，見式(8)～式(9)。

(8)

yn2=r

(9)

2.1.3 落在右半圓弧上的平衡點位置坐標計算

同理可得，構(gòu)造出計算多少個平衡點落在右半圓弧上的通項計算公式(10)。

(n2+n1)

(10)

(11)

落在右圓弧上的平衡點位置坐標計算公式為式(12)～式(13)。

(12)

(13)

2.1.4 落在下腰線上的平衡點位置坐標計算

基于已知的n=n4+n3+n2+n1，于是n4也是相當于已知的。

要計算下腰線上各道平衡點的坐標值，關(guān)鍵在于確定在第3分段上其最后一個道的平衡點的位置，可得第3分段上最后一個平衡點至D段圓弧計算公式為式(14)。

(14)

落在下腰線上的平衡點位置坐標計算公式為式(15)～式(16)。

(15)

yn4=-y

(16)

2.2 采樣點三維坐標計算

這里主要研究左右半圓弧上的各道采樣點的坐標值計算問題。

2.2.1 落在左半圓弧上的各道采樣點坐標計算

為了便于左半圓弧上各道的各采樣點坐標值計算公式推導(dǎo)，在圖2的基礎(chǔ)上繪制出相應(yīng)的幾何關(guān)系示意圖(圖3)。首先假設(shè)各振幅值走向為該道平衡點的切線F1WF2(W為切點，即平衡點)。另外，由本文相關(guān)定義可知，∠WO1F1=2(πr+d)/n，∠QO1F1=(πr+d)/n。

∴∠O1WH1=∠H1F1W

∴ΔO1H1W?ΔO1WF1

于是可得，∠F2WH2=∠F2WH1=∠WO1H1=2(πr+d)/2。為保持與圖3計算坐標系一致，并進一步假設(shè)正振幅在右側(cè)、負振幅在左側(cè)，于是當要計算振幅為正振幅(F1)、負振幅(F2)、的2種情況的坐標值，可利用三角函數(shù)關(guān)系可得具體計算公式為式(17)～式(20)。

圖3 左半圓弧測線段上采樣點二維坐標值計算原理示意圖Fig.3 Diagram of geometric relationship for sampling points on surveying line segment of the first quadrant

正振幅時：

xn1(sample,scan)=xn1+cos((n1(sample,scan)-1)*

(17)

yn1(sample,scan)=yn1-sin((n1(sample,scan)-1)*

(18)

負振幅時：

xn1(sample,scan)=xn1-cos((n1(sample,scan)-1)*

(19)

yn1(sample,scan)=yn1+sin((n1(sample,scan)-1)*

(20)

其中，xn1、yn1由式(3)～式(4)計算而得。

2.2.2 落在右半圓弧上的各道采樣點坐標計算

同理，為了便于右半圓弧上各道的各采樣點坐標值計算公式推導(dǎo)，在圖2的基礎(chǔ)上簡化繪制出相應(yīng)的幾何關(guān)系示意圖(圖4)。

圖4 右半圓弧測線段上采樣點二維坐標值計算原理示意圖Fig.4 Diagram of geometric relationship for sampling points on surveying line segment of the first quadrant

顯然，與圖3的幾何關(guān)系一致的，同樣存在ΔO2H1W?ΔO2WF1?ΔWH1F1?ΔWH2F2，且均為直角三角形。故關(guān)鍵在于計算出∠F1O2Q的值，由于已知道平衡點坐標值，通過反三角正切函數(shù)可計算得到∠F1O2Q。于是可得類似于式(17)～式(20)的計算右半圓弧上各道各采樣點計算公式。

至此，圓弧(2個半圓弧)和兩腰線(2條直線段)上所有各道各采樣點關(guān)于x1o1y1直角坐標系中的坐標值便可計算出來了。然后將x1+r和y1+r便可變換到xoy直角坐標系中。

最后,關(guān)于三維直角坐標系o(x,y,z)下的z坐標值的計算問題屬于探地雷達的時深轉(zhuǎn)換問題[13]。因為已知電磁波在真空中的傳播速度c(m·s-1)、時窗大小tw(s)、采樣點數(shù)m、相對介電常數(shù)εr，所以各道各采樣點對應(yīng)的深度可表示為式(21)。

(21)

該算法核心程序段為橢圓周測線在三維直角坐標系統(tǒng)下的x、y、z坐標值的科學(xué)計算及使用plot3()繪圖函數(shù)。

2.2.3 關(guān)健代碼示例

1)左半圓弧上各道各采樣點的振幅值對應(yīng)的x,y值的計算：

d1=zeros(samples,scans5zfd1);

h1=zeros(samples,scans5zfd1);

for ls=1:scans5zfd1 % 各分段起止測道數(shù)已確定

for lsc=1:samples

% 分為5種情況，并用line5表示橢圓周測線

if (ls-1)*ddfd==0

%第1種情況 振幅右正、左負

line5x(lsc,ls)=VVVV5(lsc,ls); line5y(lsc,ls)=-radiuszd;

elseif0<(ls-1)*ddfd

cot((n-1)*ddfd)*x-y=0

if VVVV5(lsc,ls)<0

d1(lsc,ls)=sin((ls-1)*ddfd)*abs(VVVV5(lsc,ls));

h1(lsc,ls)=cos((ls-1)*ddfd)*abs(VVVV5(lsc,ls));

line5x(lsc,ls)=s(ls,1)-h1(lsc,ls);

line5y(lsc,ls)=s(ls,2)+d1(lsc,ls);

elseif VVVV5(lsc,ls)>0

d1(lsc,ls)=sin((ls-1)*ddfd)*VVVV5(lsc,ls);

h1(lsc,ls)=cos((ls-1)*ddfd)*VVVV5(lsc,ls);

line5x(lsc,ls)=s(ls,1)+h1(lsc,ls);

line5y(lsc,ls)=s(ls,2)-d1(lsc,ls);

else

line5x=s(ls,1);

line5y=s(ls,2);

end

elseif (ls-1)*ddfd==pi/2

%第3種情況

line5x(lsc,ls)=-radiuszd;

line5y(lsc,ls)=-VVVV5(lsc,ls);

elseif pi/2<(ls-1)*ddfd

%第4種情況

if VVVV5(lsc,ls)<0

d1(lsc,ls)=sin((ls-1)*ddfd-pi/2)*abs(VVVV5(lsc,ls));

h1(lsc,ls)=cos((ls-1)*ddfd-pi/2)*abs(VVVV5(lsc,ls));

line5x(lsc,ls)=s(ls,1)+d1(lsc,ls);

line5y(lsc,ls)=s(ls,2)+h1(lsc,ls);

elseif VVVV5(lsc,ls)>0

d1(lsc,ls)=sin((ls-1)*ddfd-pi/2)*VVVV5(lsc,ls);

h1(lsc,ls)=cos((ls-1)*ddfd-pi/2)*VVVV5(lsc,ls);

line5x(lsc,ls)=s(ls,1)-d1(lsc,ls);

line5y(lsc,ls)=s(ls,2)-h1(lsc,ls);

else

line5x=s(ls,1);

line5y=s(ls,2);

end

else

% 第5種

line5x(lsc,ls)=-VVVV5(lsc,ls);

line5y(lsc,ls)=radiuszd;

end

end

end

2)下腰線上各道各采樣點的振幅值對應(yīng)的x,y值的計算：

for ls=scans5zfd3+scans5zfd2+scans5zfd1+…

1:scans5zfd4+scans5zfd3+scans5zfd2+scans5zfd1

for lsc=1:samples

if VVVV5(lsc,ls)>0

line5x(lsc,ls)=s(ls,1)+VVVV5(lsc,ls);

line5y(lsc,ls)=-radiuszd;

elseif VVVV5(lsc,ls)<0

line5x(lsc,ls)=s(ls,1)-abs(VVVV5(lsc,ls));

line5y(lsc,ls)=-radiuszd;

else

line5x(lsc,ls)=s(ls,1);

line5y(lsc,ls)=-radiuszd;

end

end

end

3)三維直角坐標系下Z坐標值的計算：

line5zz=ones(samples,scans5); % 賦初值

for jz5=1:scans5

for iz5=1:samples

line5zz(iz5,jz5)=(iz5*((timec*10^-9)/samples)*cc/(2*sqrt(cxdjd))); %由雙程走時進行時深轉(zhuǎn)換

end

end

4)3D和wiggle剖面繪制for循環(huán)段：

for j15=1:scans5 % 測線5按道數(shù)循環(huán)繪制曲線

plot3(line5xx(:,j15),line5yy(:,j15),line5zz(:,j15));

xlabel('X軸方向(AD)');

ylabel('Y軸方向(AB)');

zlabel('Z軸方向');

% set(findobj('Type','line'),'Color','k');

hold on;

axis square

end

set(gca,'YDir','reverse');

另外在編寫樁底橢圓周測線的3D wiggle剖面生成算法程序時，還運用了Matlab的Graphical user interface(GUI)功能，使該程序更有推廣使用價值。

3 試驗分析

以某鐵路橋梁橢圓樁的樁底探地雷達實測數(shù)據(jù)為例，研究其最終的對比試驗可視化顯示效果。

傳統(tǒng)的2D wiggle剖面顯示如圖5(a)所示，而由本文提出方法得到的3D wiggle剖面見圖5(b)所示。對比圖5(a)與圖5(b)分析可知，3D wiggle剖面可視化還原了探測時測線的各測道各采樣點實際空間位置，可以圖準確表述樁的哪一側(cè)面或哪一方位存在某異常，而傳統(tǒng)的2D wiggle剖面只能準確到距測線多少米處，具體位于樁的哪一面或哪個方位無法直接獲得，另外基于matlab圖形顯示窗口的圖形旋轉(zhuǎn)、縮放、移動等功能還能進行不同精度、不同視角的成圖顯示和精細化地質(zhì)解釋，甚至還可基于其“數(shù)-形”映射功能及生成的三維坐標值集及其對應(yīng)振幅值集，開展定量化地質(zhì)解釋工作。

圖5 實測數(shù)據(jù)的對比顯示Fig.5 Contrast display of measured data(a)基于radan5.0軟件生成的樁底橢圓周測線wiggle剖面;(b)本文提出的方法生成的樁底橢圓周測線3D wiggle剖面

在實際工作中，在數(shù)據(jù)處理完畢之后還可直接應(yīng)用根據(jù)本文提出的方法及之前作者的研究成果[13]進行橢圓周測線、圓周測線、直測線等組合或單獨的3D科學(xué)可視化，從而提高探地雷達樁底探測的技術(shù)水平。

4 結(jié)論

1)基于橢圓樁底的探地雷達周測線的幾何特征，構(gòu)造了其真實3D wiggle剖面的生成算法。

2)編寫了生成相應(yīng)的橢圓周測線wiggle剖面的matlab程序，同時為便于技術(shù)成果推廣應(yīng)用，并基于matlab的GUI功能開發(fā)了探地雷達檢測樁底基巖完整性3D可視化數(shù)據(jù)處理系統(tǒng)。

3)與傳統(tǒng)的橢圓周測線的2D wiggle剖面對比分析，①證實了本方法是正確、可行的；②說明了其在成果圖解譯時，可直接準確判別異常的空間位置或方位，大大地提高了地質(zhì)解釋環(huán)節(jié)的工作效率，同時基于matlab的圖形窗口的旋轉(zhuǎn)、放縮、移動等功能還可以進行精細化成果圖解譯；③給用戶較好的3D可視化視覺體驗。

4)基于標準歸一化的測道數(shù)據(jù)集，筆者巧妙應(yīng)用“四舍五入”數(shù)值修約法則與判別測道平衡點落在前一測線段還是后一測線段的問題，具有一一映射的等價關(guān)系，成功地解決了無論橢圓周測線長度或測道數(shù)怎么變化還是橢圓周測線怎么分段均能實現(xiàn)測道間距的精準等分的難題；然后利用橢圓周測線自身的幾何特征，基于“形數(shù)”的數(shù)形結(jié)合的思路成功解決了各道各采樣點空間坐標值的精確計算問題，故從應(yīng)用角度而言，具有一定的創(chuàng)新性。

5)隨著國內(nèi)、外大力建設(shè)高速鐵路等交通基礎(chǔ)設(shè)施，將有大量的挖孔橢圓樁需要開展樁底探地雷達探測工作，故本文提出的3D wiggle剖面生成方法具有較好的實用前景，值得大力推廣應(yīng)用。