衛星載荷磁場優化計算方法研究

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

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

0 引言

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

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

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

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

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

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

Hm(r)=

(1)

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

在不存在自由電流的區域中，從磁荷觀點出發，均勻磁化衛星載荷產生的標量磁位[8]可表示為

(2)

(3)

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

(4)

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

(5)

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

(6)

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

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

(7)

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

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

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

2.1 微粒群算法的基本原理

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

(8)

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

2.2 計算方法與流程

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

(9)

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

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

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

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

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

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

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

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

3 實驗方案與結果分析

3.1 實驗方案

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

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

3.2 實例計算過程說明

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

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

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

4 結束語

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

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

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

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

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

[4] 林瑞淋, 張效信, 劉四清. 高緯磁層頂位形統計分析[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]. 北京: 機械工業出版社, 1988: 65-68.

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

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