粒子群優(yōu)化算法反演無線電波導(dǎo)傳輸損耗分布

龐佳瑋 杜曉燕 張水蓮 江長蔭

(1.信息工程大學(xué)，河南 鄭州 450002；2.中國電波傳播研究所，河南 新鄉(xiāng) 453003)

引 言

無線電波導(dǎo)現(xiàn)象產(chǎn)生于對流層內(nèi)大氣折射率的異常變化，是海洋大氣環(huán)境的一種常見現(xiàn)象.修正折射率在垂直方向的逆變化可以使電磁波的傳播路徑發(fā)生改變，導(dǎo)致電磁波陷獲.全球范圍內(nèi)，沿海地區(qū)存在著豐富的無線電波導(dǎo)[1-2].我國海域是無線電波導(dǎo)發(fā)生的高概率區(qū)域(某些區(qū)域蒸發(fā)波導(dǎo)的發(fā)生概率達(dá)到了70%[3])，這就為無線電波導(dǎo)條件下實(shí)現(xiàn)電磁波的超視距傳播提供了有利條件.

目前，研究無線電波導(dǎo)傳輸特性主要是利用基于雷達(dá)海雜波反演大氣折射率(Refractivity From Clutter，RFC).目前公開的報(bào)道表明,該技術(shù)多是反演蒸發(fā)波導(dǎo)，所涉及的設(shè)備相對復(fù)雜，發(fā)射功率大.同時(shí)，海雜波出現(xiàn)的方式多數(shù)情況下存在著不確定性[4]，海表面散射截面帶來的不確定性也極大限制了其反演的精度，且其以修正折射率垂直剖面作為反演目標(biāo).而實(shí)際中的很多情況下(如鏈路設(shè)計(jì)、站址選擇等)，人們更關(guān)心的是無線電波導(dǎo)條件下的電波傳輸損耗分布.通信信號發(fā)射天線有時(shí)因地理?xiàng)l件限制，處于較高的位置，其電波實(shí)現(xiàn)超視距傳播所涉及的波導(dǎo)類型也較為復(fù)雜[5].本文基于通信鏈路，探討如何反演空間中基本傳輸損耗的分布.

1 研究方法

本節(jié)主要介紹該方法的基本思路.

首先，接收端處在電波傳播的某超視距距離處接收通信信號，利用收發(fā)通信信號功率得到該位置處的基本傳輸損耗.然后，選擇合適的波導(dǎo)模型和優(yōu)化方法，通過不斷改變修正折射率模型參數(shù)，并利用發(fā)收站基本信息(基本傳輸損耗、天線基本參數(shù)、地理位置信息等)，實(shí)現(xiàn)該位置處的電波基本傳輸損耗計(jì)算值Lbr與目標(biāo)值Lbo的擬合，從而實(shí)現(xiàn)對高度—傳播距離這個(gè)二維空間上基本傳輸損耗分布的反演.基本思路如圖1(見 20 頁)所示.

2 相關(guān)模型與算法

2.1 無線電波導(dǎo)垂直剖面模型

無線電波導(dǎo)通常被特征化為依賴于高度和距離的環(huán)境折射率.垂直方向上，表面波導(dǎo)修正折射率剖面一般通過三線模型近似，而蒸發(fā)波導(dǎo)可以通過指數(shù)函數(shù)來近似建模.考慮到實(shí)際通信鏈路中涉及到的波導(dǎo)類型復(fù)雜，本文采用五參數(shù)的混合波導(dǎo)模型來對波導(dǎo)建模.該模型的大氣修正折射率垂直分布關(guān)系為[5]

(1)

式中：z0是粗糙度因子(通常取0.000 15)；c0為中性層結(jié)蒸發(fā)波導(dǎo)參數(shù)(通常取0.13)；m3=0.118 M-unit s/m；M0為混合層斜率和z=0相截的值(取為340)；表征折射率垂直方向變化的5個(gè)主要參數(shù)如下：Md為逆變層逆差；zthick為表面波導(dǎo)的厚度；m1為混合層斜率；h1逆變層底高；hd為蒸發(fā)波導(dǎo)高度.五參數(shù)模型如圖2所示.

圖2 五參數(shù)經(jīng)驗(yàn)?zāi)Ｐ?/p>

Gerstoft等人指出大氣折射率水平微小變化會對空間場強(qiáng)產(chǎn)生較大影響，并指出這種水平非均勻性可以用馬爾科夫過程來描述[6-7].圖3為利用該理論仿真得到的逆變層波導(dǎo)底高h(yuǎn)1在100 km尺度上變化幅度控制在±20 m內(nèi)，同時(shí)又保持了良好的隨機(jī)性，與折射率水平方向上的微小隨機(jī)的變化特點(diǎn)是相符的.而混合波導(dǎo)的其余4個(gè)參數(shù)在距離尺度上保持不變[7].

