由于声波在海水中具有较小的衰减，因此水声技术目前被广泛应用于水下探测和水下通信领域^[1-3]。但水下声传播的速度低、延时大，特别是在浅海环境中，声学系统面临诸多挑战，如多径传播、来自空气界面和海底的回波、人造噪声等以及靠近海岸的复杂传播路径。此外，声波无法突破海洋界面直接传播，这也限制了水声技术的发展和应用。水下光通信的传输速率很高，但是光损耗很大^[4-5]，其对海水的清洁度要求也很高。同时只有在收发设备都精确对准的情况下才能正常通信，因此实际动态海洋环境也限制了水下光通信技术的发展。电磁波在海水中传播时会产生较大的衰减，频率越高衰减越大。但同时水下电磁波传播速度快、传输延迟低、不受海洋环境噪声影响，尤其是能轻易穿透海水界面传播。近年来，电磁波在海水中传播时的诸多优点吸引了越来越多的中外研究人员的关注^[6-7]，同时水下电磁场的相关应用也越来越广泛^[8]。

对海水中电磁波传播的研究实际上是对分层导电媒质中电磁波传播的研究。长期以来，针对海洋分层媒质中电磁波传播的建模和特性分析，开展了多项理论研究。通常人们仅考虑空气-海水-海底3层介质，分析其中电磁波的传播。Wang等^[9]基于3层媒质模型，研究了海水中的偶极子在空气中的辐射场，并采用快速傅里叶变换法将电磁场表达式中Sommerfeld形式的积分展开成离散傅里叶级数的形式，而后采用数值方法求解；该方法能够快速、准确地得到几个波长以内的结果。Qu等^[10]基于3层媒质模型研究了浅海环境中上下界面对海水中电磁波传播的影响。然而传统的3层模型并没有考虑真实海水的非均匀性。真实海水中随深度变化的电导率对水下电磁波传播，尤其是跨海-空界面传播会产生较大的影响。

本文针对海水电导率随深度变化时，水下电偶极子在空气中辐射的场展开研究。在传统空气-海水-海底3层模型的基础上，将海水视为多层媒质，基于直接全局矩阵对海水中水平电偶极子在空气中产生的电磁场进行建模计算。并与广义反射矩阵计算得到电磁场分量进行了比较分析。

1 海水电导率对电磁波传播的影响

海水作为一种导电媒质，它的传播常数与波矢量的关联方程式^[11]为：

$ \gamma=\mathrm{j} k=\sqrt{\mathrm{j} \omega \mu_{s} \sigma_{s}-\omega^{2} \mu_{s} \varepsilon_{s}}=\omega \sqrt{-\mu_{s} \varepsilon_{e}} $

(1)

式中：ε_e=ε_s-jσ_s/ω为海水中的复介电常数；ε_s=ε_rε₀；σ_s为海水中的电导率；μ_s为磁导率；ε_r为相对介电常数。由γ=α+jβ与式(1)可以得出：

$ \alpha=\omega \sqrt{\mu_{s} \varepsilon_{s}\left(\sqrt{\left(1+\left(\sigma_{s}{}^{2}\right) /\left(\varepsilon_{s}{}^{2}\right)\right)-1}\right) / 2} $

(2)

$ \beta=\omega \sqrt{\mu_{s} \varepsilon_{s}\left(\sqrt{\left(1+\left(\sigma_{s}{}^{2}\right) /\left(\varepsilon_{s}{}^{2}\right)\right)+1}\right) / 2} $

(3)

式中：α为导电媒质的衰减常数，表示电磁波沿传播方向衰减快慢的程度物理量；β称为相位常数，表示电磁波在传播过程中相位变化的物理量。由于σ/ωε>>1，所以可以认为海水为良导体，由此也可以将式(2)和式(3)进行近似简化得出：

$ \alpha \approx \beta \approx \sqrt{\pi f \mu_{s} \sigma_{s}} $

(4)

$ P L=-20 \lg \left(\mathrm{e}^{-\alpha z}\right) $

(5)

由式(5)可以看出无限大海水中电磁波振幅随传播距离的增加而呈指数衰减。衰减速度的快慢与衰减常数α的大小有着直接的关系。衰减常数越大衰减速率就越快，反之则衰减速率变慢。而电导率与衰减常数正相关，电导率的大小将直接影响电磁波在海水中传播时的衰减。

实际分层海洋环境下，电磁波跨界面传播时虽然不再满足指数衰减的规律，但电导率依然对电磁波的传播具有较大的影响。

