本發(fā)明屬于電波傳播技術(shù)領(lǐng)域，具體涉及一種高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法。

背景技術(shù)：

低頻電波以其波長較長，信號的傳播損耗小，信號的幅度與相位穩(wěn)定，被廣泛應(yīng)用于授時、導(dǎo)航、通信等領(lǐng)域。為提高低頻無線電工程的使用性能和精度，需要深入研究低頻電波在各種地域、時域、頻域范圍內(nèi)的變化規(guī)律及預(yù)測技術(shù)。目前，預(yù)測復(fù)雜環(huán)境下低頻電波傳播已有的理論方法有：積分方程方法、拋物方程方法和時域有限差分(finite-differencetime-domain，fdtd)方法。積分方程方法和拋物方程方法都可以考慮地形變換，對復(fù)雜路徑具有較高精度，但是，這兩種方法都忽略了沿路徑后向傳播的反射場影響，在地形起伏劇烈的傳播路徑上會引起較大誤差。fdtd方法可方便的考慮任何地質(zhì)和地形路徑的電波傳播問題，并且精度較高，但其存在過長的計算時間、較多的內(nèi)存資源占用等問題，不適合大區(qū)域的工程化推廣。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的是提供一種高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，既能解決fdtd方法計算量大、耗時長的缺點，又能解決積分方程方法和拋物方程方法在地形起伏劇烈路徑中由于忽略后向電波傳播影響引起的誤差大的問題。

本發(fā)明所采用的技術(shù)方案是：高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，具體按照以下步驟實施：

步驟1：輸入模型文件；

步驟2：利用平地面公式計算初始位置ρ＝ρ0處縱向橫切面上的磁場hf(ρ0,z)，并通過其求解雙向拋物方程的前向初始場uf(ρ0,z)；

步驟3：結(jié)合步驟2計算出來的前向初始場uf(ρ0,z)，基于坐標(biāo)變換模型，采用分布離散傅里葉變換算法，求解計算區(qū)域任意位置電波傳播的前向場uf(ρ,z)；

步驟4：結(jié)合步驟2計算出來的前向初始場uf(ρ0,z)，基于階梯近似模型，采用ssft算法，遞歸的求解由于地形影響對電波傳播產(chǎn)生的后向場和由多次反射產(chǎn)生的前向場τ＝1,2,3,…；

步驟5：結(jié)合步驟3和步驟4中的傳播場結(jié)果uf(ρ,z)、τ＝1,2,3,…和求出電波傳播的總磁場

本發(fā)明的特點還在于：

步驟1具體為：

計算區(qū)域大小nρ×nz，其中nρ為ρ方向的網(wǎng)格數(shù)，nz為z方向的網(wǎng)格數(shù)；空間網(wǎng)格步長分別為δρ和δz，ρ和z分別表示橫向和縱向坐標(biāo)；ρ0為初始場位置；源的參數(shù)；地面相對介電常數(shù)εr和電導(dǎo)率σ；地形模型參數(shù)。

步驟2具體為：

選取二維柱坐標(biāo)系(ρ,z)，其中ρ和z分別表示為橫向和縱向坐標(biāo)，根據(jù)實際發(fā)射天線尺寸，通過測量得到垂直電偶極子的電荷間距dl，放置在距離地面高度為d的位置，假設(shè)時諧因子為e^-iωt，利用平地面公式計算初始距離ρ＝ρ0處縱向橫切面的磁場hf(ρ0,z)：

其中，i為電流大小，k0和kg分別為真空和地面的波數(shù)，r0表示從源點到觀測點的直線距離，r1表示從源的鏡像到觀測點的直線距離，p2為中間參量：

f(z)是fresnel積分并定義為

定義沿ρ軸正方向傳播的波函數(shù)，即前向場uf(ρ0,z)為：

步驟3具體為：

在原坐標(biāo)系(ρ,z)中，地形函數(shù)為t(ρ)，以地形表面上任意一點為新的坐標(biāo)原點o′建立新的坐標(biāo)系則

在新坐標(biāo)系中，波函數(shù)已經(jīng)修正為：

其中，

采用起伏地形ssft算法求解下一步進處的前向場：

其中，n為大氣折射率，和為離散混合傅里葉變換對，p＝k0sinθ，θ為電波傳播的仰角；

