基于四叉樹的海浪磁場快速仿真算法

熊雄，楊日杰，韓建輝，郭新奇

(1.海軍航空工程學院 電子信息工程系，山東煙臺264001;2.海軍航空工程學院指揮系，山東 煙臺264001)

在潛艇聲隱身性能越來越好的條件下，航空磁探儀由于幾乎不受傳播介質特性的影響，成為一種有效的反潛探測設備，在航空反潛中得到廣泛應用。雖然航空磁探潛受海洋環境復雜傳播介質的影響較小，但是隨著航空磁探儀靈敏度越來高，潛艇磁異常信號也越來越容易受到各種背景擾動的干擾。在反潛機飛行高度較低的情況下，海洋中的風浪等海水運動切割地球磁場，激發感應電磁場，對航空磁探儀的工作性能產生重要的影響。大量的實驗表明海浪產生的感應電磁場噪聲與要檢測的目標磁場量級、頻帶基本接近，這些噪聲都是不可忽略的干擾源［1-2］。然而海浪磁噪聲很難直接測量，基于線性波浪理論對海浪磁場進行建模和仿真是研究海浪感應磁場噪聲特點及能量分布特性等因素的一種重要方法［3］。

Weaver在海浪波高為常數的假設下，給出了理想海域條件下單頻重力波產生的感應磁場的理論表達式，并形成了經典 Weaver海浪磁場模型［4］。Ochadlick對Weaver海浪磁場模型進行驗證，驗證結果表明該模型具有較強的可信性［5］。近年來，有學者分別基于Weaver海浪磁場模型給出了有限深海域磁場模型和海浪磁場矢量理論模型以及計算方法［6-7］。但是以上模型和計算方法都是基于單頻海浪重力波，而實際中海浪是一種復雜的海水運動，其浪高隨頻率變化，實際應用中描述海浪特征的方法是海浪譜分析，完整的海浪譜由頻率譜和方向函數組成［8］。文獻［9］基于海浪譜推導了航空磁探儀接收到的海浪磁噪聲功率譜的理論表達式，文獻［10-11］給出了基于海浪頻率譜等分法的海浪磁場數值仿真方法，但是該方法仿真速度慢，而且沒有考慮方向函數的影響。文獻［8］在頻率等分法的基礎上采用傅里葉反變換來模擬海浪，該方法可以減少仿真次數，提高仿真速度，但是該方法無法應用到海浪磁場仿真。

針對以上問題，本文基于Weaver單頻波海浪磁場模型，建立地理坐標系下任意方向傳播實際海浪磁場數值模型，給出快速數值仿真算法，實現海浪磁場的快速仿真，并對仿真結果進行校驗和分析。

1 單頻重力波磁場模型

1.1 Weaver海浪磁場基本理論

根據Weaver海浪磁場理論，在地磁場中海洋重力波感應電場E，磁感應強度B滿足麥克斯韋電磁理論:

式中:μ為海水磁導率，ε為海水介電常數。海水電流傳導密度為J=σ(E+V×BE)，σ為海水電導率，BE表示地磁場矢量，V為海洋重力波速度矢量。

通常認為海水傳導電流密度遠大于式(2)中等號右邊第2項的位移電流，因此可以忽略位移電流。則海浪感應磁場B可以表示為［11］

為了求解B，Weaver引入速度勢φ，且定義V=?φ，并求解得到沿X軸方向傳播單頻波海浪磁場強度。本文將建立地理坐標系下任意方向傳播海浪磁場強度數學模型。

1.2 單頻重力波感應磁場的數學模型

建立地理直角坐標系OXYZ，OXY位于平均海平面，OZ垂直向上。OW為海浪傳播方向，OW與OX軸的夾角為θ，航空磁探儀在海平面上方高度Zm沿著OM直線飛行，飛行路徑OM與OX軸的夾角為β，如圖1所示。Z＞0為空氣介質，Z＜0為海水介質。ON為磁北方向，地磁場矢量BE表示為BE=，I表示磁傾角，γ表示磁北方向與X軸的夾角，如圖2所示。

圖1 地理坐標系Fig.1 Geographic coordinate system

圖2 地磁場矢量示意圖Fig.2 Schematic of the geomagnetic vector

假設海水是不可壓縮無旋流體，根據文獻［10］，波浪沿θ方向傳播，以簡諧運動描述單頻波海浪流體運動，則流體擾動速度勢可以表示為

