陳勵華 王仕成 劉志國
第二炮兵工程大學,西安 710025
地磁場是地球的固有資源,場強隨空間與時間變化,地磁場導航正是利用這種地磁場在空間分布的各異性來確定載體所在的地理位置。由于地磁定位與地形定位很相似,都是航跡曲線與數字地圖的配準,因此目前地磁場導航技術采用的導航算法主要是借鑒地形高程導航算法中的TERCOM(Terrain Contour Matching)和SITAN(Sandia Inertial Terrain Aided Navigation),即地磁匹配與地磁濾波2種方式[1-6]。與地磁濾波相比,地磁匹配可斷續使用,不存在累積誤差,易于獲得更高的導航定位精度,因此逐漸成為地磁導航的主流方法[7]。
地磁匹配(Magnetic Field Contour Matching,MAGCOM)通過對測量序列與基準圖之間相關性進行度量完成定位。現有的MAGCOM算法在理論上可以只依據磁測設備提供信息單獨實現定位,但實際應用中往往只能作為輔助導航方式,其中基于推位運算的MAGCOM要求航跡為直線,這就最大地限制了地磁匹配的應用。本文重點研究載體行進路徑形狀未知條件下,如何利用相關性度量實現地磁場定位的問題。文中采用了在最優化問題中成功運用的遺傳算法,并對其中交叉與變異概率的選取方法做出相應的改進。仿真實驗證實結合智能優化方法的MAGCOM可以完成地磁匹配導航,相關結論亦可用于地形匹配或重力場匹配定位。
MAGCOM是在預先已知的基準地磁圖中尋找與測量序列度量指標最優的一條對應序列作為對測點位置的估計結果,其實質是數據集的配準問題。該方法可具體描述如下(參見圖1)。

