GB-SAR影像坐標(biāo)到三維地形坐標(biāo)轉(zhuǎn)換方法

， (1.蘇州科技大學(xué) 環(huán)境科學(xué)與工程學(xué)院，蘇州 215009；2.武漢大學(xué) a.精密工程與工業(yè)測量國家測繪地理信息局重點實驗室; b.測繪學(xué)院，武漢 430079)

1 研究背景

地基合成孔徑雷達(dá)(Ground-based Synthetic Aperture Radar，GB-SAR)干涉測量技術(shù)是近十多年間發(fā)展起來的地面主動微波遙感探測技術(shù)。它主要應(yīng)用3項核心技術(shù)：①脈沖壓縮技術(shù)，如步進(jìn)頻率連續(xù)波技術(shù)(SFCW)、線調(diào)頻連續(xù)波技術(shù)(LFMCW)等，在獲取距離域高分辨率的同時又確保了較高的接收機(jī)功率和回波信號強(qiáng)度；②合成孔徑雷達(dá)技術(shù)，相比真實孔徑的地基雷達(dá)，合成孔徑能在較短時間內(nèi)采集雷達(dá)輻射場區(qū)的二維影像圖；③干涉測量技術(shù)，通過比較變形前后的觀測相位提取變形信息[1-3]。星載SAR受限于衛(wèi)星重訪周期和傳感器側(cè)視視角，在地表局部小區(qū)域的監(jiān)測應(yīng)用中不夠靈活。與星載SAR相比，GB-SAR結(jié)構(gòu)簡單、安裝便捷，它可以為監(jiān)測目標(biāo)建立特定的幾何場景，并根據(jù)目標(biāo)變形特征靈活確定時間基線。在GB-SAR連續(xù)影像序列的數(shù)據(jù)處理中，由于沒有空間基線，無需考慮軌道基線誤差、地平相位以及地形相位誤差影響。經(jīng)過十多年的研究發(fā)展和應(yīng)用實踐，GB-SAR業(yè)已在滑坡、支護(hù)邊坡、危巖體、水工結(jié)構(gòu)、露天煤礦和冰川移動監(jiān)測等方面展現(xiàn)了其特點和優(yōu)勢[4-9]，在星載InSAR和傳統(tǒng)大地測量監(jiān)測技術(shù)難以滿足時空分辨率和精度需求時能夠得到較好的應(yīng)用。

GB-SAR系統(tǒng)依據(jù)目標(biāo)到雷達(dá)中心的斜距和偏離雷達(dá)天線波束中心線的角度來分辨不同的監(jiān)測目標(biāo)。其成像投影方式不同于傳統(tǒng)地形圖的正射投影以及攝影測量的中心投影，在雷達(dá)影像的二維平面坐標(biāo)系中對變形目標(biāo)及其位置的識別一般采取人工判讀的方式進(jìn)行。這種方式要求分析人員具備一定的經(jīng)驗，否則容易造成變形目標(biāo)和區(qū)域的錯誤識別。因此，將雷達(dá)二維平面坐標(biāo)系統(tǒng)和常用的三維地形坐標(biāo)系聯(lián)系起來便非常有意義,不僅有利于提高變形體識別的正確性，也有助于長期趨勢性變形監(jiān)測、多種監(jiān)測手段的融合以及監(jiān)測成果直接或間接地轉(zhuǎn)化應(yīng)用。本文深入研究了GB-SAR影像坐標(biāo)到地方三維坐標(biāo)變換方法，成果有利于實際GB-SAR監(jiān)測應(yīng)用中對變形的三維解譯和分析，以更好地理解變形發(fā)展的時空規(guī)律。

2 相似變換參數(shù)估計坐標(biāo)變換及其誤差分析