步驟2中得到初始距離ρ＝ρ0處電波傳播的前向初始場為uf(ρ0,z)，結(jié)合式(4)和(5)，即可得到新坐標(biāo)系下初始場

將代入式(6)中步進計算，得到新坐標(biāo)系下任意位置處的前向場

再考慮原坐標(biāo)系和新坐標(biāo)系前向場之間的關(guān)系，即式(4)，求出原坐標(biāo)系中前向場uf(ρ,z)。

步驟4具體為：

沿ρ軸正方向傳播的各個前向場為：

其中，上標(biāo)τ表示前向場的序號，τ＝0表示電波傳播的主傳播場，而τ＝1,2,3…則表示由后向場多次反射所產(chǎn)生的前向場，為前向磁場分量；

當(dāng)電波沿ρ軸正方向傳播時，采用ssft算法即可求得下一位置處的前向場：

將步驟2中求得的前向初始場uf(ρ0,z)賦值給若傳播路徑不含起伏地形，則計算區(qū)域僅存在τ＝0時對應(yīng)的將代入式(8)依次迭代計算出各個位置的場若存在起伏地形，需要考慮后向場影響；

定義沿ρ軸負方向傳播的各個后向波函數(shù)，即后向場：

其中，表示后向場的序號，每當(dāng)前向場遇到一個相對于其傳播方向的上升階梯面，就會產(chǎn)生一個后向磁場分量

對應(yīng)于式(9)的ssft算法為：

假設(shè)前向波沿ρ軸正方向傳播到ρ＝ρe處存在相對于傳播方向的上升階梯面，由于該上升階梯面對電波傳播的反射作用，前向波分為兩部分：繼續(xù)向前傳播的前向波和反射作用產(chǎn)生的后向波；

在上升階梯面上，根據(jù)磁場所滿足切向邊界條件可以得到后向初始場為：

修正前向場：

其中，τ和都是具體的某一個標(biāo)號，具體過程為標(biāo)記為τ的前向波，每遇到一個相對于ρ正方向的上升沿就會產(chǎn)生一個標(biāo)記為的后向波，并且更新為

假設(shè)電波沿ρ軸負方向傳播到ρ＝ρt處存在相對于傳播方向的上升階梯面，后向波分為兩部分：繼續(xù)向后傳播的后向波和由后向波反射產(chǎn)生的前向波；

在上升階梯面上，根據(jù)磁場所滿足切向邊界條件可以得到后向場多次反射產(chǎn)生的前向場為：

修正后向場：

其中，τ和都是具體的某一個標(biāo)號。具體過程為標(biāo)記為的后向波，每遇到一個相對于ρ負方向的上升沿就會產(chǎn)生一個標(biāo)記為τ的前向波，并且τ更新為τ＝τ+1；

轉(zhuǎn)到式(7)和式(9)繼續(xù)循環(huán)對所有前向場和后向場迭代計算下一步進處的場強，直至步進到計算區(qū)域邊界；其中，前向場的邊界為ρmax＝nρδρ，后向場的邊界為ρ0。

步驟4中循環(huán)時設(shè)置門限終止計算，循環(huán)計算的判定依據(jù)為：

其中，u(ρ,z)表示τ＝1,2,3…或者判決系數(shù)ξ設(shè)置為0.01，||||為矢量求模。

步驟5具體為：

將步驟3和步驟4中的前向場uf(ρ,z)和多次反射產(chǎn)生的前向場τ＝1,2,3,…分別代入式(2)和式(7)，分別求出所有存在的前向場hf(ρ,z)和

將步驟4計算的代入式(9)中，求出所有存在的后向磁場分量

電波傳播的總磁場為所有前向磁場分量和后向磁場分量的疊加：

本發(fā)明的有益效果是：

①本發(fā)明高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，是基于柱坐標(biāo)系實現(xiàn)的，因此，能夠更為方便的處理圓柱對稱的電波傳播問題，對工程指導(dǎo)比較有意義；

②本發(fā)明高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，能夠準(zhǔn)確的預(yù)測復(fù)雜環(huán)境中電波傳播的前向和后向傳播影響，克服了積分方程方法和傳統(tǒng)的拋物方程方法忽略后向波影響引起的誤差。此外，與fdtd方法相比，在相同的計算精度前提下，該方法基于高效的ssft算法步進求解場強，大大的縮短了計算時間，更適合工程推廣。

附圖說明

圖1是本發(fā)明方法的流程圖；