圖1 MAGCOM原理
預先測量地磁場特征量,并繪制成基準圖。一般基準圖都以網格形式表示,其中橫軸和縱軸為地理位置。實際匹配時,根據對載體航行區域的估計,在預存的基準圖中選擇合適的子圖,完成相關性匹配。用于匹配的子圖網格數目設為ML×NL,網格寬度表征了基準圖的分辨率,縱向與橫向分辨率可以相同,也可以不同,圖1中基準圖采用均勻網格,網格大小為d×d,表示不同方向分辨率相同。
匹配時,根據對當地地磁場特征量的實際測量值Tt={Tt(i);i=1,…,N}與基準圖中預存的地磁場數據Tm={Tm(i);i=1,…,N}序列計算相應的相關性度量函數。在所有可能的數據序列中尋找度量指標最優的數據列,其橫、縱坐標即是對真實航跡的估計。常用的度量函數指標包括距離與相似性兩類。前者如平均絕對差(MAD)算法、平均平方差(MSD)算法等;后者如積相關(prod)算法和歸一化的積相關(Nprod)算法等[4-6]。度量中期望距離指標最小或相似性指標最大。
相關性度量僅對地磁測量值完成,因此理論上MAGCOM方法可單獨利用磁測計定位。但是在實際度量中,僅僅已知地磁場測量序列是不夠的。以MAD度量為例說明如下。地磁場MAD度量指標如式(1)所示:
(1)
其中,Tt(i)為第i點實測磁場數值,d為網格寬度,Tm(u,v)為基準圖上存儲的(u,v)位置上地磁數據。由于測量航跡序列Tt與基準數字地磁圖Tm之間對應關系未知,會造成可行解的MAD指標計算數目極其巨大,具體分析如下:在已知的基準圖中搜索指標極值時,路徑起點可能是ML×NL中任意一點,假定相鄰測點間距與基準圖網格寬度d相同,即相鄰點之間存在八鄰域關系,則對于長度為N的測量序列,測量值與基準圖之間所有可能的對應關系共有ML×NL×8N-1種,對所有可行解求取MAD指標,再尋找指標最優解,則會由于運算量過大而導致難以求解。因此現有的MAGCOM方法一般假定測量航跡形狀已知,即各測點之間的相對距離和方向信息都是已知的。在這種假設條件下,可行解的個數為ML×NL,有效地減小了搜索空間。
MAGCOM實現一般需借助慣導系統。慣性導航系統經初始對準后,可提供測點之間距離與方向的信息,地磁匹配運算認為真實航跡和慣導航跡形狀相同(如圖1所示),以減小可行解數目。這樣MAGCOM方法就只能作為輔助導航方式,而無法獨立定位。
采用慣性器件導致導航成本較高,因此有時也會采用里程計,利用推位計算為MAGCOM提供航跡信息[8]。由于里程計只能提供測點之間的距離信息,而無法給出方向信息,所以推位航跡一般設為直線,因此該方法只適用于載體行進路徑為直線的場合,從而限制了地磁匹配的應用。
MAGCOM的實質是在可能的解空間中搜索度量指標最優,遺傳算法是解決搜索問題的一種通用方法,因此本文考慮將遺傳算法應用于地磁匹配定位問題。遺傳算法的引入使得MAGCOM不再關注測點之間的相對方向信息,只需磁測計和里程計即可實現測量,無需慣性器件的參與,減小了匹配成本。與現有基于推位計算的地磁匹配算法相比,載體行進路徑不存在直線約束,可以為任意形狀,因此擴大了匹配方法的應用范圍。
遺傳算法是借鑒生物界中進化與遺傳的機理實現的智能計算,它模擬孟德爾遺傳變異理論和達爾文優勝劣汰思想,通過“遺傳-競爭-再遺傳-再競爭”的方式一步步逼近問題的最優解。遺傳算法求解是從代表問題可能潛在的解集的一個種群開始的,在每一代根據問題域中個體的適應度大小選擇個體,并借助于自然遺傳學的遺傳算子進行組合交叉和變異,產生出代表新的解集的種群。
遺傳算法基本上不用搜索空間的知識或其它輔助信息,而僅用適應度函數值來評估個體,在此基礎上進行遺傳操作。適應度函數不僅不受連續可微的約束,而且其定義域可以任意設定。這一特點使得遺傳算法的應用范圍大大擴展。
以下算法的約束條件為測量序列的各測點之間間距已知,該信息可由里程計測量得到。按照測點間距調整基準地磁圖網格的大小,相應的序列比對中,某點的相鄰測點位置只需在基準圖中該點的八鄰域中搜索。
2.2.1 初始種群
基于遺傳算法的MAGCOM中,基準圖上每一條長度為N的序列都對應一個可行解,即種群中的一個染色體,可表示為有序點連接而成的折線:
path=[x1,y1;x2,y2;...;xN,yN]
(2)
其中xi,yi表示路徑上第i點在基準圖上的坐標,設xi,yi取值為1~ML和1~NL范圍的實數,且相鄰點滿足八鄰域關系:
(xi-xi-1,yi-yi-1)∈{(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)}
根據相鄰點的鄰域關系,染色體亦可用相對坐標表示為
path1=[0,0;x2-x1,y2-y1;...;
(3)
相鄰點的鄰域關系僅有8種,因此path1的表達形式比path簡單,但缺少起始點的信息。
初始解path的產生可以通過以下3個步驟完成:
1) 序列的起始點(x1,y1)產生。起始點可能在已知基準圖中任何一點,因此可以隨機的產生取值為1~ML的任意實數x1和1~NL的任意實數y1;
2) 相對坐標序列path1的產生。隨機產生長度為N-1的相鄰點八鄰域關系,并按式(3)的規則產生path1;
3) 根據式(2)和(3)關系可計算得到初始解path。
按照以上步驟可以產生m個初始解,形成初始種群。
2.2.2 適應度函數
其中, bij(0-t)表示在時間段 (0-t) 內, i地區j產業的增長速度,表示在時間段(0-t)內j產業在全國的增長速度。
解的優劣可以通過該路徑與測量值的相似性程度來體現。適應度函數選為相似性度量指標如MAD,NPROD等。MAD指標對應的適應度函數如下:
NPROD指標對應的適應度函數如下:
適應度函數根據某一染色體path中每一個點對應的基準圖中地磁場強度大小Tm與實際測量序列中對應的磁測值Tt進行比較、計算得到。對某染色體而言,以上2個適應度指標越大,則意味著該染色體適應能力越強。
2.2.3 遺傳算子
遺傳算法的選擇算子采用輪盤賭的方式,適應度越高的個體被選中的概率越大。在交叉與變異中為保證運算后能滿足相鄰點的八鄰域要求,執行運算時區分起點和其余各點。起點的交叉和變異對path中的第1組坐標(x1,y1)進行,其余各點則利用path1表達式,即改變點之間的鄰域關系。經過交叉與變異后,路徑上點間的鄰域關系發生變化,獲得新個體。
遺傳算法中交叉概率Pc和變異概率Pm的選擇是影響算法性能的關鍵。交叉概率Pc過小會導致搜索過程變慢,過大會破壞優秀個體的結構;變異用于表現基因的突變,可使求解過程跳出局部搜索空間,變異概率Pm過小不易產生新的全局極值,過大則遺傳算法會成為純粹的隨機搜索。因此Pc和Pm需要反復試驗確定。Srinvivas等人提出根據適應度大小自適應調整Pc和Pm的方法[9],如式(6)和(7)所示:
(6)
(7)
式中,fmax為群體中最大的適應度,favg為每代群體平均適應度,f′為交叉的兩個體中較大適應度,f為變異個體的適應度。
本文研究的長度為N的有序點集配準問題,實際是個多級尋優問題。結合多級尋優具有可加性特點,提出以局部適應度函數對以上自適應調整Pc和Pm的方法作進一步改進。最優解從起點開始,以小于等于N的任何一級作為終點,得到的序列仍應是該長度下的優解。為保護優秀基因不在交叉和變異中遭到破壞,設置局部適應度函數。以MAD指標為例,第n級局部適應度函數如式(8)所示,其中n為[1,N]中的任意整數。
(8)
與式(4)相同,式(8)的分母部分表示的也是被搜索序列與測量序列相比平均每點的距離偏差,因此可以將局部適應度與整條序列的適應度比較。
交叉可單點進行,也可多點進行。本文采用單點交叉的方式,隨機選擇交叉點,依據式(8)計算擬交叉的2條序列在該點的局部適應度d1和d2,令f′=max(d1,d2),由式(6)計算可得Pc,同理以局部適應度取代序列適應度計算Pm。
局部適應度的作用是針對那些包含優秀基因的劣解,能在減小優秀基因段被改變概率的同時盡量增大不良基因交叉和變異的可能。
2.2.4 迭代終止條件
遺傳算法循環迭代搜索,迭代次數是影響計算量的重要因素。只要滿足以下任一條件,即停止迭代:1)設定最大遺傳代數Gn,當迭代次數到達Gn,即停止迭代,以當前最佳個體作為匹配結果; 2)迭代到一定程度可能已無法得到適應度更優的結果,若連續g代最優適應度未有改進,即停止迭代。
2.2.5 算法的復雜度分析
算法效率可用時間復雜度分析。以比較運算作為基本操作,在本文的約束條件下,現有MAGCOM方法時間復雜度為o(ML×NL×8N-1),上述遺傳算法中,種群大小m序列長度為N的染色體,最大迭代次數為Gn,整個計算過程的時間復雜度為o(m.N.Gn)。當序列長度較大時,遺傳算法的計算量優勢明顯。
為了驗證本文提出算法的有效性,在實測地磁圖上進行仿真試驗。該圖來自我國西部某地區實地測量結果。地磁圖網格數為30×30,網格間距為20m×20m,對應測點間距亦為20m。仿真中適應度函數取MAD 指標,序列長度取N=10,交叉概率Pc與變異概率Pm如前所述自適應獲得,初始種群大小m=50,遺傳代數Gn=100。
在不考慮噪聲的條件下,對某測量序列在整個基準圖中搜索,定位結果如圖2所示。圖中序列的10個測點中有8個完全匹配,整個序列的總誤差為2.828網格,且誤差隨種群規模和迭代次數的增大、參數設置的改進還有望減小。這表示遺傳算法能夠在已知測點間距的條件下,利用相似性比較獲得測量路徑的形狀,因此不必依賴慣性系統信息,即可完成定位。

