李 恪，姚崇斌，徐紅新，謝寶蓉

(上海航天電子技術(shù)研究所，上海 201109)

0 引言

地磁測量是衛(wèi)星弱磁探測技術(shù)的重要應用之一，對整個地球的磁場環(huán)境描述有著重要的意義[1-4]。因衛(wèi)星與地球相距非常遠，故衛(wèi)星所處空間的地磁感應強度很小。但衛(wèi)星在軌工作時，整星內(nèi)部會形成一個基本穩(wěn)定的磁場。衛(wèi)星內(nèi)部磁場與地磁場相互作用形成磁場干擾力矩，嚴重影響姿態(tài)控制及地磁測量的精度。為減少星上磁場噪聲的影響，需要在地面上測量載荷的磁場，并做消磁處理。

現(xiàn)有的磁場測量方法大多采用數(shù)值計算方法，將被測物體劃分為許多小的磁場單元，對各個單元在不同方向上的磁場分量進行積分，從而得到被測物體總的磁場，如有限元法、積分方程法等。有限元法要求在整個求解區(qū)域內(nèi)進行離散剖分，計算數(shù)據(jù)的前處理較繁雜，且計算開域問題時需要截斷邊界，故有限元法的使用范圍受限。積分方程法中系數(shù)矩陣元素的計算為三重體積分，計算較復雜，且當源點和場點重合時，系數(shù)矩陣的計算帶有奇異積分。該方法的優(yōu)點是剖分簡單，待求解未知量少，但非對稱滿陣導致占用資源大，計算時間長。

地磁場是一個弱磁場，對于處于不同地磁場中的載荷，可認為其磁化率基本不變。故針對數(shù)值計算方法的不足，根據(jù)實際應用的需求，采用三維靜磁場逆計算方法的思路[5-7]，對傳統(tǒng)的數(shù)值計算方法加以優(yōu)化。

本算法利用某處地磁場下衛(wèi)星載荷的磁場測量數(shù)據(jù)，逆計算衛(wèi)星載荷的等效磁化率，根據(jù)等效磁化率計算衛(wèi)星載荷產(chǎn)生的感應磁場，根據(jù)均勻磁化衛(wèi)星載荷的特點，將積分方程法中系數(shù)矩陣元素的體積分形式化為面積分形式，并消除積分計算的奇異性，建立基于單元表面積分的感應磁場計算模型。最后，以實際載荷為實驗對象，驗證算法的有效性。相比傳統(tǒng)算法，在保證精度的前提下，極大地提高了工作效率。

1 基于表面積分的均勻磁化體磁場計算模型

基于體積分的積分方程法磁場計算公式為[2]

Hm(r)=

(1)

式中：Hm為鐵磁物體產(chǎn)生的磁場；M為鐵磁物體內(nèi)部磁化強度；r為場點的坐標；r′為源點的坐標；V為鐵磁物體的體積。由式(1)計算鐵磁物體的磁場時需要計算三重積分，計算較復雜，且當場點和源點重合時，積分計算具有奇異性。

在不存在自由電流的區(qū)域中，從磁荷觀點出發(fā)，均勻磁化衛(wèi)星載荷產(chǎn)生的標量磁位[8]可表示為

(2)

(3)

式中：dS為單位面元的投影面積；n為曲面S在r處的法向矢量。故均勻磁化衛(wèi)星載荷在空間產(chǎn)生的磁場強度可表示為

(4)

一般情況下，不規(guī)則衛(wèi)星載荷在外磁場中未被均勻磁化。若將衛(wèi)星載荷剖分成N個小單元，當剖分單元體積足夠小時，則可認為每個小單元是等效均勻磁化體。由磁場疊加原理可知，衛(wèi)星載荷在空間產(chǎn)生的磁場強度為所有等效均勻磁化體的疊加，即式(4)可表示為

(5)

考慮到每個剖分單元衛(wèi)星載荷的線性度和均勻性，設每個剖分單元滿足Mi=χiHi， 式(5)可進一步表示為

(6)

式中:χ=[χ1,χ2,…,χN]為待優(yōu)化的參數(shù)組。一般情況下，該問題很難直接求解，但通過迭代過程，即多次計算不同的χ所對應的目標函數(shù)值，即可逼近問題的真實解。當已知每個剖分單元的磁化率和磁場強度時，根據(jù)式(6)就可計算出衛(wèi)星載荷在空間任一點產(chǎn)生的感應磁場。