圖2(a)是高度為0.5km時本發(fā)明方法與fdtd方法的地面場強比較；

圖2(b)是高度為1.5km時本發(fā)明方法與fdtd方法的地面場強比較；

圖3是在多個高斯山脈地形情況下，本發(fā)明方法與fdtd方法的地面場強比較圖；

圖4(a)是fdtd方法的空間場強偽彩色圖比較；

圖4(b)是本發(fā)明方法的空間場強偽彩色圖比較。

具體實施方式

下面結(jié)合附圖和具體實施方式對本發(fā)明進行詳細說明。

本發(fā)明提供了高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，改進了傳統(tǒng)拋物方程方法在復(fù)雜地形傳播路徑上由于忽略后向傳播影響引起的預(yù)測誤差，具體實現(xiàn)過程是：首先利用平地面公式計算出前向初始場；接著，基于坐標(biāo)變化模型，采用拋物方程方法預(yù)測整個區(qū)域在忽略后向傳播影響情況下的前向場的場值；然后，基于階梯近似模型，通過雙向拋物方程方法預(yù)測由于電波在復(fù)雜地形作用下引起的反射場(包括后向場和多次反射產(chǎn)生的前向場)影響。最后，通過拋物方程方法計算的前向場和雙向拋物方程方法計算的后向場和多次反射產(chǎn)生的前向場疊加求解總磁場，具體實施步驟流程圖見圖1。

在采用雙向拋物方程方法求圓柱對稱模型問題中的前向場和后向場時，首先需要推導(dǎo)柱坐標(biāo)系下的雙向拋物方程及其對應(yīng)的ssft算法：

假設(shè)電磁波的時諧因子是e^-iωt，取二維的柱坐標(biāo)系(ρ,z)，其中ρ和z分別表示為橫向和縱向坐標(biāo)，并以標(biāo)量φ表示任意標(biāo)量場分量，假設(shè)φ與方位角無關(guān)，則對于水平極化波，而對于垂直極化波，φ滿足二維波動方程

其中k0是真空中的波數(shù)，n為大氣折射率。首先，作如下變量代換：

其中，ψ(ρ,z)僅考慮了電磁波的柱面?zhèn)鞑ビ绊?，沒有涉及相位修正。將式(18)代入(17)中，可以得到：

然后，對上式作遠場近似，即假設(shè)ρ>>1且k0ρ>>1，式(19)變?yōu)椋?/p>

當(dāng)折射指數(shù)n隨ρ的變化很小時，式(20)可分解為

由于式(21)中的兩個因式互不相關(guān)，故可以得到關(guān)于ρ的兩個獨立的式子

ψf和ψb分別對應(yīng)于前向波和后向波。接下來，考慮相位修正，引入合適的相位因子，對ψf和ψb分別作如下修正：

其中，φf和φb分別代表前向和后向場分量。將式(18)代入式(24)和(25)中，可以分別得到

其中，φf和φb分別為前向和后向磁場分量。

下面將式(24)和(25)分別代入(22)和(23)，可以得到雙向拋物方程：

上面兩式分別稱為前向和后向拋物方程。采用ssft算法求解上面兩式，可得：

式中和為離散混合傅里葉變換對，為大氣折射率，δρ為ρ方向的網(wǎng)格大小。實際計算中，φf和φb可能分別為一系列的和其中τ＝0,1,2,3,…和

本發(fā)明高精度預(yù)測低頻電波傳播特性的雙向拋物方程方法，具體按照以下步驟實施：

步驟1：輸入模型文件，具體為：

計算區(qū)域大小nρ×nz，其中nρ為ρ方向的網(wǎng)格數(shù)，nz為z方向的網(wǎng)格數(shù)；空間網(wǎng)格步長分別為δρ和δz，ρ和z分別表示橫向和縱向坐標(biāo)；ρ0為初始場位置；源的參數(shù)；地面相對介電常數(shù)εr和電導(dǎo)率σ；地形模型參數(shù)。

步驟2：利用平地面公式計算初始位置ρ＝ρ0處縱向橫切面上的磁場hf(ρ0,z)，并通過其求解雙向拋物方程的前向初始場uf(ρ0,z)，具體為：