星載SAR傳感器的目標(biāo)區(qū)域滿足天線遠(yuǎn)場近似條件，SAR影像坐標(biāo)系可視作規(guī)則格網(wǎng)。而GB-SAR作用距離相對較近且輻射視場距離跨度大，形成了特殊的扇形格網(wǎng)坐標(biāo)系。GB-SAR所探測的變形是監(jiān)測目標(biāo)實際變形在雷達(dá)視線向(Light of Sight，LOS)的投影，各像元以雷達(dá)中心點為中心，以到雷達(dá)中心的斜距為半徑沿圓弧線投影到成像投影面上。GB-SAR的投影方式不僅決定了距離向和方位向分辨單元的劃分形式，同時也決定了GB-SAR影像坐標(biāo)系統(tǒng)變換的方式方法。

為簡化分析，進(jìn)行如下假設(shè)：①小區(qū)域內(nèi)高斯-克呂格投影不存在投影變形；②雷達(dá)二維平面坐標(biāo)系和高斯平面坐標(biāo)系在同一個平面內(nèi)且軸系轉(zhuǎn)角為0，僅存在平移參數(shù)X0和Y0；③高程與高斯平面坐標(biāo)系構(gòu)成地面三維坐標(biāo)系，與XY軸構(gòu)成右手坐標(biāo)系。對于GB-SAR內(nèi)部的幾何關(guān)系，選定與高斯平面坐標(biāo)系平行的平面作為雷達(dá)成像投影面，則目標(biāo)點到雷達(dá)中心的相對高度就是目標(biāo)點到該平面的高度，三維地形坐標(biāo)系與雷達(dá)平面坐標(biāo)系的幾何關(guān)系如圖1所示。

圖1 三維地形坐標(biāo)系與GB-SAR影像平面坐標(biāo)系的幾何關(guān)系Fig.1 Geometric relationship between 3D terrain coordinate system and coordinate system of GB-SAR imaging plane

根據(jù)上述假設(shè)，可以計算得到雷達(dá)二維平面坐標(biāo)系中第i行、第j列的像元Tij在高斯平面坐標(biāo)系下的坐標(biāo)為

(1)

式中：θij為相對成像投影面的高度角；(Sij,βj)為目標(biāo)點在雷達(dá)二維平面坐標(biāo)(xij,yij)的極坐標(biāo)形式，即

利用相似變換參數(shù)估計的模型，由雷達(dá)二維平面坐標(biāo)計算高斯平面坐標(biāo)，即

(3)

經(jīng)過一定變形，計算縮放因子m。

(4)

代入式(1)，得

(5)

假設(shè)雷達(dá)中心線平行于高斯平面，雷達(dá)對同樣平行于高斯平面的一個目標(biāo)平面進(jìn)行觀測，雷達(dá)中心高出該目標(biāo)平面距離Hij，給定不同的相對高度Hij，計算縮放因子m。

從圖2中可以看出，隨著距離的增加，縮放因子m逐漸減小，減小的速率由快到慢，直至趨近于0。由上述分析可知，由相似變換參數(shù)估計方法算得的m值大小固定，并不適用于整個雷達(dá)輻射區(qū)域的坐標(biāo)變換。該方法僅適用于距離雷達(dá)較遠(yuǎn)、范圍較小的局部區(qū)域，并且相對于雷達(dá)中心的高度變化不能過于劇烈。

圖2 縮放因子m的變化規(guī)律Fig.2 Variation rule of scaling factor m

圖3 坐標(biāo)直接變換的誤差分析Fig.3 Error analysis for direct coordinate transformation

如果不考慮雷達(dá)影像投影方式，只進(jìn)行平面轉(zhuǎn)角旋轉(zhuǎn)和原點平移，將雷達(dá)二維平面坐標(biāo)系轉(zhuǎn)換至高斯平面坐標(biāo)系中，坐標(biāo)值產(chǎn)生的點位誤差實際上就是雷達(dá)沿圓弧線的投影位置與向高斯平面的正射投影之間的距離差值ΔSij，ΔSij與斜距和高度角之間的關(guān)系如圖3所示。易知，該方法在目標(biāo)區(qū)域分布距離雷達(dá)中心線較近時產(chǎn)生的點位誤差較小，例如高度角在5°以內(nèi)，直接坐標(biāo)變換方法產(chǎn)生的點位誤差在4 m以內(nèi)。但使目標(biāo)區(qū)域均在與雷達(dá)中心線高度角[-5°, 5°]之內(nèi)的要求是十分苛刻的，在實際中往往難以滿足。

