王亞偉 鐘選明 周亮 廖成

(西南交通大學(xué)電磁場(chǎng)與微波技術(shù)研究所,成都 610031)

引 言

隨著信息化社會(huì)的發(fā)展,人們對(duì)通信質(zhì)量的要求越來(lái)越高,無(wú)線(xiàn)電波作為信息傳遞的一種重要載體,其傳播特性越來(lái)越受到人們的重視,準(zhǔn)確地預(yù)測(cè)電波傳播特性能夠?yàn)榛镜倪x址、網(wǎng)絡(luò)的優(yōu)化、電臺(tái)的配置等提供有效的理論依據(jù)和工程指導(dǎo)[1]. 從20世紀(jì)70年代開(kāi)始,國(guó)內(nèi)外相關(guān)學(xué)者對(duì)電波傳播特性預(yù)測(cè)展開(kāi)了大量的研究工作,建立了很多的電波預(yù)測(cè)模型,如射線(xiàn)跟蹤法、時(shí)域有限差分法等,但這些模型在處理復(fù)雜地理環(huán)境中的電波傳播問(wèn)題時(shí)難以在精度和效率之間達(dá)到一個(gè)很好的平衡. 由波動(dòng)方程推導(dǎo)而來(lái)的拋物方程方法(Parabolic Equation Method, PEM)能夠快速準(zhǔn)確地處理復(fù)雜地表邊界和復(fù)雜大氣結(jié)構(gòu)對(duì)電波傳播的影響[2-3],在求解大區(qū)域復(fù)雜環(huán)境中的電波傳播問(wèn)題時(shí)得到了廣泛的應(yīng)用.

不規(guī)則的地形,多樣化的地表覆蓋物等是影響電波傳播特性的主要因素. 目前,已有許多國(guó)內(nèi)外學(xué)者基于拋物方程對(duì)地理環(huán)境中的電波傳播問(wèn)題展開(kāi)了研究,Holm[4]、郭建炎[5]等人將森林視為有耗的介質(zhì)層,采用二維拋物方程研究了不規(guī)則地形中森林覆蓋環(huán)境下的電波傳播問(wèn)題,并與實(shí)測(cè)結(jié)果對(duì)比,驗(yàn)證了模型的正確性;張青洪[6]等人將森林看成空氣和水的混合物,建立了森林等效介電常數(shù)模型,對(duì)拋物方程的森林模型進(jìn)行了改進(jìn),這些研究一般考慮的地表環(huán)境比較單一,電磁建模也比較粗糙. 白瑞杰[7]等人基于數(shù)字高程模型對(duì)傳播區(qū)域的電波特性做了預(yù)測(cè)分析,由于數(shù)字高程模型數(shù)據(jù)描述的是地形起伏的狀況,缺少森林、植被、建筑等地表覆蓋物的信息,難以滿(mǎn)足城市小區(qū)域電波傳播預(yù)測(cè)的電磁建模要求. 因此,研究更加接近真實(shí)場(chǎng)景下的電磁建模很有必要.

LiDAR點(diǎn)云是由機(jī)載激光雷達(dá)系統(tǒng)通過(guò)發(fā)射和接收激光脈沖獲取的地表密集的高精度的三維點(diǎn)坐標(biāo),是對(duì)地形、地物特征的一種高精度的描述[8]. 通過(guò)對(duì)其分析處理,不僅可以獲得地形起伏的信息,還可以獲得植被、建筑等地物信息,將其應(yīng)用于復(fù)雜地理環(huán)境的電磁建模,將會(huì)使電波傳播的預(yù)測(cè)更加準(zhǔn)確,因此本文研究了基于LiDAR點(diǎn)云的地理電磁環(huán)境建模. 采用數(shù)學(xué)形態(tài)學(xué)方法對(duì)格網(wǎng)化后的LiDAR點(diǎn)云做了濾波等處理,獲取地形、地物等信息,結(jié)合相應(yīng)物質(zhì)的電磁參數(shù),進(jìn)行了地理環(huán)境的電磁建模,在此基礎(chǔ)上,利用三維拋物方程模型對(duì)地理環(huán)境中的電波傳播特性進(jìn)行了仿真分析.

1 三維拋物方程模型