為了更好地解決圓柱對稱模型問題，在此選取二維柱坐標(biāo)系(ρ,z)，其中ρ和z分別表示為橫向和縱向坐標(biāo)，根據(jù)實際發(fā)射天線尺寸，通過測量得到垂直電偶極子的電荷間距dl，放置在距離地面高度為d的位置，假設(shè)時諧因子為e^-iωt，利用平地面公式計算初始距離ρ＝ρ0處縱向橫切面的磁場hf(ρ0,z)：

其中，i為電流大小，k0和kg分別為真空和地面的波數(shù)，r0表示從源點到觀測點的直線距離，r1表示從源的鏡像到觀測點的直線距離，p2為中間參量：

f(z)是fresnel積分并定義為

定義沿ρ軸正方向傳播的波函數(shù)，即前向場uf(ρ0,z)為：

步驟3：結(jié)合步驟2計算出來的前向初始場uf(ρ0,z)，基于坐標(biāo)變換模型，采用分布離散傅里葉變換(split-stepfouriertransform，ssft)算法，求解計算區(qū)域任意位置電波傳播的前向場uf(ρ,z)，具體為：

在原坐標(biāo)系(ρ,z)中，地形函數(shù)為t(ρ)，以地形表面上任意一點為新的坐標(biāo)原點o′建立新的坐標(biāo)系則

在新坐標(biāo)系中，波函數(shù)已經(jīng)修正為：

其中，

采用起伏地形ssft算法求解下一步進處的前向場：

其中，n為大氣折射率，和為離散混合傅里葉變換對，p＝k0sinθ，θ為電波傳播的仰角；

步驟2中得到初始距離ρ＝ρ0處電波傳播的前向初始場為uf(ρ0,z)，結(jié)合式(4)和(5)，即可得到新坐標(biāo)系下初始場

將代入式(6)中步進計算，得到新坐標(biāo)系下任意位置處的前向場

再考慮原坐標(biāo)系和新坐標(biāo)系前向場之間的關(guān)系，即式(4)，求出原坐標(biāo)系中前向場uf(ρ,z)。

步驟4：結(jié)合步驟2計算出來的前向場uf(ρ0,z)，基于階梯近似模型，采用ssft算法，遞歸的求解由于地形影響對電波傳播產(chǎn)生的后向場和由多次反射產(chǎn)生的前向場τ＝1,2,3,…，具體為：

沿ρ軸正方向傳播的各個前向場為：

其中，上標(biāo)τ表示前向場的序號，τ＝0表示電波傳播的主傳播場，而τ＝1,2,3…則表示由后向場多次反射所產(chǎn)生的前向場，為前向磁場分量；與式(2)相比，可以看出，前向場的磁場求解具有相同的表達形式。

當(dāng)電波沿ρ軸正方向傳播時，采用ssft算法即可求得下一位置處的前向場：

將步驟2中求得的前向初始場uf(ρ0,z)賦值給若傳播路徑不含起伏地形，則計算區(qū)域僅存在τ＝0時對應(yīng)的將代入式(8)依次迭代計算出各個位置的場若存在起伏地形，就需要考慮后向場影響；

下面定義沿ρ軸負方向傳播的各個后向波函數(shù)，即后向場：

其中，表示后向場的序號，每當(dāng)前向場遇到一個相對于其傳播方向的上升階梯面，就會產(chǎn)生一個后向磁場分量

對應(yīng)于式(9)的ssft算法為：

比較式(8)和(10)可以得出，前向場和后向場的ssft算法迭代形式基本一致。

假設(shè)前向波沿ρ軸正方向傳播到ρ＝ρe處存在相對于傳播方向的上升階梯面，由于該上升階梯面對電波傳播的反射作用，前向波分為兩部分：繼續(xù)向前傳播的前向波和反射作用產(chǎn)生的后向波。在上升階梯面上，根據(jù)磁場所滿足切向邊界條件可以得到后向初始場為：

修正前向場：

這里的τ和都是具體的某一個標(biāo)號。具體過程為標(biāo)記為τ的前向波，每遇到一個相對于ρ正方向的上升沿就會產(chǎn)生一個標(biāo)記為的后向波，并且更新為

同理，假設(shè)電波沿ρ軸負方向傳播到ρ＝ρt處存在相對于傳播方向的上升階梯面，后向波分為兩部分：繼續(xù)向后傳播的后向波和由后向波反射產(chǎn)生的前向波。在上升階梯面上，根據(jù)磁場所滿足切向邊界條件可以得到后向場多次反射產(chǎn)生的前向場為：