3 外部高程信息輔助GB-SAR影像坐標(biāo)到地方三維坐標(biāo)轉(zhuǎn)換

由上述分析可知，GB-SAR成像特殊的投影方式導(dǎo)致影像坐標(biāo)相對實際平面坐標(biāo)發(fā)生偏移，且不同高度和位置的地形點成像后偏移的大小是不一樣的。實現(xiàn)GB-SAR影像坐標(biāo)到地方三維地形坐標(biāo)高精度的轉(zhuǎn)換計算，需要利用外部高程數(shù)據(jù)(如測區(qū)數(shù)字地面模型(DTM)/數(shù)字高程模型(DEM)、三維地形點云數(shù)據(jù)等)進(jìn)行。在三維地形坐標(biāo)系下的高程數(shù)據(jù)可依據(jù)GB-SAR投影原理計算為雷達(dá)投影面的平面坐標(biāo)，進(jìn)而利用其與GB-SAR影像像元坐標(biāo)的匹配關(guān)系實現(xiàn)逆向求解像元三維坐標(biāo)，主要包括以下幾個步驟。

3.1 軸系水平旋轉(zhuǎn)

在確定雷達(dá)中心的三維地形坐標(biāo)后，坐標(biāo)軸系的朝向不一致，因此首先需要進(jìn)行GB-SAR二維影像坐標(biāo)系到三維地形坐標(biāo)系的水平旋轉(zhuǎn)角度的計算。圖4顯示了GB-SAR影像平面坐標(biāo)系和地方平面坐標(biāo)系之間的相互關(guān)系。

圖4 坐標(biāo)系軸系關(guān)系與水平角度旋轉(zhuǎn)Fig.4 Horizontal angle rotation of coordinate axes

圖4中x-y平面是GB-SAR的雷達(dá)二維平面坐標(biāo)系；E-N平面是地方高斯平面坐標(biāo)系；θ是雷達(dá)二維平面坐標(biāo)系下目標(biāo)點到雷達(dá)中心線方向的偏角；αx-y是x-y坐標(biāo)系下目標(biāo)點的極角，αE-N是E-N坐標(biāo)系下目標(biāo)點的極角，而βE-N即為我們需要計算的軸系水平轉(zhuǎn)角。通常情況下，軸系水平轉(zhuǎn)角的計算利用公共點進(jìn)行[10]。但實際監(jiān)測工作中難有足夠小的自然地表可以作為控制點。一般需要布設(shè)角反射器并測量出該角反射器的三維空間坐標(biāo)，同時在雷達(dá)影像中識別出該角反射器相應(yīng)像元的雷達(dá)平面坐標(biāo)。

3.2 三維地形坐標(biāo)的極坐標(biāo)化

假設(shè)雷達(dá)影像中關(guān)注區(qū)域內(nèi)一像元點P，其雷達(dá)平面坐標(biāo)為(xP,yP)，相應(yīng)地表三維坐標(biāo)系坐標(biāo)為(EP,NP,ZP)。P點的雷達(dá)平面坐標(biāo)(xP,yP)也可以表示為極坐標(biāo)形式，如式(6)。

(6)

像元三維空間關(guān)系的恢復(fù)實際上就是求取雷達(dá)影像中的P點在三維地形點中的對應(yīng)點。為了達(dá)到目的，通常先將三維地形數(shù)據(jù)按照式(7)轉(zhuǎn)化為二維平面上的極坐標(biāo)形式，從而可以將三維地形數(shù)據(jù)的極坐標(biāo)數(shù)據(jù)與雷達(dá)像元的極坐標(biāo)直接進(jìn)行比較，這一過程稱為三維坐標(biāo)的極坐標(biāo)化。

(7)

式中：rTerrain，θTerrain分別為地形點極坐標(biāo)的極徑和極角；(E,N,Z)為地形點三維坐標(biāo)；(Eradar,Nradar,Zradar)為雷達(dá)中心在三維地形坐標(biāo)系下的坐標(biāo)。

3.3 GB-SAR成像投影面上像元與地形點的匹配