在笛卡爾直角坐標(biāo)系中,假定電磁場(chǎng)的時(shí)諧因子為e-iωt,電場(chǎng)或磁場(chǎng)分量ψ滿(mǎn)足如下形式的標(biāo)量Helmholtz方程:

(1)

式中:n為媒質(zhì)的折射指數(shù);k0為自由空間中的波數(shù).定義波函數(shù)u(x,y,z)=ψe-ik0x(x,y,z),將其代入式(1),并進(jìn)行因式分解,可以得到

(2)

式(2)中,Q為偽微分算子,

(3)

僅考慮前向傳播,對(duì)Q作Feit-Fleck型近似,可得到如下形式的三維拋物方程:

(4)

式(4)的數(shù)值解法主要為有限差分 (Finite Difference, FD[9]) 法和分步傅里葉變換 (Split-step Fourier Transform, SSFT[10]) 法,它們都是步進(jìn)迭代算法. 其中FD算法的計(jì)算步長(zhǎng)受波長(zhǎng)嚴(yán)格限制,主要用于計(jì)算目標(biāo)的雷達(dá)散射截面(Radar Cross Section, RCS)以及城鎮(zhèn)小尺度電波傳播特性的預(yù)測(cè). SSFT算法引入了快速傅里葉變換(Fast Fourier Transform, FFT)技術(shù),在傳播方向上可以取較大的步長(zhǎng),能夠快速地預(yù)測(cè)復(fù)雜地理環(huán)境中的電波傳播特性,其解的形式如下[10]:

u(x+Δx,y,z)=eik0Δx(n-2)F-1{eikxΔxF[u(x,y,z)]}.

(5)

式中:kx表示k0在x方向上的分量;F和F-1分別表示二維傅里葉正逆變換.知道初始場(chǎng),利用式(5)結(jié)合邊界條件可迭代遞推求得整個(gè)空間中的場(chǎng)強(qiáng)分布.

電磁波在空間遠(yuǎn)距離傳輸時(shí),森林的影響是一個(gè)不容忽視的因素. 目前在研究森林對(duì)電波傳播的影響時(shí),通常把森林等效為均勻有耗的介質(zhì)層[4-5],其等效的折射系數(shù)可表示成

(6)

式中:εr,σ分別表示森林的等效相對(duì)介電常數(shù)和電導(dǎo)率. 在迭代式(5)的過(guò)程中,通過(guò)式(6)不斷地修正森林中的折射系數(shù),則可處理森林中的電波傳播問(wèn)題. 通常森林中的等效折射系數(shù)趨近1,滿(mǎn)足拋物方程的近似條件.

2 地理電磁環(huán)境建模

機(jī)載LiDAR是一種新型的主動(dòng)式遙感技術(shù),它通過(guò)測(cè)量激光脈沖的傳播時(shí)間,結(jié)合POS系統(tǒng)提供的定位定姿數(shù)據(jù),能夠直接、快速獲取地面的三維坐標(biāo)[11]. LiDAR獲取的點(diǎn)云文件記錄了一系列被掃描物體表面的三維離散點(diǎn)坐標(biāo),其不規(guī)則地分布在三維空間中. 不同的地形、地物對(duì)應(yīng)不同的激光腳點(diǎn)分布,反應(yīng)到點(diǎn)云數(shù)據(jù)就會(huì)呈現(xiàn)出不同的幾何拓?fù)湫螒B(tài). 通常情況下,地面點(diǎn)云的高程變化比較平緩,點(diǎn)的分布比較均勻;樹(shù)木點(diǎn)云分布比較散亂,一般在邊緣處會(huì)有高程突變;建筑物點(diǎn)云整體上幾何形態(tài)較為規(guī)整,高程變化也較為一致. 水體對(duì)激光信號(hào)(波長(zhǎng)為1 064 nm)的吸收,則會(huì)造成大面積的水體空洞區(qū). 依據(jù)這些地物呈現(xiàn)出的不同特征,采用一定的技術(shù)手段可以對(duì)LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行分類(lèi).

2.1 點(diǎn)云濾波

