粒子群優(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ā)生概率達到了70%[3])，這就為無線電波導(dǎo)條件下實現(xiàn)電磁波的超視距傳播提供了有利條件.

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

1 研究方法

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

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

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

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

無線電波導(dǎo)通常被特征化為依賴于高度和距離的環(huán)境折射率.垂直方向上，表面波導(dǎo)修正折射率剖面一般通過三線模型近似，而蒸發(fā)波導(dǎo)可以通過指數(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個主要參數(shù)如下：Md為逆變層逆差；zthick為表面波導(dǎo)的厚度；m1為混合層斜率；h1逆變層底高；hd為蒸發(fā)波導(dǎo)高度.五參數(shù)模型如圖2所示.

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

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

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

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

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

(2)

此處5個水平參數(shù)的變化范圍應(yīng)小于對應(yīng)特征值的方根，這樣才能與大氣折射率的微小變化相符[7].至此，通過10個參數(shù)(水平方向上的5個參數(shù)和垂直方向上的5個參數(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)垂直剖面模型即可計算電磁波基本傳輸損耗.

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

F[u(x,z)]}，

(4)

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

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

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)初始模型的選取遠離真實模型時將很容易陷入局部極小值中得出錯誤的求解；同時該方法還需要計算偏導(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ù)難于控制、計算工作量大、收斂速度慢、計算時間長，在實際反演應(yīng)用中受到了限制[9-10].最近發(fā)展的貝葉斯濾波的方法理論相對復(fù)雜，實現(xiàn)較難.

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

(7)

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

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

fitness=|Lbo-Lbr(x)|.

(8)

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

圖5 粒子群算法流程

3 仿真分析

3.1 仿真條件和主要步驟

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

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

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

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

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

步驟2反演計算

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

3.2 反演結(jié)果分析

在相同無線電波導(dǎo)和發(fā)射天線情況下，分別利用傳播距離上90 km處的幾個不同高度反演位置的基本傳輸損耗進行該方法的反演計算，其中，最小損耗位于88.5 m高度處.每種情況均進行20次反演計算，并對反演結(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)計結(jié)果比較

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

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

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

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

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

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

4 結(jié) 論

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

從反問題信息量的角度看，不同的接收位置處得到的關(guān)于電波空間損耗的信息不同，如何對目標(biāo)實現(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]. 電波科學(xué)學(xué)報, 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]. 電波科學(xué)學(xué)報, 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.