對三維地形數(shù)據(jù)極坐標(biāo)轉(zhuǎn)化后，在其極坐標(biāo)數(shù)據(jù)中按照最小距離準(zhǔn)則尋找P(rP,θP)的對應(yīng)點，并將該點的三維坐標(biāo)作為像元點P(xP,yP)對應(yīng)的三維空間坐標(biāo)系下的坐標(biāo)，如圖5。為保證結(jié)果的可靠性，在提取像元的三維信息時，亦可在目標(biāo)點云中選擇距離最小的前幾個點，計算其三維坐標(biāo)的平均值作為像元的三維坐標(biāo)。

圖5 極坐標(biāo)化點云與雷達(dá)像元的匹配Fig.5 Coordinate polarization of point cloud and matching with GB-SAR pixel

4 試驗案例分析

隔河巖水利樞紐位于湖北省長陽縣境內(nèi)長江支流清江干流上，上距恩施市207 km，下距高壩洲水電站50 km。工程開發(fā)任務(wù)以發(fā)電為主，兼有防洪、航運等綜合效益。樞紐工程由混凝土重力拱壩、泄水建筑物、右岸引水式水電站和左岸垂直升船機(jī)組成。水電站位于右岸，引水式地面廠房，4條直徑為9.5 m的隧洞接直徑8 m的壓力鋼管，分別接至4臺300 MW的混流式水輪機(jī)組。隧洞采用預(yù)應(yīng)力混凝土襯砌。廠房和壓力鋼管開挖形成170 m的高邊坡，采用混凝土局部置換、設(shè)置預(yù)應(yīng)力錨束及加強(qiáng)山體排水等措施處理。

本試驗主要對壩體與壩區(qū)右岸發(fā)電廠房南側(cè)高邊坡進(jìn)行監(jiān)測。雷達(dá)設(shè)備安裝于大壩下游左岸、距壩體約1 km處，導(dǎo)軌安裝基本處于水平狀態(tài)，雷達(dá)視線中心朝向壩體右側(cè)，壩體與右岸高邊坡均處于雷達(dá)最優(yōu)視場范圍內(nèi)，如圖6。在進(jìn)行連續(xù)變形監(jiān)測試驗之前，首先利用高精度全站儀獲取了導(dǎo)軌空間坐標(biāo)。在導(dǎo)軌兩端固定刻畫位置處安置固定高度的棱鏡基座，分別測量出兩端棱鏡中心的三維坐標(biāo)。經(jīng)計算雷達(dá)中心三維坐標(biāo)為(2 200.608 1，1 109.037 4，106.044 6)(單位：m)。

圖6 壩區(qū)平面坐標(biāo)系與雷達(dá)平面坐標(biāo)系Fig.6 Plane coordinate system of dam area and GB-SAR plane coordinate system

表1 測站坐標(biāo)系與地方坐標(biāo)系控制點對應(yīng)關(guān)系Table 1 Coordinates of control points in scanner coordinate system and dam area coordinate system

4.1 三維激光掃描儀邊坡高程提取

試驗中利用地面激光三維掃描儀RIEGL VZ400采集了高邊坡及其鄰近地物的點云數(shù)據(jù)，如圖7。

圖7 三維激光掃描儀地形點云數(shù)據(jù)采集Fig.7 Acquisition of terrain point cloud data using 3D laser scanner

直接獲取的三維點云數(shù)據(jù)仍然處于掃描儀三維坐標(biāo)系下，為將地面三維激光掃描儀的點云坐標(biāo)轉(zhuǎn)換到地方三維坐標(biāo)系，在邊坡上的幾個控制點上布設(shè)反射片作為控制點，如圖7(b)，本次試驗中共獲取了4個較為理想的控制點，其掃描儀三維坐標(biāo)和壩區(qū)三維坐標(biāo)列于表1。

點云數(shù)據(jù)預(yù)處理時，將邊坡表面的植被剔除，僅保留邊坡坡面三維坐標(biāo)數(shù)據(jù)。在得到控制點2套坐標(biāo)之后，通過應(yīng)用大轉(zhuǎn)角的七參數(shù)坐標(biāo)變換方法，得到邊坡點云在地方壩區(qū)三維坐標(biāo)系下的坐標(biāo)，如圖8所示。