圖3 逆變層底高隨距離變化的模型

首先，利用計(jì)算機(jī)模擬實(shí)現(xiàn)106個(gè)馬爾科夫過程樣本，每隔1 km利用高斯隨機(jī)數(shù)對h1更新一次，如圖3(a)所示為其中20個(gè)樣本.

其次，K-L變換是最優(yōu)正交線性變換，通過特征降維，可以提取隨機(jī)過程中能量最大的特征向量和對應(yīng)特征值并進(jìn)行主成分分析.利用該方法提取出馬爾科夫過程的5個(gè)最主要的主特征向量eigenvectori和對應(yīng)系數(shù)[1].利用5個(gè)系數(shù)來表征折射率在水平方向的隨機(jī)變化趨勢.圖3(b)為前10個(gè)特征值，圖3(c)為5個(gè)主要的特征向量.圖3(d)為基于5個(gè)主特征向量和相應(yīng)特征值疊加組合的20個(gè)馬爾科夫過程樣本.對比圖3(a)可以看出，利用5個(gè)主特征向量和對應(yīng)特征值擬合所得曲線可以準(zhǔn)確地描述圖3(a)曲線的變化趨勢.因而可以用5個(gè)水平參數(shù)coei,i=1,2,3,4,5，與5個(gè)特征向量重新擬合得到h1的馬爾科夫過程.h1在水平方向上的變化，可以表示為[7]

(2)

此處5個(gè)水平參數(shù)的變化范圍應(yīng)小于對應(yīng)特征值的方根，這樣才能與大氣折射率的微小變化相符[7].至此，通過10個(gè)參數(shù)(水平方向上的5個(gè)參數(shù)和垂直方向上的5個(gè)參數(shù))即可描述大氣修正折射率剖面的二維模型.選取馬爾科夫過程樣本，由水平參數(shù)coei即可得到一組h1的水平方向數(shù)據(jù)，圖4為由此構(gòu)建出的修正折射率分布圖.

2.2 拋物方程法

通過地球展平變換和保角變化，結(jié)合尺度分析，可將Helmholtz方程簡化為電磁波在大氣中傳播的拋物型方程(Parabolic Equation，PE)

圖4 修正折射率分布示意圖

(3)

式中：u為電場或磁場的函數(shù)(依賴于輻射場的極化方式)；m為大氣修正折射指數(shù)；x和z分別為電磁波傳播的距離和高度.采用拋物型方程法，結(jié)合上述無線電波導(dǎo)垂直剖面模型即可計(jì)算電磁波基本傳輸損耗.

假設(shè)下邊界z=0是平滑的海面，由天線參數(shù)確定初始場，就可以利用分步傅里葉數(shù)值解法的拋物型方程法(FFT-PE)循環(huán)計(jì)算得到全空間的電波強(qiáng)度[7-8].解的形式為

F[u(x,z)]}，

(4)

式中： Δx是選擇的步長；k是真空中的波數(shù)；F(·)是傅里葉正弦變換；F-1(·)是傅里葉正弦逆變換.

眾所周知，若忽略介質(zhì)對電磁波吸收、散射等作用引起的損耗，實(shí)際中的電波傳播傳輸損耗為

L= 32.45+20lgdkm+20lgfMHz

-A-Gt-Gr，

(5)

Lb=32.45+20lgdkm+20lgfMHz-A.

(6)

本文主要研究基本傳輸損耗.

2.3 反演方法

由于波導(dǎo)參數(shù)與基本傳輸損耗間存在著復(fù)雜的非線性關(guān)系，這就要求反演算法必須能夠克服局部極值影響，并適應(yīng)復(fù)雜的非線性特性.

目前現(xiàn)有的線性化反演方法，其反演結(jié)果的可靠性嚴(yán)重取決于初始模型，當(dāng)初始模型的選取遠(yuǎn)離真實(shí)模型時(shí)將很容易陷入局部極小值中得出錯(cuò)誤的求解；同時(shí)該方法還需要計(jì)算偏導(dǎo)數(shù)信息，雅可比矩陣偏導(dǎo)數(shù)的求取精度也將直接影響反演結(jié)果的質(zhì)量.現(xiàn)有的非線性反演優(yōu)化算法，如遺傳算法(Genetic Algorithm,GA)、模擬退火(Simulated Annealing, SA)、人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network, ANN)等，大大減小了對初始模型的依賴程度，用于無線電波導(dǎo)反演是相當(dāng)有用的.但由于它們編程相對復(fù)雜、參數(shù)難于控制、計(jì)算工作量大、收斂速度慢、計(jì)算時(shí)間長，在實(shí)際反演應(yīng)用中受到了限制[9-10].最近發(fā)展的貝葉斯濾波的方法理論相對復(fù)雜，實(shí)現(xiàn)較難.