圖2 仿真結果
遺傳算法中解的生成具有不確定性,而且地磁匹配定位還需考慮測量噪聲與基準圖制備等誤差因素影響,因此很難保證定位結果足夠理想。本節從基準圖大小、序列長度、噪聲條件等方面對地磁匹配進行仿真,尋找以上因素對算法影響的規律以及改進措施。以下仿真中,如未作特殊說明,仿真條件同2.3節。
在不考慮噪聲的條件下,選取序列長度為N=10,在不同基準圖大小條件下完成匹配。基準圖網格數分別為10×10,20×20,30×30,經過500次匹配,計算所得的平均誤差如表1所示。表中誤差亦指整個序列誤差。

表1 基準圖網格數的影響
由于遺傳算法參數選取一致,因此表1中3種情況搜索的運算量是相同的。然而基準圖越大,會導致匹配空間越大,因此遺傳算法越容易落入指標的局部極小值點而無法實現序列的準確比對,造成誤差隨基準圖增大而增大。對于本例而言,如表1所示,即使在基準圖為30×30的仿真條件下,整個序列的定位誤差仍可達到5個網格左右,相應的平均單點誤差僅有0.5個網格,定位精度比較理想。

綜合分析序列長度與噪聲的影響。在N=5,N=10兩種序列長度下,匹配基準圖選為20×20,分別疊加如上所述高斯噪聲和均勻噪聲,其它參數設置不變。對各種情況測試500次,平均結果如表2所示。