式中:Ω=xcos θ+ysin θ，a、ω 、k分別表示單頻波幅度、頻率、波數，g為重力加速度，k和ω的散布關系可以表示為

將式(4)代入式(3)求解得到坐標點(x，y，z)處t時刻單頻波海浪磁場B(x，y，z，t)為

式中:hB(z，θ)為磁場幅度矢量。

海浪磁場傳播經過海水和空氣兩層介質。根據邊界條件z=0處海浪感應磁場的垂直分量連續，并且結合式(7)可以得到磁場幅度標量

其中:

則t時刻海平面上方坐標點(x，y，z)處標量磁探儀探測到單頻重力波標量海浪磁場為

式中:ε為海浪初始相位，在(0，2π)上均勻分布。

2 基于線性理論的海浪磁場模型

海浪是一種復雜的隨機運動過程，在海洋學研究中利用隨機過程來描述海浪是進行海浪研究的主要途徑之一。為了模擬實際的海洋環境，根據Longuet-higgins線性波浪理論，t時刻位于坐標點(x，y)處波面 ζ(x，y，t)可以表示為無限個隨機相位正弦波的疊加:

式中:εn為第n個組成波的初始相位，an、ωn、kn第n個組成波的幅度、頻率、波數，θ為波浪主傳播方向。

根據線性理論，海浪產生的磁場也可以表示為無限個單頻重力波海浪感應磁場的疊加。由式(9)、(10)可以得到磁探儀靜止條件下探測海浪磁場為

其中:

在航空磁探測過程中，磁探儀是隨著飛機運動的，因此其接收到的海浪磁噪聲不僅隨時間變化而且隨觀測位置變化。根據式(11)可以得到以速度v按照航向β飛行的航空磁探儀探測海浪磁場為

其中:

3 海浪磁場快速數值仿真算法

3.1 基本仿真算法

式(11)、(12)給出了實際海浪磁場的數值計算模型，但是模擬實際的海浪磁場需要根據數值模型給出有效的數值仿真算法。由式(11)、(12)可知，海浪磁場與波浪頻率以及相應頻率的波高有關，而且海浪具有三維不規則性，不僅波高不同、頻率不同，而且會從各個方向傳到某一點，除沿主風向產生的主浪以外，在主浪向兩側±π/2角度范圍內都有諧波的擴散。描述海浪三維不規則性常用的方法是海浪譜［12-13］。

海浪譜定義為單位頻率間隔和方向間隔內的海浪平均能量密度，二維海浪譜又叫方向譜，方向譜S(ω，θ)是由頻率和角度相關的2個函數組成，可表示為

式中:S(ω)為海浪頻譜，G(ω，θ)為海浪方向分布函數。

海浪頻譜比較容易觀測，國外根據大量海浪觀測資料，提出了許多的海浪頻譜模型。Person-Moscowitz譜模型簡稱P-M譜，能較好地描述風速為0～20 m/s之間的海浪譜。P-M譜模型表示為

式中:S(ω)為能量頻譜，ω為頻率，α=0.008 1，β=0.74，g為重力加速度，U為海面上19.5 m處風速，譜峰頻率為ωn=8.565/U。

根據ITTC的觀測資料方向分布函數可以簡化表示為

根據方向譜可以得到在 ωi-Δωi/2～ωi+Δωi/2頻段和θi-Δθi/2 ～ θi+Δθi/2角度內海浪波高ai，j，可以表示為［12-13］

海浪磁場仿真算法具體步驟如下:

1)海浪頻段的選擇。為了提高仿真速度和仿真時間，需要對海浪頻段進行估計。設定風速，根據海浪譜表達式(13)對海浪譜的頻段進行估計，選擇有限的頻段ω1～ωn來數值計算。

2)進行頻段和方向的離散化采樣。根據海浪譜密度函數，對頻譜和方向進行離散化采樣，頻率采樣間隔為Δω，方向的采樣間隔為Δθ。

3)計算每個離散網格上海浪波高。根據式(16)可以得到ωi和θj對應網格下的海浪的波高ai，j。

4)產生隨機相位εn。利用隨機數生成原理產生[0，2π)之間均勻分布的隨機數。

5)單頻波海浪磁場的計算。設定坐標位置點(x，y，z)，根據式(9)計算 ωi和 θj對應網格下的單頻重力波產生的磁場信號模值。