假設衛(wèi)星載荷處于外磁場為He的空間，則空間任一點的磁場H=He+Hm。當把場點配置在第j個剖分單元的中心時，可得到第j個剖分單元中心磁場強度的表達式為

(7)

寫成矩陣形式為AH=He，其中系數(shù)矩陣A由場點坐標、剖分單元中心坐標及單元磁化率決定，H為每個剖分單元中心的磁場強度。

已知衛(wèi)星載荷的磁性參數(shù)，根據(jù)式(7)建立每個單元中心處磁場強度的線性方程組，求解該方程組即可得到每個剖分單元的磁場強度，再由式(6)就可計算衛(wèi)星載荷在空間任一點產(chǎn)生的感應磁場。由此可看出，系數(shù)矩陣元素的計算為面積分的形式，較傳統(tǒng)的體積分形式簡單；由于積分是在剖分單元表面進行，而場點選取在剖分單元的中心，故當源點和場點重合時，積分計算無奇異性。

2 基于微粒群算法的等效磁化率逆計算方法

2.1 微粒群算法的基本原理

微粒群算法中，所有粒子都有一個由被優(yōu)化函數(shù)決定的適應值，每個粒子還有一個速度用于決定他們搜索的方向和距離，粒子們追隨當前的最優(yōu)粒子在解空間中搜索。在求解過程中，微粒群初始化為一群隨機粒子(隨機解)，并通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自身，第一個極值是粒子本身所找到的最優(yōu)解，即個體極值；另一個極值是整個種群目前找到的最優(yōu)解，即全局極值。找到這兩個最優(yōu)值后,粒子根據(jù)如下的公式來更新自己的速度和新的位置[9]

(8)

式中：vi為粒子的速度；xi為當前粒子的位置；pbest為個體極值；gbest為全局極值；xrand為介于0與1之間的隨機數(shù)；c1、c2為學習因子；w為慣性權(quán)重；t為迭代的次數(shù)。每一維粒子的速度都會受限于最大速度vmax，如果粒子某一維更新后的速度超過vmax，那么這一維的速度就被限定為vmax。

2.2 計算方法與流程

由于外加磁場的影響，衛(wèi)星載荷的磁性狀態(tài)會根據(jù)外磁場的大小，沿著其磁化曲線工作于相應的磁化點。電磁場逆問題是由果求因的問題[10]，在靜磁場計算中即為已知磁場測量數(shù)據(jù)來反演衛(wèi)星載荷的磁性參數(shù)。建立的逆問題數(shù)學模型可表示為

(9)

由式(9)可知，模型中有N個未知數(shù)，故將粒子的解空間設置成N維空間，每個粒子的坐標為[χ1,χ2,…,χN]，粒子的適應度函數(shù)取為式(9)中的J，每個粒子具有的速度v=[v1,v2,…,vN]，第t+1代粒子在第t代粒子的基礎上按式(8)更新自己的位置和速度，其優(yōu)化過程如下：

1) 設置相關初始參數(shù)，如粒子群規(guī)模m、慣性權(quán)重w、學習因子c1和c2、粒子最大飛行速度vmax、迭代結(jié)束條件；

2) 初始化粒子的位置及速度；

3) 評價各粒子的適應度函數(shù)；

4) 對每個微粒，將其適應度值與歷史最好位置pbest相比較，若當前適應度值更優(yōu)，則用當前適應度值更新pbest；

5) 將每個微粒的適應度值與群體經(jīng)歷過的歷史最佳位置gbest相比較，若當前群體中最好的適應度值較好,則將其設置為新的gbest；

6) 根據(jù)式(8)更新各微粒的位置和速度；

7) 計算優(yōu)化目標函數(shù)，若達到精度要求或迭代次數(shù)，則返回當前最佳微粒的結(jié)果，算法結(jié)束，否則返回步驟3)。

3 實驗方案與結(jié)果分析

3.1 實驗方案

圖1 實驗方案及現(xiàn)場示意圖Fig.1 Schematic diagrams of experimental scheme and experimental site

為考察不同場點、不同外磁場及不同磁化方向下計算的感應磁場的準確性，測量時，被測物體先沿南北方向滑動測量，再轉(zhuǎn)動測量架，使被測物體沿東西方向滑動測量。地磁平面的選取標準原則上是離被測物體越近越好，但是要考慮到其他設備的干擾，因此需要綜合考慮取測量偏差較小的距離。用上述相同的方法測量縱向載荷上端高度為h2、載荷下端高度為h3，以及橫向載荷下端高度為h4時的平面磁場數(shù)據(jù)，以用于與計算值比較。