修正后向場：

這里的τ和都是具體的某一個標(biāo)號。具體過程為標(biāo)記為的后向波，每遇到一個相對于ρ負方向的上升沿就會產(chǎn)生一個標(biāo)記為τ的前向波，并且τ更新為τ＝τ+1。

然后，采用式(7)和(9)繼續(xù)對所有前向場和后向場迭代計算下一步進處的場強，接著判斷是否存相對于傳播方向上的上升階梯面，具體執(zhí)行過程見上述假設(shè)，循環(huán)此過程直至步進到計算區(qū)域邊界(對于前向場，邊界為ρmax＝nρδρ，對于后向場，邊界為ρ0)。

由此可以看出對于孤立山峰，電波向前傳播遇到山峰反射會產(chǎn)生許多后向波，而后向波不存在反射，即不會產(chǎn)生前向波；而對于多個山峰，當(dāng)電波向前傳播遇到第一個山峰時，一部分電波會由于反射作用而產(chǎn)生后向波，另一部分電波則會繼續(xù)向前傳播，當(dāng)遇到第二個山峰時，又會出現(xiàn)后向波，并且，后向波會在兩個山峰之間往復(fù)反射，以此類推，從而產(chǎn)生一系列的反射波(其中包括多次反射產(chǎn)生的前向波和后向波)。由于受空間散射及介質(zhì)的吸收影響，循環(huán)計算的反射波場強逐漸減小，繼而不會對全域場產(chǎn)生影響。因此，有必要設(shè)置門限來終止計算。為了方便計算，給出遞歸計算的判定依據(jù)為：

其中，u(ρ,z)表示τ＝1,2,3…或者判決系數(shù)ξ設(shè)置為0.01，||||為矢量求模。

步驟5：結(jié)合步驟3和步驟4中的傳播場結(jié)果uf(ρ,z)、和τ＝1,2,3,…，求出電波傳播的總磁場具體為：

將步驟3和步驟4中的前向場uf(ρ,z)和多次反射產(chǎn)生的前向場τ＝1,2,3,…分別代入式(2)和式(7)，分別求出所有存在的前向場hf(ρ,z)和注意：這里沒有采用而是用步驟3求出來的hf(ρ,z)替代，這是因為采用階梯近似模型求出來的前向場分量精度在山后不夠。

將步驟4計算的代入式(9)中，求出所有存在的后向磁場分量

電波傳播的總磁場為所有前向磁場分量和后向磁場分量的疊加：

實施例1

單個高斯山脈地形地面場強預(yù)測

垂直電偶極子天線的輻射功率均采用1kw，信號源頻率為100khz。計算區(qū)域總大小為ρmax：100km×zmax：102.4km，網(wǎng)格剖分大小均分別為dρ＝200m和dz＝100m，初始距離ρ0＝10km；地面電參數(shù)為εr＝13，σ＝3×10^-3s/m(陸地)；在中心位置ρc＝50km處有一孤立高斯山峰，其高度函數(shù)為

l為山峰寬度取2km，h為山峰高度且分別取0.5km和1.5km。圖2(a)和圖2(b)分別為在相同寬度不同高度的單個高斯山脈地形情況下，本發(fā)明方法與fdtd方法的地面場強比較。由圖2(a)和圖2(b)可以看出，兩種方法的計算結(jié)果吻合一致，并且都可以預(yù)測到山前電波反射影響的振蕩存在。但是，相比于fdtd方法，雙向拋物方程方法計算更快。

實施例2

多個高斯山脈地形地面場強預(yù)測

將實施例1中的山峰改為多個高斯山脈地形，第一個山峰的高度為1km，寬度為4km，中心位置為40km；第二個山峰的高度為1.5km，寬度為2km，中心位置為60km，其他參數(shù)不變。圖3給出了在多個高斯山脈地形情況下，本發(fā)明方法與fdtd方法的地面場強比較。由圖3可以看出，對于多個高斯山脈地形，雙向拋物方程方法也能精確的預(yù)測山峰之間波的往復(fù)反射影響。并且經(jīng)統(tǒng)計，fdtd方法的計算時間是雙向拋物方程方法的21倍。圖4(a)和圖4(b)分別為fdtd方法與本發(fā)明方法的空間場強偽彩色圖比較，從圖中可以直觀的看到由于山峰的影響，電波在山前和山峰之間的反射影響。