圖8 壩區(qū)坐標(biāo)系下邊坡點云的三維坐標(biāo)Fig.8 Three-dimensional coordinates of point cloud of a slope in dam area

4.2 軸系水平轉(zhuǎn)角的計算與三維地形坐標(biāo)的極坐標(biāo)化

在完成空間坐標(biāo)系的變換之后，三維地形坐標(biāo)歸算至以雷達(dá)中心為原點的水平面上，如圖9(a)。由于壩區(qū)坐標(biāo)系的橫軸和縱軸定義剛好與雷達(dá)平面坐標(biāo)系的橫軸和縱軸反向，所以需要對平移計算后的x,y坐標(biāo)值做反號處理，便得到局部坐標(biāo)系下的地形點云，如圖9(b)。

圖9 壩區(qū)坐標(biāo)系、局部坐標(biāo)系和雷達(dá)坐標(biāo)系下的地形點云坐標(biāo)Fig.9 Coordinates of terrain point cloud in dam area coordinate system, local coordinate system,and GB-SAR coordinate system

壩區(qū)坐標(biāo)系和局部坐標(biāo)系中，電廠廠房結(jié)構(gòu)清晰可辨。此時局部坐標(biāo)系與雷達(dá)平面坐標(biāo)系的縱軸之間只相差一個數(shù)值為銳角的水平轉(zhuǎn)角。三維坐標(biāo)的極坐標(biāo)化后的位置只與目標(biāo)點到原點的斜距和相對于過原點的水平面的高度有關(guān)，而與y軸的朝向沒有關(guān)系。點云數(shù)據(jù)中特征線由大量三維坐標(biāo)形成，因此可以較為精確地確定直線極坐標(biāo)化后的方位角，而建筑結(jié)構(gòu)在雷達(dá)影像中所成的像同樣占據(jù)多個分辨單元，也能夠較為精確地確定其方位角。因此本文利用場區(qū)特征直線計算軸系水平轉(zhuǎn)角，為0.177 142 rad。利用坐標(biāo)旋轉(zhuǎn)公式計算極坐標(biāo)化后的坐標(biāo)順時針旋轉(zhuǎn)0.177 142 rad后的坐標(biāo)，如圖9(c)。

4.3 GB-SAR影像像元高程提取

在將三維坐標(biāo)極坐標(biāo)化并完成軸系水平角旋轉(zhuǎn)之后，便可以從得到的雷達(dá)平面坐標(biāo)系中尋找GB-SAR影像各像元的對應(yīng)點。如前文所述，試驗中采用到待求像元距離最小作為對應(yīng)三維地形點目標(biāo)選擇的準(zhǔn)則，直接在極坐標(biāo)化后的坐標(biāo)系中尋找與雷達(dá)影像像元坐標(biāo)距離最近的目標(biāo)點，將該點相應(yīng)的極坐標(biāo)化前的三維坐標(biāo)作為像元的新坐標(biāo)。至此便完成了GB-SAR影像像元的高程提取工作，只要導(dǎo)軌未發(fā)生變動，監(jiān)測獲取的所有影像(例如信號強(qiáng)度圖、相干圖和變形分布圖等)都可以按照上述工作求得的關(guān)系映射到局部三維表面模型上。

選取高邊坡區(qū)域的雷達(dá)影像，熱信噪比(Thermal Signal-to-Noise Ratio，TSNR)閾值設(shè)為10 dB，剔除影像中部分信號較弱的區(qū)域。對經(jīng)閾值處理后的TSNR影像坐標(biāo)值計算相應(yīng)的壩區(qū)三維坐標(biāo)XYH并繪制到高邊坡的三維表面模型上。為去除不在模型區(qū)域的像元，在計算前設(shè)定淘汰距離為3 m，超過該距離的像元不會出現(xiàn)在三維模型上，如圖10。由于試驗中高邊坡整體偏離雷達(dá)中心線的角度有些大，只有高邊坡右側(cè)部分區(qū)域的雷達(dá)反射信號較為理想。第三級馬道邊緣較為明顯，并直接連接4個引水隧道外壁，具有相對較高的TSNR值。該區(qū)域GB-SAR影像像元與實際地表結(jié)構(gòu)符合較好。