點(diǎn)云分類(lèi)中重要的一步是從點(diǎn)云數(shù)據(jù)分離出地面點(diǎn)和地物點(diǎn),也就是點(diǎn)云濾波. 目前國(guó)內(nèi)外學(xué)者提出許多種的濾波方法,其中具有代表性的算法有移動(dòng)曲面擬合法、迭代三角網(wǎng)加密法、數(shù)學(xué)形態(tài)學(xué)方法等[10]. 為了適應(yīng)拋物方程計(jì)算傳播迭代步進(jìn)的需要,本文采用了一種基于網(wǎng)格操作的數(shù)學(xué)形態(tài)學(xué)濾波算法[13].

數(shù)學(xué)形態(tài)學(xué)是一門(mén)建立在集合論的學(xué)科,目前已被廣泛應(yīng)用于各類(lèi)圖像分析和處理中,其基本操作有四種,分別為腐蝕、膨脹、開(kāi)和閉運(yùn)算. 對(duì)于LiDAR觀測(cè)點(diǎn)p(x,y,z),膨脹(Dilation)運(yùn)算和腐蝕(Erosion)運(yùn)算分別定義為:

(7)

(8)

點(diǎn)(xp,yp,zp)表示在窗口(也稱(chēng)結(jié)構(gòu)元)大小為w時(shí)點(diǎn)p的鄰域點(diǎn),理論上,窗口形狀可以是任意的,實(shí)際處理中,我們通常選擇簡(jiǎn)單的形狀,如矩形、圓形、菱形等. 由式(7)、(8)可以推斷出,膨脹運(yùn)算和腐蝕運(yùn)算分別能夠獲得鄰域窗口內(nèi)高程最大值和最小值. 基于膨脹和腐蝕運(yùn)算,可以得到用于處理LiDAR數(shù)據(jù)的開(kāi)運(yùn)算和閉運(yùn)算,開(kāi)運(yùn)算是先對(duì)給定的數(shù)據(jù)執(zhí)行腐蝕運(yùn)算,然后執(zhí)行膨脹運(yùn)算,閉運(yùn)算則剛好相反.

數(shù)學(xué)形態(tài)學(xué)濾波的主要思想是將地物(樹(shù)木、建筑物等)通過(guò)腐蝕運(yùn)算腐蝕至地表,然后通過(guò)膨脹運(yùn)算進(jìn)行恢復(fù)[14]. 圖1是采用線(xiàn)形結(jié)構(gòu)元對(duì)點(diǎn)云數(shù)據(jù)的某一縱剖面執(zhí)行開(kāi)運(yùn)算濾波的示意圖. 當(dāng)有多個(gè)點(diǎn)落入窗口中時(shí),執(zhí)行開(kāi)運(yùn)算后,將保留高程值最小的數(shù)據(jù)點(diǎn). 由于樹(shù)木的尺寸小于結(jié)構(gòu)元尺寸,經(jīng)過(guò)腐蝕后首先被去除,對(duì)于尺寸大于結(jié)構(gòu)元的建筑物,則會(huì)在膨脹運(yùn)算中得到重建.如果某點(diǎn)濾波前后的高程差在閾值范圍內(nèi),則判定該點(diǎn)為地面點(diǎn),否則判定為地物點(diǎn),并標(biāo)記. 可以看出,開(kāi)運(yùn)算能夠很好分離出地物點(diǎn)以及地面點(diǎn),并且能夠獲得較為平滑的地形.

圖1 LiDAR點(diǎn)云開(kāi)運(yùn)算濾波示意圖

結(jié)構(gòu)元尺寸的選擇決定了最終的濾波效果,當(dāng)采用較小的濾波窗口時(shí),尺寸大于窗口的地物(如建筑物)得不到很好的分離,采用較大的濾波窗口在分離地物的同時(shí),也會(huì)過(guò)度地平整更多的地形特征,造成錯(cuò)誤的分類(lèi). 由此可見(jiàn),采用固定大小的濾波窗口很難對(duì)點(diǎn)云進(jìn)行很好的分類(lèi),因此本文采用了線(xiàn)性漸進(jìn)增大的濾波窗口.

2.2 點(diǎn)云分類(lèi)