6)海浪磁場信號的合成。根據式(11)或式(12)進行單頻波海浪磁場信號的線性疊加，計算磁探儀靜止或者運動條件下接收到的磁噪聲信號。

3.2 基于四叉樹分解仿真算法優化

一種簡單的采樣方法就是區間等分法:將頻率區間和方向區間分別進行M、N等分，取固定大小的采樣子區域Δω×Δθ，使，將每個采樣子區域中心對應的頻率和方向角作為單元波的頻率和方向，按照3.1節仿真算法將不同振幅和頻率的單頻波合成得到海浪磁場。區間等分法算法簡單、容易實現，仿真過程中為了盡可能的精確，采樣數相對要大，時空消耗大，不適合在線計算。

四叉樹分解合成算法具有節省存儲空間，提高運算速度等優點，適合于快速計算，廣泛應用于地形學圖形繪制和圖像處理［14-17］。相對于其他多叉樹算法，四叉樹具有結構簡單，檢索效率高的優點。四叉樹分解的基本思想是將二維平面按4個象限進行遞歸分割，直到子象限的數值符合設定的條件，從而得到一棵四分叉的倒向樹。四叉樹分解的示意圖如圖3所示。在四叉樹分解中，一個根結點有4個子結點，這4個子節點按順序標為東北(NE)、西北(NW)、西南(SW)、東南(SE)4個子區域，4個子區域將原圖形區域四等分。依此判斷4個子區域是否滿足進一步分解的條件，如果不滿足分解條件則，子圖形成為葉子節點并存儲該節點，如果滿足分解條件，則子圖形成為根節點進一步分解為4個節點，依此遞歸循環直至分解結束。

按照四叉樹分解的思想，提出基于四叉樹分解的海浪磁場快速優化仿真算法。其優化方法是在3.1節海浪磁場數值仿真算法步驟2)中采用等能量四叉樹分解方法，其步驟2)可以分解為以下步驟:

1)設定每個網格最小采樣能量與總能量的比例PE。

2)根據步驟1)中選擇的頻段和角度范圍，將S(ω，θ)進行四叉樹遞歸分解。

3)判斷每個子節點是否滿足該網格的能量小于或者等于設定比例，若不滿足條件則繼續分解該網格，若滿足條件，則結束分解并記錄葉子節點網格信息。

4)將每個網格葉子節點按照樹形鏈表結構記錄，在后續的仿真過程中采用樹形鏈表遍歷的方法快速遍歷每個葉子節點，得到步驟3)中所需的每個分解單元的信息。

圖3 四叉樹分解示意圖Fig.3 Schematic of quadtree division

4 海浪磁場數值仿真及分析

基于海浪磁場快速仿真算法，對不同條件下海浪磁場進行數值仿真計算，并分析仿真性能?；镜姆抡鏃l件為:地磁場BE=50 000 nT，地磁傾角為60°，地磁偏角為 10°，重力加速度為 9.8 m/s2，海水的磁導率為4π×10-7H/m，海水電導率為5 S/m，采樣頻率為10 Hz。

4.1 仿真速度比較分析

設定海況等級為3級，對應標準風速為8.23 m/s，得到海浪波譜圖如圖4所示。從圖4可以看出，海浪波譜能量主要分布在0.5～1.5 Hz頻帶范圍和-π/2～π/2角度范圍內。

圖4 海況等級為3級海浪方向譜Fig.4 Ocean wave directional spectrum under sea state level 3

對圖4中3級海況下海浪方向譜進行等能量四叉樹離散化分解，圖5為網格能量為總能量的0.5%時方向和頻段的離散化結果。從圖5離散化結果可以看出，四叉樹等能量分解法在能量密度低的區域采樣稀疏，在能量密度高的區域采樣密集。

圖5 基于四叉樹方向譜等能量離散化Fig.5 Equal energy division of directional spectrum based on quadtree

為了客觀比較區間等分法和四叉樹分解法的仿真速度，對不同的能量百分比條件下區間等分法和四叉樹分解法的計算次數進行比較，結果如表1?？梢钥闯瞿芰康确值谋壤叫?，四叉樹分解算法與等間隔分解算法所需要的計算次數比例越小。

表1 仿真速度比較Table 1 Simulation speed comparison

4.2 時間統計特性及頻域特性分析

