任麗麗,李家存,張 迪,鐘若飛,曾凡洋
(首都師范大學資源環境與旅游學院三維信息獲取與應用教育部重點實驗室,北京100048)
探地雷達作為一種新型的無損探測工具,由于操作簡單及較好的分辨率,被廣泛應用于近地表下目標物或介質層的探測。探地雷達獲取的圖像為水平方向和豎直方向成比例的二維灰度圖,當地形比較平坦時,探地雷達圖像可以將真實的地下構造反映出來。但在實際野外探測中,經常會遇到地形起伏變化的情況。由于地形變化影響及圖像顯示軟件的限制,獲取的雷達圖像無法將地下界面中的真實形態反映出來,從而影響到圖像的解譯和目標物的精確定位。因此,需對探地雷達圖像做地形校正處理[1]。
地形校正的過程,即探地雷達圖像與地形數據匹配的過程。對于地形數據的獲取,目前多采用激光水平儀[2]、全站儀[3]、GPS 和 DGPS[4]、傾角儀和里程計[5]等傳統的測量方法實現。這些傳統測量方法存在著數據采集工作量大,地形數據和探地雷達圖像匹配易受人為因素干擾等缺點。本文結合新型的激光測繪技術,提出一種采用激光點云實現探地雷達圖像地形校正的方法。與傳統測量方法相比,激光掃描技術可以快速、高效地獲取地形數據,通過測量線上均勻分布的離散點可以快速實現探地雷達圖像與地形數據的精確匹配,從而實現對探地雷達圖像的地形校正。
探地雷達地形校正的方法借鑒地震波靜校正的思想,基本原理是將探地雷達各道數據的雙程傳播時間統一校正到距空氣-大地交界面上方一定距離處的一個水平基準面上,這個基準參考面通常為測量線的最高或者最低點所在的水平面[6]。圖1為探地雷達在地形起伏情況下采集數據,假設地面下某深度處存在均勻水平分布的介質層。探地雷達獲取的二維時間剖面如圖2(a)所示,由于地形不斷起伏變化,電磁波在介質中的傳播時間發生變化,從而導致探地雷達圖像上水平分布的介質層發生扭曲。根據時間移位原理,選擇地形最高點所在的水平面為參考面,將探地雷達在高程上的傳播時間都校正到這一參考平面上。運用公式(1)結合線性插值的方法可計算出探地雷達每道數據相對于基準參考面的時間差,然后對探地雷達圖像進行校正,校正后的時間剖面如圖2(b)所示,探地雷達圖像與真實的地形變化相一致。


圖1 探地雷達在起伏地形下采集數據

圖2 地形校正前后的二維時間剖面
測量區的地形數據是由車載激光掃描系統獲取的,該系統主要由360°激光掃描儀、組合導航系統(IMU/GPS)和全景相機等傳感器組成[7]。激光掃描儀獲取空間物體表面的距離、角度和強度信息;組合導航系統記錄載體實時的位置和姿態信息;全景相機獲取空間物體的彩色紋理信息。激光點云數據本身是沒有彩色紋理信息的,這樣就很難提取地面上對應的標識點坐標信息。因此要對激光點云數據進行賦相機彩色紋理處理,以便實現標識點坐標信息的提取。
根據曝光時刻全景相機中心、全景球面上像點、物點三點共線的原理,通過一系列的坐標轉換,實現全景影像與對應的點云匹配,使激光點云賦上彩色紋理[8]。其步驟如下:首先將在大地坐標系下的激光數據轉換到以全景球心為中心坐標系下:

其中,(XG,YG,ZG)為大地坐標系下的坐標;(dx,dy,dz)為當前全景球球心的大地坐標;(a1,b1,c1,a2,b2,c2,a3,b3,c3)為旋轉矩陣的系數,由全景影像的三個姿態角φ(橫滾角)、ω(俯仰角)、κ(航向角)確定。
其次,由坐標(XC,YC,ZC)可以計算對應的像點在全景影像上的像素坐標(m,n)


其中,Width為全景影像的寬;Height為全景影像的高。將對應像點的顏色值賦給該點,實現激光點云賦彩色:

其中,RGB(XS,YS,ZS)表示(XS,YS,ZS)的 RGB 顏色值;N為全景影像編號;RGB(m,n,N)表示全景影像上像素(m,n)的RGB顏色值。
測量區域的地形數據和探地雷達獲取圖像之間的同步,是通過測量線上均勻分布的反光標識點實現的。數據采集之前,沿探地雷達測線方向每隔20 cm的距離布置標識點,共46個,如圖3(a)所示。測量線上標志點分布越多,探地雷達圖像與地形數據匹配的精度越高。采集數據時,當天線的中心與測量線的標識點重合時,通過探地雷達采集軟件記錄此時雷達圖像上的水平位置。然后從彩色點云數據上提取出所有標識點的精確高程信息,通過均勻分布的標識點最終實現探地雷達圖像與地形數據精確匹配,如圖3(b)所示。選擇測量線上的最高點所在的平面作為參考平面,根據已知目標換算的方法求出電磁波的傳播速度,那么離散點相對于參考面在順直方向上的電磁波傳播雙程時間差可以計算出來,如表1所示。最后運用分段插值方法(如公式(6)所示)計算出探地雷達每一道數據相對于參考平面的時間差,從而實現探地雷達圖像在豎直方向上的校正處理:


表1 離散點相對于基準參考面的時間差

圖3 實驗區的彩色點云數據
圖4(a)為經過濾波處理后的雷達剖面,為水平方向和豎直方向成比例的矩形剖面圖。水平距離3~4 m和7~8 m區域內存在明顯的雙曲線反射,8~10 m區域為土質疏松區,含水量比較大,反射波波形較混亂,呈高頻狀。根據電磁波的振幅和相位變化特征,用紅色的虛線繪制出地面下介質層的分界線,隨著水平距離的增加深度不斷加深。圖4(b)為經過地形校正后的探地雷達剖面。結合實驗區的情況,可以確定出沿測線方向3~4 m處雙曲線反射為部分裸露在地表上的樹根,7~8 m處雙曲線反射為埋在地面下的金屬材質自來水管。跟地形校正前的剖面圖相比,地形校正后的探地雷達圖像將實驗區實地形變化情況真實反映出來。通過圖像上對地形變化的補償校正,使地下目標物到地面的距離得到修正,定位的精度進一步提高;同時使地面下介質層分界線發生變化,與地面下真實的介質分布相一致。

圖4 地形校正前后的雷達剖面
本文利用激光測量技術具有快速、高效獲取數據及數據分辨率高的優點,通過獲取實驗區的彩色點云數據實現了探地雷達圖像的地形校正。實驗結果證明,地形校正效果明顯。但此方法僅適用于地形變化傾斜度較小情況,當地形變化傾斜度較大時,需要在地形校正基礎上做天線傾斜校正處理。
[1] ZHOU Hui,WANG Zhao - lei,HAN Bo,et al.Terrain correction and migration of GPR profile fulfilled simultaneously using reverse - time migration[J].Journal of Jilin University:Earth Science Edition,2004,34(3):459 -463.(in Chinese)周輝,王兆磊,韓波,等.同時實現地質雷達數據地形校正和偏移成像方法[J].吉林大學學報:地球科學版,2004,34(3):459 -463.
[2] D Percy,C Peterson.Rapid acquisition of ground penetrating radar enabled by LIDAR[J].Digital Mapping Techniques,2006,8(10):183 -185.
[3] Mercedes Solla,Henrique Lorenzo,Alexandre Novo,et al.Evaluation of ancient structures by GPR(ground penetrating radar):The arch bridges of Galicia(Spain)[J].Scientific Research and Essays,2011,6(8):1877 - 1884.
[4] S Urbini,J A Baskaradas.GPR as an effective tool for safety and glacier characterization:experiences and future development[J].Proc of the XIII Int.Conf.on Ground Penetrating Radar,2010:489 -494.
[5] Vitalii P,Volodymyr I,Sergiy K,et al.Topographic correction of GPR profile based on odometer and inclinometer data[C].14th International Conference on Ground Penetrating Radar,2012.
[6] Yilmaz O.Seismicdata processing[M].Tulsa:Society of Exploration Geophysicists,1987.
[7] ZHANG Di,ZHONG Ruo - Fei,LU Xu - Wei.Time synchronization method of laser scanner without external trigger[J].Laser & Infrared,2013,43(6):618 - 621.(in Chinese)張迪,鐘若飛,魯旭偉.無外觸發的激光掃描儀的時間同步方法研究[J].激光與紅外,2013,43(6):618-621.
[8] ZHANG Yu - xiang,ZHANG Xing - jun.Study on reconstruction property of free surface fitting for point clouds with laser scanning[J].Laser & Infrared,2011,41(3):351 -356.(in Chinese)張玉香,張興軍.采用激光掃描點云擬合自由曲面的重構特征研究[J].激光與紅外,2011,41(3):351-356.