2 传播模型的建立与求解 2.1 电磁波跨海—空界面传播模型

通常，海水中的辐射源在空气中产生电磁场时的3层模型结构^[12]如图 1所示(这里以水平电偶极子为例)，其中各层媒质的电磁参数如表 1所示。

实际上，由于海水并非空间均匀媒质，其电磁参数会随时间、空间变化，这点在前期的实验中也有所发现。这些参数的变化，尤其是电导率，会影响海水中及跨界面电磁波的传播。由于海水中电导率随深度变化尤为明显，如图 2所示，根据不同深度电导率变化，可以将变电导率的海水分为多层。这样，模型才更加符合实际海洋环境，得到的电磁场计算结果和传播特性分析才更加准确。

2.2 传播模型的求解

当电偶极子位于多层媒质中时，可以将其在媒质中产生的辐射场看作许多TE和TM波分量的线性叠加。设定所有分层区域中的介质都是各向同性的，其中偶极子天线位于第m层海水中，将柱坐标原点设置在偶极子的中心位置，则模型中任一接收点处电磁场的坐标为(ρ, φ, z)，其各分量强度计算公式为^[13]：

$ \begin{aligned} &E_{l \varphi}=-\int_{-\infty}^{\infty} \mathrm{d} k_{\rho} \frac{\mathrm{i} k_{l z}}{k_{\rho}^{2} \rho}\left(A_{lz} \mathrm{e}^{\mathrm{i} k_{l z} z}-B_{l} \mathrm{e}^{-\mathrm{i} k_{lz}{z}}\right) \mathrm{H}_{n}^{1}\left(k_{\rho} \rho\right) \sin \varphi+ \\ &\ \ \int_{-\infty}^{\infty} \mathrm{d} k_{\rho} \frac{-\mathrm{i} \omega \mu_{l}}{k_{\rho}}\left(C_{l} \mathrm{e}^{\mathrm{i}k_{l z}{z}}-D_{l} \mathrm{e}^{-\mathrm{i} k_{lz}{z}}\right) \mathrm{H}_{n}^{(1)^{\prime}}\left(k_{\rho} \rho\right) \sin \varphi \end{aligned} $

(6)

$ \begin{aligned} &H_{l \varphi}=\int_{-\infty}^{\infty} \mathrm{d} k_{\rho} \frac{\mathrm{i} k_{l z}}{k_{\rho}^{2} \rho}\left(C_{l} \mathrm{e}^{\mathrm{i} k_{lz}{z}}-D_{l} \mathrm{e}^{-\mathrm{i} k_{l z} z}\right) \mathrm{H}_{n}^{1}\left(k_{\rho} \rho\right) \cos \varphi+ \\ &\ \ \int_{-\infty}^{\infty} \mathrm{d} k_{\rho} \frac{\mathrm{i} \omega \varepsilon_{l}}{k_{\rho}}\left(A_{l} \mathrm{e}^{\mathrm{i} k_{lz} z}+B_{l} \mathrm{e}^{-\mathrm{i} k_{l z} z}\right) \mathrm{H}_{n}^{(1)^{\prime}}\left(k_{\rho} \rho\right) \cos \varphi \end{aligned} $

(7)

式中：$k_{l}=\sqrt{k_{\rho}{ }^{2}+k_{l z}{ }^{2}}$；k_l为传播常数；μ_l为磁导率；ρ为接收点与天线之间水平距离；ω为发射信号角频率；H_n¹(k_ρρ)为第1类n阶汉克尔函数；H_n^(1)′(k_ρρ)为H_n¹(k_ρρ)对k_ρρ的导数；A_l、B_l分别表示沿正z和负z方向TM波的幅度系数；C_l、D_l分别表示沿正z和负z方向TE波的幅度系数。从中可以看出，只要求出了每层的幅度系数，就得到了每层的电场和磁场。

当分层有N层时，根据介质交界处的边界条件上要求的切向电场分量在所有ρ和ϕ处连续可以得到4(N-1)个方程，由于没有波来自无限远处，所以未知量个数也是4(N-1)。选取A_l、B_l分量，根据相邻两层之间电磁场分量幅度系数的递推关系式，使用矩阵形式表示并简化为：