波浪傳播主方向為45°，磁探儀高度50 m，基于四叉樹分解的優化仿真算法仿真不同海況等級下磁探儀靜止條件下采樣磁異常時間歷程信號如圖6所示，并利用Welch功率譜計算方法進行譜分析得到功率譜如圖7所示。從圖7譜分析可知，海浪磁場能量隨海況的增長而迅速的增加。隨著海況的增長，中心頻率是逐漸向低頻方向移動的。這與海浪譜的分布特征是吻合的。

圖6 不同海況下靜止磁探儀采樣海浪磁場信號仿真Fig.6 Simulation on ocean wave generated magnetic field signals sampled by a staying magnetometer under different sea state levels

圖7 不同海況等級靜止磁探儀采樣信號功率譜Fig.7 Power spectrum of ocean wave generated magnetic field signals sampled by a staying magnetometer under different sea state levels

根據線性海浪理論，磁探儀靜止條件下采樣仿真結果應該能反映線性隨機海浪磁場外觀上和統計上的特征，整體統計特征表現為上下對稱，均值為零，其正態性偏度和峰度應為0和3。對圖6仿真數據進行統計，得到均值、偏度和峰度如表2所示。根據表2統計數據可見，海浪磁場仿真結果與線性隨機海浪外觀上和統計上的理論特征吻合。

表2 靜止磁探儀采樣時域統計檢驗Table 2 Time-domain statistical analysis of samples by stationary magnetometer

海況等級為4級條件，海浪傳播主方向60°。磁探儀飛行方向為45°，飛行高度50 m，基于四叉樹分解的優化仿真算法仿真不同速度下航空磁探儀探測的海浪磁場時間歷程信號如圖8，并利用Welch法得到頻譜如圖9所示。

根據圖7、9可以看出，隨著飛行速度的增加，磁探儀接收到的海浪磁場信號存在明顯的多普勒效應。磁探儀靜止條件下采樣的信號能量主要集中在0～0.4 Hz，而航空磁探儀運動速度50 m/s時接收信號的主要頻段集中在0.2～1.6 Hz，存在明顯的頻帶擴展和頻率移動。隨著航空磁探儀運動速度的增加多普勒效應越明顯，這與文獻［9］中實驗分析的結果是一致的。

圖8 不同飛行速度下運動磁探儀采樣海浪磁場信號仿真Fig.8 Simulation on ocean wave generated magnetic field signals sampled by a magnetometer moving with different speeds

圖9 運動磁探儀不同飛行速度采樣海浪磁場信號功率譜Fig.9 Power spectrum of ocean wave generated magnetic field signals sampled by a magnetometer moving with different speeds

圖10 運動磁探儀不同飛行角度采樣海浪磁場信號功率譜Fig.10 Power spectrum of ocean wave generated magnetic field signals sampled by a magnetometer moving with different angles

海況等級為4級，海浪傳播主方向60°，航空磁探儀運動速度80 m/s，基于四叉樹分解的優化仿真算法仿真運動磁探儀不同運動角度采樣信號的功率譜分析結果如圖10所示。從圖10分析可知，當飛行方向與海浪傳播方向接近時，海浪磁場信號頻率更加集中，頻帶范圍更窄，隨著飛向方向與海浪傳播方向夾角的增大，海浪磁場信號帶寬變窄，并且明顯向低頻方向擴展。

5 結束語

基于Weaver海浪磁場模型推導了地理坐標系下沿任意方向傳播海浪單頻重力波感應磁場數學模型，在此基礎上基于線性理論推導了實際海浪磁場數值模型。基于觀測海浪譜給出了實際海浪磁場數值仿真算法，在海浪頻率和海浪方向上更加真實地描述實際海浪，并且采用四叉樹理論對仿真算法進行了仿真速度優化。時頻域驗證和分析結果表明，該仿真算法仿真速度快，仿真結果與理論和實際情況吻合。該仿真算法能夠仿真遠離海岸任意精度的連續海洋波浪產生的磁場，可用于航空磁探儀海浪磁噪聲背景消除研究，可為進一步的海浪磁場實驗研究提供指導。

［1］NELSON J B.Aeromagnetic noise during low-altitude flights over the Scotian shelf［R］.Dartmouth:Defence Research ＆Development，2002:15-18.