圖2為某一地區(qū)的地理影像,圖3為對(duì)應(yīng)的原始LiDAR點(diǎn)云高程值圖,區(qū)域中存在著大量的樹(shù)木,右半部分有一個(gè)明顯的河流空洞區(qū).區(qū)域大小約為694 m×2100 m,點(diǎn)云高程差約為172 m,點(diǎn)云間距約1.2 m,共有1 717 762個(gè)離散點(diǎn).

圖2 地理影像圖

圖3 原始LiDAR點(diǎn)云高程值圖

針對(duì)該區(qū)域的特點(diǎn),采用以下流程將區(qū)域中的物質(zhì)分為樹(shù)木、地面、水三種類(lèi)型.

1) 剔除噪點(diǎn). 采用k鄰近查詢(xún)球[11]方法,剔除極低點(diǎn)、孤立點(diǎn)等高程異常值.

2) 點(diǎn)云格網(wǎng)化. 使用一個(gè)二維格網(wǎng)覆蓋在LiDAR數(shù)據(jù)點(diǎn)上,格網(wǎng)間距通常略小于LiDAR點(diǎn)云密度,如果一個(gè)格網(wǎng)單元里有多個(gè)數(shù)據(jù)點(diǎn),取格網(wǎng)單元中高程的最小數(shù)據(jù)點(diǎn)代表該格網(wǎng)單元的值,若網(wǎng)格單元中沒(méi)有數(shù)據(jù)點(diǎn)的則不做處理[14].

3) 用二值圖像表示格網(wǎng)后的點(diǎn)云數(shù)據(jù),有值的網(wǎng)格單元高程值用“1”表示,否則用“0”表示. 用合適的“結(jié)構(gòu)元”對(duì)二值圖像執(zhí)行閉運(yùn)算. 經(jīng)過(guò)閉運(yùn)算后,小面積的空白區(qū)域?qū)?huì)被填充,剩下的空白就是大面積的水體造成空洞區(qū),并標(biāo)記.

4) 對(duì)非水體區(qū)域執(zhí)行數(shù)學(xué)形態(tài)學(xué)濾波,判斷濾波前后高程差閾值是否滿(mǎn)足條件,分離出樹(shù)木點(diǎn)和地面點(diǎn).

5) 由地面點(diǎn)通過(guò)線(xiàn)性插值獲得地形.

圖4和圖5分別為最終獲得的地形以及分類(lèi)結(jié)果.從圖4中可以看出,在整個(gè)區(qū)域中,樹(shù)木被很好地分離,而地形特征也保持相對(duì)圓滿(mǎn). 將圖5對(duì)照原始地理影像圖,可以看出,各種地物得到了很好的分類(lèi). 將二者結(jié)合起來(lái),對(duì)各種地物賦上相匹配的電磁參數(shù),就可以對(duì)上述地理環(huán)境進(jìn)行電磁建模.

圖4 地形圖

圖5 點(diǎn)云分類(lèi)結(jié)果

3 電磁環(huán)境仿真

利用三維拋物方程模型預(yù)測(cè)上節(jié)所述地理環(huán)境下的電波傳播特性.仿真參數(shù)設(shè)置如下:天線(xiàn)頻率為900 MHz;高斯方向圖3 dB帶寬為10°;水平極化波;天線(xiàn)架設(shè)在橫向距離為0 m,縱向距離為347 m,距離地表高度為45 m的地理環(huán)境處;大氣條件為理想大氣條件. 受限于數(shù)據(jù)源,最大傳播距離為2 100 m,dy=dz=0.5 m,傳播方向上的步長(zhǎng)取dx=1 m.

圖6顯示了源所在的地形剖面上傳播因子的分布情況,分別選取了(a)單一地表為黏性干土(介電常數(shù)為5,電導(dǎo)率為0.021 S/m)、(b) 真實(shí)的地理環(huán)境兩種情況,其中地理環(huán)境中水的相對(duì)介電常數(shù)和電導(dǎo)率為81和0.22 S/m,樹(shù)木的等效介電常數(shù)和電導(dǎo)率分別為1.004和180 μS/m[5]. 圖6(a)中由于地形的存在,在地形后方出現(xiàn)了電磁波的盲區(qū)以及繞射現(xiàn)象,同時(shí)也可以看到由于地形的反射出現(xiàn)了明顯的干涉條紋. 圖6(b)中由于樹(shù)木引起的遮擋,散射和吸收會(huì)產(chǎn)生較大的路徑損耗,導(dǎo)致空間電磁波的整體損耗要大于沒(méi)有樹(shù)木時(shí)的情形,特別是在森林的內(nèi)部,電磁波的損耗明顯增大.