$ \begin{aligned} &{\left[\begin{array}{c} A_{l-1} \\ B_{l-1} \end{array}\right]=\left[\begin{array}{l} \mathrm{e}^{\mathrm{j}\left(k_{(l-1) z}-k_{l z}\right) d_{(l-1)}} R_{l(l-1)}^{\mathrm{TM}} \mathrm{e}^{\mathrm{j}\left(k_{l z}-k_{(l-1) z}\right) d_{(l-1)}} \\ R_{(l-1) l}^{\mathrm{TM}} \mathrm{e}^{\mathrm{j}\left(k_{(l-1) z}+k_{l z}\right) d_{(l-1)}} \mathrm{e}^{-\mathrm{j}\left(k_{l z}-k_{(l-1) z}\right) d_{(l-1)}} \end{array}\right] \times} \\ &\ \ \ \ \ \ \frac{1}{2}\left[\frac{\varepsilon_{l}}{\varepsilon_{l-1}}+\frac{k_{l z}}{k_{(l-1) z}}\right]\left[\begin{array}{c} A_{l} \\ B_{l} \end{array}\right]=\left[\begin{array}{cc} Q_{l}^{-} & P_{l}^{-} \\ Q_{l}^{+} & P_{l}^{+} \end{array}\right]\left[\begin{array}{c} A_{l} \\ B_{l} \end{array}\right] \end{aligned} $

(8)

在式(8)中，各层分界面的菲涅尔反射系数为：R_l(l-1)^TM=(ε_l-1k_lz-ε_lk_(l-1)z)/(ε_l-1k_lz+ε_lk_(l-1)z)，图 2中根据电导率随深度的变化将海水沿波传播的z方向进行了分层，这里将每层涉及的幅度系数在图 3中给出。

2.2.1 广义反射矩阵

通常，根据式(8)中相邻层幅度系数计算公式，可以逐层计算得到各层的幅度系数^[14]:

$ \left[\begin{array}{l} A_{1} \\ B_{1} \end{array}\right]=\left[\begin{array}{ll} Q_{1}^{-} & P_{1}^{-} \\ Q_{1}^{+} & P_{1}^{+} \end{array}\right]\left[\begin{array}{cc} Q_{2}^{-} & P_{2}^{-} \\ Q_{2}^{+} & P_{2}^{+} \end{array}\right] \cdots\left[\begin{array}{cc} Q_{m-1}^{-} & P_{m-1}^{-} \\ Q_{m-1}^{+} & P_{m-1}^{+} \end{array}\right]\left[\begin{array}{c} A_{m} \\ B_{m} \end{array}\right] $

(9)

由于水平电偶极子在第m层中，所以m层中电磁波的幅度系数还与电偶极子的激励幅度和结构相关，因此第m层中幅度系数应为：

$ A_{m}^{+}=A_{m}^{-}+E_{\mathrm{hed}} ; B_{m}^{+}=B_{m}^{-}+E_{\mathrm{hed}} $

(10)

式中：E_hed= iIlk_ρ²/(8πε_m)，A_m⁺和B_m⁺为第m层中z≥0区域中的幅度系数，A_m^-和B_m^-为第m层中z≤0区域中的幅度系数。结合式(8)和式(10)，可以得到第m层中所有幅度系数。

当使用数值计算方法对幅度系数进行求解时不可避免会引入一定量计算误差，设定计算误差为η，由式(8)中P_l^-/Q_/⁺=-e^{-j2k_(l-1)zd_(l-1)}>1，P_l⁺/Q_/^-=e^{-j2k_(l)zd_(l-1)}>1可见矩阵中P_l^-、P_l⁺绝对值较大，Q_l^-、Q_l⁺较小，计算得空气中电磁场的幅度系数为：

$ \begin{gathered} {\left[\begin{array}{c} \hat{\boldsymbol{A}}_{1} \\ \hat{\boldsymbol{B}}_{1} \end{array}\right]=\left[\begin{array}{cc} Q_{1}^{-}+\eta & P_{1}^{-}+\eta \\ Q_{1}^{+}+\eta & P_{1}^{+}+\eta \end{array}\right] \ldots\left[\begin{array}{cc} Q_{m-1}^{-}+\eta & P_{m-1}^{-}+\eta \\ Q_{m-1}^{+}+\eta & P_{m-1}^{+}+\eta \end{array}\right] \times} \\ {\left[\begin{array}{c} \boldsymbol{A}_{m}^{+} \\ \boldsymbol{B}_{m}^{+} \end{array}\right] \approx\left[\begin{array}{c} \left(P_{1}^{-} \cdots P_{m-1}^{-}\right) \eta & \left(P_{1}^{-}+P_{2}^{-}\right) \cdots\left(P_{m-2}^{-}+P_{m-1}^{-}\right) \eta \\ \left(P_{1}^{+} \cdots P_{m-1}^{+}\right) \eta & \left(P_{2}^{-} \cdots P_{m-1}^{-}\right) \eta \end{array}\right] \times} \\ {\left[\begin{array}{c} A_{m}^{+} \\ B_{m}^{+} \end{array}\right]+\left[\begin{array}{l} A_{1} \\ B_{1} \end{array}\right]} \end{gathered} $