［2］LILLEY F E，HITCHMAN A P，MILLIGAN P R，et al.Sea-surface observations of the magnetic signals of ocean swells［J］.Geophysical Journal International，2004，159(2):565-572.

［3］李洪平，張海濱.海洋背景磁場模擬計算及東中國海表層磁場分布［J］.海洋技術，2008，27(3):70-74.LI Hongping，ZHANG Haibin.Simulation of ocean background magnetic field and its distribution in the East China Sea［J］.Ocean Technology，2008，27(3):70-74.

［4］WEAVER J T.Magnetic variations associated with ocean waves and swell［J］.Journal of Geophysical Research，1965，70(8):1921-1929.

［5］OCHADLICK A R.Experimental demonstration of weaver’s model of magnetic fields from ocean waves［R］.Washington，DC:Naval Air Development Center，1980:10-12.

［6］閆曉偉，閆輝，肖昌漢.海浪感應磁場矢量的模型研究［J］.海洋測繪，2011，31(6):8-11.YAN Xiaowei，YAN Hui，XIAO Changhan.Research on model of induced magnetic vector of ocean waves［J］.Hydrographic Surveying and Charting，2011，31(6):8-11.

［7］張揚，康崇，呂金庫.有限深海域的海浪感應磁場［J］.中國慣性技術學報，2012，20(5):594-595.ZHANG Yang，KANG Chong，LYU Jinku.Inductive magnetic field with deep ocean waves［J］.Journal of Chinese Inertial Technology，2012，20(5):594-595.

［8］JOCELYN F.Realistic simulation of ocean surface using wave spectra［C］//Proceedings of the First International Conference on Computer Graphics Theory and Applications，Portugal:CCSD Press，2006:76-83.

［9］唐勁飛，龔沈光，林春生.經典海浪譜下運動飛行器接收到的海浪磁場的噪聲［J］.海軍工程大學學報，2002，14(6):54-58.TANG Jinfei，GONG Shenguang，LIN CHunsheng.Wave generated magnetic noise received by moving airborne magnetometers under classical spectrum assumptions［J］.Journal of Naval University of Engineering，2002，14(6):54-58.

［10］YAAKOBI O，ZILMAN J，MILOH T.Detection of the electromagnetic filed included by the wake of a ship moving in a moderate sea state of finite depth［J］.Journal of Engineering Mathematics，2011，70(3):17-27.

［11］SEMKIN S V，SMAGIN V P.The effect of self-induction on magnetic field generated by sea surface waves［J］.Atmospheric and Oceanic Physics，2012，48(2):207-213.

［12］STEELE K E，EARLE M D.Directional ocean wave spectra using buoy azimuth，pitch and roll derived from magnetic field components［J］.IEEE Journal of Oceanic Engineering，1991，16(4):427-433.

［13］STEELE K E，TENG C C，WANG D W.Wave direction measurements using pitch and roll buoys［J］.IEEE Journal of Oceanic Engineering，1992，19(4):349-375.

［14］郭利進，師五喜，李穎，等.基于四叉樹的自適應柵格地圖創建算法［J］.控制與決策，2011，26(11):1690-1693.GUO Lijin，SHI Wuxi，LI Ying，et al.Mapping algorithm using adaptive size of occupancy grids based on quad-tree［J］.Control and Decision，2011，26(11):1690-1693.

［15］劉揚，宮阿都，李京.基于數據分層分塊的海量三維地形四叉樹簡化模型［J］.測繪學報，2010，39(4):410-413.LIU Yang，GONG Adu，LI Jing.A model for massive 3D terrain simplification based on data block partition and quad-tree［J］.Acta Geodaetica et Cartographica Sinica，2010，39(4):410-413.

［16］曾凱，楊華，翟月，等.光電成像干擾圖像質量評估［J］.電子與信息學報，2011，33(9):2164-2166.ZENG Kai，YANG Hua，ZHAI Yue，et al.Quality assessment of photoelectric image interference［J］.Journal of E-lectronics＆ Information Technology，2011，33(9):2164-2166.

［17］尹長林，詹慶明，許文強，等.大規模三維地形實時繪制的簡化技術研究［J］.武漢大學學報:信息科學版，2012，37(5):556-558.YIN Changlin，ZHAN Qingming，XU Wenqiang，et al.Simplification technology for real-time rendering of largescale 3D Terrain［J］.Geomatics and Information Science of Wuhan University，2012，37(5):556-558.