基于四叉树的海浪磁场快速仿真算法

熊雄，杨日杰，韩建辉，郭新奇

(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 海浪磁场数值仿真及分析

基于海浪磁场快速仿真算法，对不同条件下海浪磁场进行数值仿真计算，并分析仿真性能。基本的仿真条件为:地磁场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.