(11)

式中：${\widehat A_1}$、${\widehat B_1}$为经过数值计算得到的幅度系数；A₁、B₁为真实的幅度系数。经过误差的层层累积，层数较多时计算得到矩阵与真实值相差较大。同时当d较大时，矩阵中P_l^-、P_l⁺与Q_l^-、Q_l⁺相差较大，在数值运算时经过多次相互乘积相加，不可避免会引起存在大数吃小数的现象，这将会导致求解得到电磁场分量强度不稳定。尤其在本文的研究中，将海水分成多层后，传统的基于广义反射矩阵的幅度系数计算方法不再适用。

2.2.2 直接全局矩阵

直接全局矩阵多用于声场研究中多层传播媒质时对每层介质的声反射系数进行求解^[15-16]，其主要原理是：根据相邻层之间的边界条件，构造矩阵将其结合起来组成方程组，对方程组进行求解后得到反射系数。这样保证了数值计算的稳定性，避免逐层误差传递。

由于界面处电场切向分量相等，相邻2层幅度系数之间关系如式(12)。其中，当分层厚度比较大并且所有层共用一个原点时，直接全局矩阵中的系数相差几个数量级。这将直接导致求解矩阵的奇异。但是通过将层的上界面作为下行波的原点，层的下界面作为上行波的原点，就确保了所有指数在介质层中都为负。且相互之间相差较小，避免了矩阵奇异的产生。

$ \left[\begin{array}{cccc} -k_{l z} \varPhi_{l}^{-} & k_{l z} & k_{(l+1) z} & -k_{(l+1) z} \varPhi_{l}^{+} \\ \varepsilon_{l z} \varPhi_{l}^{-} & \varepsilon_{l z} & \varepsilon_{(l+1) z} & -\varepsilon_{(l+1) z} \varPhi_{l}^{+} \end{array}\right]\left[\begin{array}{c} \bar{B}_{l} \\ \bar{A}_{l} \\ \bar{B}_{l+1} \\ \bar{A}_{l+1} \end{array}\right]=\left[\begin{array}{c} 0 \\ 0 \end{array}\right] $

(12)

式中：

$ \begin{gathered} \bar{A}_{l}=A_{l} \mathrm{e}^{\mathrm{j} k_{l z} d_{l}}, \quad \bar{B}_{l}=B_{l} \mathrm{e}^{-\mathrm{j} k_{lz} d_{l-1}} \\ \bar{A}_{l+1}=A_{l+1} \mathrm{e}^{\mathrm{j} k(l+1) z^{d}(l+1)}, \quad \bar{B}_{l+1}=B_{l+1} \mathrm{e}^{-\mathrm{j} k_{(l+1) z} d_{l}} \\ \varPhi_{l}^{+}=\mathrm{e}^{\mathrm{j} k_{(l+1) z}\left(d_{l}-d_{l+1}\right)}, \varPhi_{l}^{-}=\mathrm{e}^{-\mathrm{j} k_{l z}\left(d_{l}-d_{l-1}\right)} \end{gathered} $

将水平电偶极子所在位置看作一个虚拟分界面，由式(10)得到A_m⁺、B_m⁺和A_m^-、B_m^-的求解矩阵为：

$ \left[\begin{array}{ccc} \mathrm{e}^{-\mathrm{j} k_{m e} d_{m-1}} & 0 & -1 & 0 \\ 0 & 1 & 0 & -\mathrm{e}^{-\mathrm{j} k_{m e} d_{m}} \end{array}\right]\left[\begin{array}{c} \bar{B}_{m}^{+} \\ \bar{A}_{m}^{+} \\ \bar{B}_{m}^{-} \\ \bar{A}_{m}^{-} \end{array}\right]=\left[\begin{array}{c} E_{\text {hed }} \\ E_{\text {hed }} \end{array}\right] $