(a) 地表為干土

(b) 真實(shí)地理環(huán)境下圖6 傳播因子分布偽彩圖

圖8 兩種算法下傳播因子隨高度的變化

圖7為采用FD算法[15]計(jì)算得到的上述地理環(huán)境中傳播因子的分布情況.與圖6(b)對(duì)比,可以看出兩者總體趨勢(shì)大致相同,采用三維拋物方程能夠考慮到電波的橫向繞射,更加真實(shí)地反應(yīng)出地理環(huán)境中電波傳播特性. 圖8顯示了主輻射方向上,傳播距離為2 100 m處,采用FD算法和SSFT算法得到的傳播因子隨高度的變化,兩者吻合較好,說(shuō)明了本文算法的正確性及有效性,SSFT算法相對(duì)于FD算法擁有更高的效率,能夠快速地預(yù)測(cè)環(huán)境中的電波傳播特性.

圖9中顯示的是主輻射方向上,傳播距離一定時(shí),有樹(shù)木和無(wú)樹(shù)木、有水無(wú)水的場(chǎng)景下傳播因子隨高度的變化情況.圖9(a)是傳播距離為2 000 m處,對(duì)比有樹(shù)木和無(wú)樹(shù)木情況下的傳播因子.可以看出:在高度小于50 m的樹(shù)木內(nèi)部,電波的衰減較大;在遠(yuǎn)離樹(shù)木的區(qū)域,電波的衰減相對(duì)較小. 圖9(b)是傳播距離為1 920 m處,有無(wú)水時(shí)傳播因子隨高度的變化曲線(xiàn).從圖中可以看出,兩者相差很小,其主要原因是,水體的范圍只有幾十米,遠(yuǎn)離輻射源,電磁波成掠入射接觸水體,并且水體前面存在著樹(shù)木的遮擋作用,總體上水體對(duì)電磁波的傳播影響較小.

(a) 傳播距離為2 000 m

(b) 傳播距離為1 920 m圖9 不同傳播距離處傳播因子隨高度的變化

圖10是在上述地理環(huán)境中,發(fā)射天線(xiàn)分別架設(shè)在距離地面高度35 m、45 m和55 m的情況下,距離發(fā)射天線(xiàn)2 100 m處的傳播因子隨高度的變化情況. 通過(guò)比較可以看出,在一定接收高度下,電波的衰減隨著發(fā)射天線(xiàn)架設(shè)高度的增加而減小,因此適當(dāng)把天線(xiàn)架設(shè)在高處,可以減小到達(dá)接收點(diǎn)的信號(hào)衰減,擴(kuò)大無(wú)線(xiàn)信號(hào)的覆蓋范圍.

圖10 不同天線(xiàn)高度下傳播因子隨高度的變化

4 結(jié) 論

本文研究了基于LiDAR點(diǎn)云的PEM電波傳播模型,針對(duì)傳統(tǒng)地理電磁環(huán)境建模的不足,研究了LiDAR點(diǎn)云數(shù)據(jù)中地物分類(lèi),對(duì)傳播區(qū)域的地理電磁環(huán)境進(jìn)行了更為準(zhǔn)確的建模,并采用拋物方程模型對(duì)有地表覆蓋物的環(huán)境下電波傳播特性進(jìn)行了仿真分析. 相對(duì)于傳統(tǒng)的電磁建模方法,本方法能夠考慮到地表覆蓋物對(duì)電波的影響,從而更真實(shí)地反映出電波的傳播特性. 下一步的工作重點(diǎn)是在現(xiàn)有的基礎(chǔ)上,研究包含建筑物、樹(shù)木的復(fù)雜城區(qū)的LiDAR點(diǎn)云分類(lèi),并結(jié)合大氣環(huán)境對(duì)區(qū)域中的電波傳播特性進(jìn)行預(yù)測(cè).