表2 序列長度與噪聲的影響
將表1和表2中序列長度N=10,基準圖20×20的結果對比可以發現,均勻噪聲的影響非常小,有無噪聲條件下誤差差別不大;高斯噪聲則影響明顯,會引起匹配精度的顯著下降。另一方面,序列長度會引起整條序列誤差的增加,但平均每個點所引起的誤差隨序列長度增長并無明顯增大的表現。
本文提出利用遺傳算法實現的地磁匹配。該方法有效解決解空間過大的難題,計算量較小。具體實現中以相似性度量作為適應度函數,采用相對路徑完成交叉和變異,并對交叉變異的概率自適應調整方法實現了改進。通過仿真討論了基準圖大小、噪聲、匹配序列長度對算法的影響,發現基準圖增大、高斯噪聲會引起匹配誤差的明顯增加,而序列長度影響不明顯。結果表明,該方法降低了現有MAGCOM方法在使用中的約束條件,在僅提供測點地磁場信息和相對距離信息的條件下,能依靠算法獲得路徑形狀,以較小誤差實現定位。
參 考 文 獻
[1] Mu Hua, Wu Mei-ping. Geomagnetic Surface Navigation Using Adaptive EKF[C].2007 IEEE Conference on Industrial Electronics and Application,Harbin,2007.
[2] E R Benson,T S Stombaugh,N Noguchi, et al. An Evaluation of a Geomagnetic Direction Sensor for Vehicle Guidance in Precision Agriculture Applications[C]. ASAE Paper 983203, ASAE.St. Joseph,MI,1998.
[3] Mohammad Nizam Filipski,Ermira Junita Abduliah. Nanosatellite Navigation with the WMM2005 Geomagnetic Field Model[J].Turkish Journal of Environmental Sciences, 2006,30(1): 43-55.
[4] 劉穎,吳美平.地磁匹配中的匹配長度估計方法[J].中國慣性技術學報,2010,18(4):439-443.(Liu Ying, Wu Mei-ping.Match Length Estimation in Geomagnetic Matching System [J]. Journal of Chinese Inertial Technology, 2010,18(4):439-443.)
[5] 劉飛,周賢高,等.相關地磁匹配定位技術[J].中國慣性技術學報,2007,15(1):59-62.(Liu Fei, Zhou Xian-gao,et al. Geomagnetic Matching Location Using Correlation Method[J]. Journal of Chinese Inertial Technology, 2007,15(1):59-62.)
[6] 任治新,楊功流,李士心.基于向量搜索的地磁導航匹配算法[J].中國慣性技術學報,2008,16(4):424-427.(Ren Zhixin,Yang Gongliu,Li Shixin.Geomagnetic Matching Algorithm Based on Vector Search[J]. Journal of Chinese Inertial Technology, 2008,16(4):424-427.)
[7] 周軍,葛致磊,施桂國,等.地磁導航發展與關鍵技術[J].宇航學報,2008,29(5):1467-1742.(Zhou Jun, Ge Zhilei, Shi Guiguo, et al. Key Technique and Development for Geomagnetic Navigation[J]. Journal of Astronautics, 2008, 29(5):1467-1742.)
[8] 喬玉坤.地磁導航基準圖制備與數學仿真方法研究[D].西安:第二炮兵工程學院,2011.(Qiao Yukun.Study on the Preparation and Numerical Simulation of Reference Maps for Geomagnetic Navigation[D]. Doctor Dissertation of Second Artillery Engineering Institute,Xi-an,2011.)
[9] Srinivas M,Patnaik L M.Adaptive Probabilities of Crossover and Mutations in Gas[J]. IEEE Transactions on SMC,1994,24(4):656-667.