粒子群(Particle Swarm Optimization, PSO)算法兼具全局和局部尋優(yōu)能力，操作簡單，收斂較快[10].粒子在自身記憶、個(gè)體經(jīng)驗(yàn)和社會經(jīng)驗(yàn)的共同作用下，逐漸收斂于目標(biāo)函數(shù)的全局最優(yōu)位置.其位置和速度更新公式為

(7)

式中：ω為慣性權(quán)重，權(quán)衡算法局部最優(yōu)與全局最優(yōu)的能力；r1和r2為[0,1]之間的隨機(jī)數(shù)，用來保持群體多樣性；c1和c2為學(xué)習(xí)因子，表示粒子自我總結(jié)和向群體中優(yōu)秀個(gè)體學(xué)習(xí)的能力.

實(shí)際應(yīng)用中合理設(shè)置慣性權(quán)重、最大速度、種群數(shù)和搜索次數(shù)，可以得到較好的反演效果.由此，建立關(guān)于接收位置基本傳輸損耗計(jì)算值Lbr(x)與目標(biāo)值Lbo的適應(yīng)度函數(shù)為

fitness=|Lbo-Lbr(x)|.

(8)

式中，x為波導(dǎo)參數(shù)組成的向量，分別是修正折射率剖面垂直方向的5個(gè)參數(shù)和水平方向的5個(gè)參數(shù).利用粒子群算法進(jìn)行損耗計(jì)算值與正演值的尋優(yōu)搜索，并確定波導(dǎo)參數(shù)，最終得出基本傳輸損耗空間分布.本文使用的PSO算法流程如圖5所示.

圖5 粒子群算法流程

3 仿真分析

3.1 仿真條件和主要步驟

發(fā)射天線和大氣邊界參數(shù)如表1所示，反演位置選取在距離發(fā)射天線90 km處的超視距位置，反演具體步驟如下：

步驟1目標(biāo)值Lbo的計(jì)算

本文采用文獻(xiàn)[1]中的實(shí)測大氣修正折射率垂直廓線(如圖6(a)所示).下邊界條件和發(fā)射天線基本參數(shù)如表1所示.利用2.1節(jié)中產(chǎn)生水平非均勻性的方法以及2.2節(jié)中的FFT-PE方法可計(jì)算出電波基本傳輸損耗分布，將反演位置的損耗值作為目標(biāo)值，如圖6(b)所示，以便于驗(yàn)證所提方法的可行性.

表1 發(fā)射天線和環(huán)境參數(shù)

(a) 實(shí)測大氣修正折射率垂直分布圖6 正演方法

步驟2反演計(jì)算

假設(shè)圖6(a)中波導(dǎo)模型參數(shù)未知，已知收發(fā)天線基本參數(shù)和環(huán)境參數(shù)，以及空間中某點(diǎn)的基本傳輸損耗(實(shí)際中傳輸損耗可由收/發(fā)站的收/發(fā)功率得到).本文中，在圖6(b)(見20頁)90 km處選擇幾個(gè)不同高度(分別為5 m，88.5 m，99 m，123 m，160 m)反演位置的基本傳輸損耗，利用粒子群算法反演一定范圍內(nèi)高度—距離二維空間中的基本傳輸損耗分布.

3.2 反演結(jié)果分析

在相同無線電波導(dǎo)和發(fā)射天線情況下，分別利用傳播距離上90 km處的幾個(gè)不同高度反演位置的基本傳輸損耗進(jìn)行該方法的反演計(jì)算，其中，最小損耗位于88.5 m高度處.每種情況均進(jìn)行20次反演計(jì)算，并對反演結(jié)果求均值.圖7為不同高度反演位置反演所得90 km處的垂直方向基本傳輸損耗分布與正演結(jié)果的比較.圖8為88.5 m高度反演位置處反演所得該高度上的水平路徑損耗分布與正演結(jié)果的比較.圖9(見20頁)為88.5 m反演位置處反演所得二維損耗分布與正演結(jié)果的比較.

(a)

(c)

(d)

(e)圖7 不同高度反演位置反演所得90 km處的垂直方向基本傳輸損耗分布與正演結(jié)果的比較

表2為90 km各反演位置處的基本傳輸損耗大小，以及各反演結(jié)果在垂直方向上的損耗分布相對于圖6正演結(jié)果的平均誤差.可以得出：