[1] DUMONT N, WATSON R J, PENNOCK S R. Use of the parabolic equation propagation model to predict TV white space availability[C]//Antennas and Propagation Conference. Loughborough, November 8-9, 2010: 353-356.

[2] LEONTOVICH M, FOCK V. Solution of the problem of propagation of electromagnetic waves along the earth’s surface by the method of parabolic equation[J]. Journal of physics USSR, 1946(7): 557-573.

[3] DONOHUE D J, KUTTLER J R. Propagation modeling over terrain using the parabolic wave equation[J]. IEEE transactions on antennas & propagation, 2000, 48(2): 260-277.

[4] HOLM P, WAERN A. Wave propagation over a forest edge parabolic equation modelling vs. GTD modelling[C]//IEEE International Symposium on Electromagnetic Compatibility. Istanbul, May 11-16, 2003: 764-767.

[5] 郭建炎, 王劍瑩, 龍?jiān)屏? 等. 基于拋物方程法的部分森林覆蓋山區(qū)電波傳播分析[J]. 電波科學(xué)學(xué)報(bào), 2008, 23(6): 1045-1050.

GUO J Y, WANG J Y, LONG Y L, et al. Analysis of radio propagation in partly forested terrain environment using parabolic equation approach[J]. Chinese journal of radio science, 2008, 23(6): 1045-1050. (in Chinese)

[6] 張青洪, 廖成, 盛楠, 等. 森林環(huán)境電波傳播拋物方程模型的改進(jìn)研究[J]. 物理學(xué)報(bào), 2013, 62(20): 204101.

ZHANG Q H, LIAO C, SHENG N, et al. Improved study on parabolic equation model for radio wave propagation in forest[J]. Acta physica sinica, 2013, 62(20):204101. (in Chinese)

[7] 白瑞杰, 廖成, 張青洪,等. 復(fù)雜地理環(huán)境的圖像分割及電波傳播特性[J]. 強(qiáng)激光與粒子束, 2015, 27(10): 63-67.

BAI R J, LIAO C, ZHANG Q H, et al. Image segmentation of complex geographical environment and wave propagation characteristics[J]. High power laser and particle beams, 2015, 27(10): 63-67. (in Chinese)

[8] 陳松堯, 程新文. 機(jī)載LIDAR系統(tǒng)原理及應(yīng)用綜述[J]. 測(cè)繪工程, 2007, 16(1): 27-31.

CHEN S Y, CHENG X W. The principle and application of airborne LIDAR[J]. Engineering of surveying and mapping, 2007, 16(1): 27-31. (in Chinese)

[9] LEVY M. Parabolic equation methods for electromagnetic wave propagation[M]. London: IEE Press, 2000: 35-40.

[10] HARDIN R H, TAPPERT F D. Application of the split-step Fourier method to the numerical solution of nonlinear and variable coefficient wave equations[J]. SIAM review, 1973, 15(423): 423.

[11] 成曉倩, 馬洪超, 趙紅強(qiáng),等. 一種基于區(qū)域增長(zhǎng)的機(jī)載LIDAR濾波[J]. 測(cè)繪科學(xué), 2010, 35(1): 61-63.

CHEN X Q, MA H C, ZHAO H Q, et al. A filtering of LIDAR based on region growing[J]. Science of surveying and mapping, 2010, 35(1): 61-63. (in Chinese)

[12] ZHANG W, QI J, WAN P, et al. An easy-to-use airborne LiDAR data filtering method based on cloth simulation[J]. Remote sensing, 2016, 8(6): 501.

[13] PINGEL T J, CLARKE K C, MCBRIDE W A. An improved simple morphological filter for the terrain classification of airborne LIDAR data[J]. ISPRS journal of photogrammetry & remote sensing, 2013, 77(1): 21-30.

[14] ZHANG K, CHEN S C, WHITMAN D, et al. A progressive morphological filter for removing nonground measurements from airborne LIDAR data[J]. IEEE transactions on geoscience & remote sensing, 2003, 41(4): 872-882.

[15] DING W, WANG K, LONG Y. Study of the electromagnetic waves propagation over the improved fractal sea surface based on parabolic equation method[J]. International journal of antennas and propagation, 2016(3): 1-7.