3.2 實例計算過程說明

應用本文算法計算衛(wèi)星載荷感應磁場過程中，兼顧計算時間和計算精度，將載荷剖分為20×5個長方體單元，并假定每個剖分單元為等效均勻磁化體。應用微粒群優(yōu)化算法反演計算磁化率的過程中，考慮到載荷材料的均勻性，認為每個剖分單元在地磁場的作用下的磁化率相同，故式(8)可簡化為只含一個未知數(shù)的優(yōu)化問題，相應的粒子的解空間為1維空間，即待優(yōu)化的參變量χ=[χ,χ,…,χ]。

用微粒群優(yōu)化算法進行優(yōu)化求解時，參數(shù)設定為：微粒粒子數(shù)m=40，慣性權(quán)重w=0.9，學習因子c1=c2=2，最大速度vmax=10，迭代次數(shù)50。為考察微粒群算法反演磁化率的有效性， 逆計算15次。圖2為15次逆計算所得的等效磁化率和最終最優(yōu)粒子的適應度曲線和一次逆計算時粒子適應度的收斂曲線。由圖2可知，用微粒群優(yōu)化算法反演磁化率是有效的，且收斂速度較快。計算衛(wèi)星載荷感應磁場時所用磁化率取15次逆計算得到的等效磁化率的平均值為140.8。圖3為衛(wèi)星載荷受地磁場水平分量縱向磁化和橫向磁化時，產(chǎn)生的感應磁場三分量測量值和計算值對比曲線。由圖3可看出，測量值和計算值吻合較好，其最大相對誤差約為10.79%。

圖2 磁化率逆計算參數(shù)圖和收斂曲線Fig.2 Parameters and convergence curve of magnetic susceptibility inverse calculation

4 結(jié)束語

由均勻磁化體的磁化特點出發(fā)，建立了基于單元表面積分的衛(wèi)星載荷感應磁場求解方法。與傳統(tǒng)基于體積分的磁場求解方法相比，該方法簡化了系數(shù)矩陣運算，消除了源點和場點重合時的積分奇異性。以衛(wèi)星載荷的磁化率為未知量，建立了求解磁化率的數(shù)學模型，并運用隨機類微粒群優(yōu)化算法對其求解。最后，以衛(wèi)星載荷為計算實例，驗證了本文算法的可行性，計算值與測量值具有較好的一致性。結(jié)果表明，該方法完全可用于衛(wèi)星載荷感應磁場的計算。

圖3 不同磁化方向、不同場點的感應磁場Fig.3 Induction magnetic field at different magnetization directions and field points

[1] 吳衛(wèi)權(quán). 地磁場中衛(wèi)星不側(cè)置態(tài)下磁矩測試方法[J]. 上海航天, 2014, 31(1): 50-55.

[2] 徐文耀. 地磁場的三維巡測和綜合建模[J]. 地球物理學進展, 2007, 22(4): 1035-1039.

[3] 傅綏燕， 宗秋剛, 濮祖蔭. 地磁活動與太陽活動對環(huán)電流離子成分與位置的影響[J]. 地球物理學報, 2003, 46(6): 725-730.

[4] 林瑞淋, 張效信, 劉四清. 高緯磁層頂位形統(tǒng)計分析[J]. 地球物理學報, 2003, 53(1): 1-9.

[5] PECHENKOV A N . A mathematical algorithm for solving the inverse problem of magnetostatic flaw detection[J]. Russian Journal of Nondestructive Testing, 2005, 41(11): 714-718.

[6] PECHENKOV A N. Numerical simulation of an inverse problem of magnetostatic flaw detection[J]. Russian Journal of Nondestructive Testing, 2005, 41(11): 719-723.

[7] 周克定. 工程電磁場專論[M]. 武漢: 華中工學院出版社, 1986: 28-29.

[8] 樊明武, 顏威利. 電磁場積分方程法[M]. 北京: 機械工業(yè)出版社, 1988: 65-68.

[9] 連麗婷, 楊明明, 孫開江. 基于改進微粒群算法的磁靶電流參數(shù)優(yōu)化調(diào)整[J]. 艦船科學技術(shù), 2015, 37(3): 73-76.

[10] 黃卡瑪, 趙翔. 電磁場中的逆問題及應用[M]. 北京: 科學出版社, 2005: 13-19.