表2 90 km各位置處反演統(tǒng)計(jì)結(jié)果比較

1) 不同反演位置反演得到的損耗分布在其相應(yīng)位置附近與目標(biāo)值均較接近，由此得出算法在優(yōu)化迭代下已經(jīng)有效收斂.

2) 從損耗分布的大致趨勢上看，由最小損耗(88.5 m處)反演所得結(jié)果與目標(biāo)值一致性最好.即當(dāng)接收點(diǎn)選擇在基本傳輸損耗較小處時(shí)，所獲得的關(guān)于無線電波導(dǎo)分布的信息相對較多，反演效果較好.

圖8 88.5 m高度反演位置處反演所得該高度上的水平路徑損耗分布與正演結(jié)果的比較

3) 由圖8可知，電波傳播距離越遠(yuǎn)，即反演計(jì)算過程中反演位置距發(fā)射天線越遠(yuǎn)，反演計(jì)算結(jié)果的誤差越大.

4) 由圖9可以看出，反演和正演偽彩圖中均存在著明顯的電波陷獲層，分布具有較好的一致性.

綜上所述，驗(yàn)證了該方法的可行性.上述結(jié)論中存在的誤差可以通過增加對反演求解的約束來減小.如通過減小搜索區(qū)間、增加反演損耗信息的輸入、選擇更為合適的粒子群算法參數(shù)來實(shí)現(xiàn).

4 結(jié) 論

本文綜合運(yùn)用PSO、PE等方法，初步探討了如何利用通信鏈路反演無線電波導(dǎo)條件下電波傳輸特性，并通過計(jì)算仿真驗(yàn)證了其可行性.

從反問題信息量的角度看，不同的接收位置處得到的關(guān)于電波空間損耗的信息不同，如何對目標(biāo)實(shí)現(xiàn)更有效的反演仍待探索，如可以通過利用更多接收位置處的功率損耗信息增加對反演的約束.此外，在后續(xù)工作中，還應(yīng)考慮不同波導(dǎo)模型、環(huán)境噪聲等條件對反演的影響.

[1] KARIMIAN A, YARDIM C, WILLIAM S, et al. Estimation of radio refractivity using a multiple angle clutter model[J]. Radio Sci, 2012, 47(3): RS0M07. doi: 10.1029/2011RS004701.

[2] 藺發(fā)軍, 劉成國, 成 思, 等.海上大氣波導(dǎo)的統(tǒng)計(jì)分析[J]. 電波科學(xué)學(xué)報(bào), 2005, 20(1): 64-68.

LIN Fajun, LIU Chengguo, CHENG Si, et al. Statistical analysis of marine atmospheric duct[J]. Chinese Journal of Radio Science, 2005, 20(1): 64-68. (in Chinese)

[3] 成印河, 趙振維, 張玉生. 季風(fēng)期間南海低空大氣波導(dǎo)統(tǒng)計(jì)分析[J]. 電波科學(xué)學(xué)報(bào), 2010, 27(2): 268-273.

CHENG Yinhe, ZHAO Zhenwei, ZHANG Yusheng. Statistical analysis of the lower atmospheric ducts during monsoon period over the South China Sea[J]. Chinese Journal of Radio Science, 2012, 27(2): 268-274. (in Chinese)

[4] SKOLNIK M. Radar Handbook[M]. 2nd ed. New York：McGraw-Hill Publishing Company, 1990.

[5] KARIMIAN A, YARDIM C, GERSTOFT P. Refractivity estimation from sea clutter: an invited review[J]. Radio Sci, 2010, 46(3): RS6013. doi: 10.1029/2011RS004818.

[6]YARDIM C, GERSTOFT P, HODGKISS W S. Tracking Refractivity from clutter using kalman and particle filters[J]. IEEE Trans Antennas Propag, 2008, 56(4): 1058-1070.

[7] GERSTOFT P, ROGERS L T, KROLIK J L, et al. Inversion for refractivity parameters from radar sea clutter[J]. Radio Sci, 2003, 38(3): 8053. doi:10.1029/2002RS002640.

[8] DOUVENOT R, FABBRO1 V, GERSTOFT P, et al. Real time refractivity from clutter using a best fit approach improved with physical information[J]. Radio Sci, 2010, 45(1): RS1007. doi: 10.1029/2009rs004137.

[9] ZHAO Xiaofeng. Evaporation duct height estimation and source location from field measurements at an array of radio receivers[J]. IEEE Trans on Antennas and Propag, 2012, 60(2): 1020-1025.

[10] TRELEA I C. The particle swarm optimization algorithm : convergence analysis and parameter selection[J]. Information Processing Letters, 2003, 85(6): 317-325.