(13)

将(12)式和(13)式联立即可得到一个全局矩阵。求解后得到各层的幅度系数。在这期间也会产生一定计算误差η，其中空气中电磁场幅值系数为：

$ \hat{A}_{1}=\left(\bar{A}_{1}+\eta\right) \mathrm{e}^{-\mathrm{j} k_{l z}{d_{1}}}=A_{1}+\eta \mathrm{e}^{-\mathrm{j} k_{l z}{d}_1}, \hat{B}_{1}=0=B_{1} $

(14)

式中A₁、B₁为真实的幅度系数。对比式(11)与式(14)，直接全局矩阵计算得到幅值系数的误差相对较小，且不随海水分层累积传递。

3 仿真与结果分析

利用直接全局矩阵的方法，本文对真实海洋环境参数下水平电偶极子在空气中产生的电磁场进行仿真计算，其中设定海水中水平电偶极子的电偶极矩为1 Am，频率为15 kHz，发射深度15 m，仿真中空气接收高度20 m、接收角度均为φ=45°。

3.1 非均匀电导率对电磁波跨海—空界面传播的影响

图 4中，海水垂直深度上的电导率为2018年夏季在南海某站位使用conductivity-temperature-depth system(CTD)实测数据。其中，在整个深度范围内，该站位海水电导率并非一个恒定值。在10 m深度范围内电导率基本不变，超过10 m后基本保持为负梯度变化，在30 m左右时再次出现拐点。其实测最大电导率为5.59 S/m，最小电导率为4.93 S/m。

根据不同深度的电导率将海水进行分层，并使用直接全局矩阵求解得到多层介质下海水中水平电偶极子在空气中产生电磁场分量如图 5所示。除此之外，图 5还给出了假设海水电导率恒定为4.93 S/m和5.59 S/m时的结果。经比较可以看出：当海水电导率非恒定时，得到的电磁场分量强度介于海水电导率恒定为最大和最小值时得到电磁场强度之间，与两者最大相差4.8 dB。由于本站位10 m深度范围内电导率基本保持为5.59 S/m，所以计算结果与恒定时近似。

3.2 直接全局矩阵与广义反射矩阵对比

在2.2.1节中，分析了广义反射矩阵在计算多层导电媒质中电磁场时存在误差传递的现象，本节将通过仿真进行验证。

根据电导率变化，将海水分别设置为5层和10层，同时使用广义反射矩阵和直接全局矩阵对分层海水情况下电磁波跨海—空界面传播模型进行仿真。通过对图 6中仿真结果对比，可以发现使用广义反射矩阵计算得到的电场分量在5层时与直接全局矩阵计算得到结果相近，但当层数升高到10层时电场分量曲线之间差距变大。这说明广义反射矩阵计算方法的确会随着介质分层数量的增加而误差增大。

为了进一步说明直接全局矩阵相比广义反射矩阵计算时的稳健性，对电磁场分量中的幅度系数B₁随积分变量k_ρ的变化曲线画出如图 7。理论上来说，由于本文给出的传播模型中空气为无限大上半空间，幅度系数B₁=0。从图 7中可以看出，直接全局矩阵得到的B₁在所有k_ρ处都为0，与理论值相符；但广义反射矩阵得到B₁不全为0，而且随着层数的增加，B₁也在变大即误差也在变大。所以使用广义反射矩阵计算得到幅度系数相对于直接全局矩阵有较大误差。

4 结论

1) 本文建立了垂直深度上电导率非均匀的海水中电磁波跨海—空界面传播模型，提出了基于直接全局矩阵的多层媒质中电磁波传播求解方法。并对变电导率海水中水平电偶极子在空气中产生的电磁场进行了建模和分析。结果表明，垂直深度海水电导率的非均匀性将对水下电磁波跨海—空界面传播带来影响，因此在理论建模过程中，需要将海水视为多层导电媒质。

2) 文中定性对比了广义反射矩阵和直接全局矩阵的求解误差，分析了多层媒质中广义反射矩阵的误差累积过程。通过仿真获得了不同分层数量下电磁场强度和幅度系数，对比发现了采用直接全局矩阵求解多层媒质中电磁波传播将比广义反射矩阵更加准确和稳健。

3) 文中提出的求解模型和方法对于利用电磁波实现水下目标探测和水下数据无线传输具有重要的理论价值，该模型和方法还可以应用到深地探测等涉及分层媒质中电磁波传播的领域。