圖10 GB-SAR TSNR圖到三維表面模型的映射Fig.10 Mapping of GB-SAR TSNR graph to 3D surface model

5 結(jié) 論

(1)由于GB-SAR成像過程投影方式的特殊性，各像元在影像坐標(biāo)系中的位置相比于實際地形坐標(biāo)發(fā)生徑向偏移，使得僅靠人工識圖方式對變形區(qū)域進(jìn)行識別與劃分容易產(chǎn)生偏差或錯誤。為更加直接和有效地解譯變形監(jiān)測成果，有必要將雷達(dá)影像坐標(biāo)系統(tǒng)與常用坐標(biāo)系統(tǒng)統(tǒng)一起來。

(2)利用相似變換參數(shù)估計方法進(jìn)行GB-SAR影像坐標(biāo)變換模型不嚴(yán)密，在雷達(dá)監(jiān)測視場有一定跨度時易產(chǎn)生較大誤差。精確的GB-SAR影像坐標(biāo)到地方三維坐標(biāo)變換必須顧及各像元相對于雷達(dá)中心的高度信息。本文研究的顧及外部高程信息的GB-SAR影像坐標(biāo)三維變換方法在理論上能夠降低坐標(biāo)轉(zhuǎn)換誤差，但目前缺少對轉(zhuǎn)換精度的直接定量分析。后期筆者考慮引入角反射器陣列以定量分析和評價GB-SAR影像坐標(biāo)系統(tǒng)的轉(zhuǎn)換精度。

參考文獻(xiàn)：

[1] PIERACCINI M, FRATINI M, PARRINI F,etal. High-speed CW Step-frequency Coherent Radar for Dynamic Monitoring of Civil Engineering Structures[J]. Electronics Letters，2004,40(14):907-908.

[2] 王 鵬,周 校.地基SAR干涉測量原理及其形變監(jiān)測應(yīng)用研究[J].測繪信息與工程,2012,37(4):22-25.

[3] R?DELSPERGER S. Real-time Processing of Ground Based Synthetic Aperture Radar (GB-SAR) Measurements [D]. Darmstadt, Germany: Technische Universit?t Darmstadt, 2011.

[4] MARGOTTINI C, ANTIDZE N, COROMINAS J,etal. Landslide Hazard, Monitoring and Conservation Strategy for the Safeguard of Vardzia Byzantine Monastery Complex, Georgia[J]. Landslides, 2015, 12(1): 193-204.

[5] ATZENI C, BARLA M, PIERACCINI M,etal. Early Warning Monitoring of Natural and Engineered Slopes with Ground-based Synthetic-aperture Radar[J]. Rock Mechanics and Rock Engineering, 2015, 48(1): 235-246.

[6] 邱志偉,岳建平,汪學(xué)琴.地基雷達(dá)系統(tǒng)IBIS-L在大壩變形監(jiān)測中的應(yīng)用[J].長江科學(xué)院院報,2014,31(10):104-107.

[7] TARCHI D, ANTONELLO G, CASAGLI N,etal. On the Use of Ground-based SAR Interferometry for Slope Failure Early Warning: the Cortenova Rock Slide (Italy) [M]∥Landslides. Berlin：Springer Berlin Heidelberg，2005：337-342.

[8] TRAGLIA F D, NOLESINI T, INTRIERI E,etal. Review of Ten Years of Volcano Deformations Recorded by the Ground-based InSAR Monitoring System at Stromboli Volcano: A Tool to Mitigate Volcano Flank Dynamics and Intense Volcanic Activity[J]. Earth-Science Reviews, 2014,139:317-335.

[9] SEVERIN J, EBERHARDT E, LEONI L,etal. Development and Application of a Pseudo-3D Pit Slope Displacement Map Derived from Ground-based Radar[J]. Engineering Geology, 2014, 181: 202-211.