兩步SVD在極低頻電場透地通信中的應用

吳云具，孫 暢，提文猛，宋新昌

（1.清江創新中心，湖北 武漢 430076；2.中國船舶集團有限公司第七一〇研究所，湖北 宜昌 433003）

0 引言

極低頻電場透地通信[1]（Through-the-earth ELF Communications）同陸地可控源電磁法類似，利用人工交流電源通過接地發射電極以一定頻率對地或在海水中供電，產生時變電磁場，在地表或海底布設電磁場測量設備測量該電磁場，進而實現通信的目的。

各接收點為了獲得準確的信號響應以及傳遞的信息，需要對所采集的電場數據進行處理分析。然而，在開展極低頻電場透地通信時，各接收點采集的電場數據中包含各種復雜的噪聲，如采集設備的本底噪聲、人文噪聲以及其它干擾源產生的噪聲等[2]。當收發距較大時，采集獲得的電場數據具有較低的信噪比，并且有效信號可能淹沒于復雜的噪聲中，嚴重干擾所傳遞的信息。因此，需對大收發距、低信噪比的極低頻電場透地通信數據進行濾波處理。

目前對極低頻電場透地通信數據進行濾波處理的主要方法有頻率域濾波、小波閾值去噪等方法。然而，在有效信號與噪聲成分頻率相近或者噪聲能量較強時，此時頻率域濾波或者小波閾值濾波等不易區分有效信號和噪聲，需要多次嘗試不同的濾波參數，不易達到理想的噪聲壓制效果。奇異值分解[3-4]是一種重要的矩陣分解方法，在數據分離[5-7]、圖像處理[8-10]及信號處理[11-15]等方面具有廣泛的應用，并且在UUV集群反水雷搜索中也得到了應用[16]。奇異值分解應用到信號處理時，通過將數據構造成Hankel 矩陣[17]，分解得到對角奇異值矩陣和2 個正交矩陣。不同的奇異值分別對應數據中不同的成分，一般情況下，信號所對應的奇異值較大，噪聲所對應的奇異值較小。因此，將噪聲對應的奇異值與相應的奇異向量賦零值，而保留信號對應的元素，重構數據后便達到了噪聲壓制的目的。然而，在有效信號與噪聲成分頻率相近或者噪聲能量較強時，經過一次奇異值分解濾波處理，往往只能壓制部分噪聲，同樣無法達到理想的噪聲壓制效果。

本文提出了一種兩步奇異值分解濾波方法，同樣假設有效信號與噪聲相互獨立，并且有效信號具有較大的奇異值。通過2 次奇異值分解，分別刪除較小的奇異值且保留較大的奇異值來提取電場數據中的信號成分，從而壓制其噪聲成分。將該方法應用到某水庫實測數據中，結果表明該方法可以明顯壓制數據中的噪聲成分，極大提高了數據的信噪比。

1 奇異值分解

奇異值分解是一種矩陣正交化方法，在一個行或列線性相關的矩陣左右分別乘上一個正交矩陣，便可將原矩陣中的行或列變為線性獨立的[18]。對任意大小為m×n的矩陣A，其秩為r，那么存在2 個大小分別為m×m和n×n標準化的正交矩陣U和V，以及大小為m×n對角矩陣D，滿足[13]：

其中：

式中：Σ=（λ1λ2···λr），λ1≥λ2≥···≥λk≥···≥λr，λk為奇異值；對角矩陣D又稱為奇異值矩陣。

矩陣的奇異值具有良好的穩定性，不同的奇異值可反映矩陣不同的特征。

2 奇異值分解濾波算法

在對信號進行濾波處理時，可以通過奇異值分解篩選出信號所對應的奇異值，從而達到濾波去噪的目的。然而，通過各種傳感器采集得到的信號響應并不是矩陣而是時間序列，因此需要先將信號的時間序列轉換成矩陣，然后才能進行奇異值分解。

假設有時間序列X=（x1x2···xN），將其構造成Hankel 矩陣：

式中：1＜m＜N；1＜n＜N。

考慮到時間序列包含元素的數量很大，為了避免Hankel 矩陣中每行的元素過多，在構造Hankel矩陣時相鄰2 行重疊的元素可以少于n- 1個。

對Hankel 矩陣進行奇異值分解，便可得到正交矩陣U和V以及奇異值矩陣D。實測信號中往往包含各種噪聲，考慮到信號和噪聲具有不同的特性，因此它們具有不同的奇異值。一般情況下，信號具有較大的奇異值，而噪聲具有較小的奇異值。因此，保留奇異值矩陣D中較大的奇異值，將較小的奇異值賦零，然后按照奇異值分解的逆過程對其進行重構，便可得到濾波后的信號。通常，將特征值由大到小排列，較大特征值之和占特征值總和達到80%，保留這部分特征值并將剩余較小特征值賦零。

然而，在信號和噪聲的特征相近時，它們具有相近的奇異值，進行一次奇異值分解濾波并不能達到令人滿意的噪聲壓制的效果。這時，將濾波后的信號構造成與第一次奇異值分解濾波時大小不同的Hankel 矩陣，再進行一次奇異值分解濾波，便可進一步壓制與信號特征相近的噪聲。

3 應用

3.1 試驗情況

2022年9月，在某水庫開展了極低頻電場透地通信試驗。該試驗觀測系統布設如圖1所示，其中發射源為電偶極源，布設在岸邊；3 個測點均設置2 個通道，可采集2 個互相垂直的水平電場分量，布設在水庫底部。發射源和3 個測點的參數信息如表1所示，其中測點的布設角度為通道1 的方向角。

表1 發射源和測點參數信息Table 1 Layout information of a transmitter and three receivers

圖1 觀測系統布設示意圖Fig.1 Layout of observation system

試驗期間，分別發射8 Hz、8 Hz 和4 Hz 的正弦信號，每次發射持續5 min，采集系統的采樣率為1 024 Hz，為了處理方便，我們將其重采樣為32 Hz。選取信號發射期間3 個測點采集的電場數據，來驗證兩步奇異值分解濾波方法的有效性。

3.2 數據處理

通過兩步奇異值分解濾波方法對3 個測點進行濾波處理，以A 測點為例，說明該方法的濾波效果。A 測點通道1 的時間序列如圖2（a）所示，紅色曲線為原始電場數據，藍色曲線為1 次奇異值分解濾波后的電場數據，青色曲線為2 次奇異值分解濾波后的電場數據，可以看出1 次奇異值分解濾波后噪聲大部分被壓制，2 次奇異值分解濾波后噪聲被進一步壓制，從時間序列中已可分辨出3 次信號發射時段。

圖2 兩步奇異值分解對A 測點通道1 的噪聲壓制效果Fig.2 Data processing results of two-step SVD filter of channel 1 of receiver A

圖2（b）為A 測點通道1 電場數據濾波前的時頻振幅譜，可知在7～9 Hz 之間持續存在條帶狀異常，這與舒曼共振的頻率（7.83 Hz）吻合。此外，A 測點通道1 原始數據中還存在貫穿全頻率的線狀異常，以及能量較強的背景噪聲等。圖2（c）為一次奇異值分解濾波后的時頻振幅譜，7～9 Hz之間持續存在條帶狀異常以及背景噪聲得到壓制，但貫穿全頻率的線狀異常仍然存在。圖2（d）為二次奇異值分解濾波后的時頻振幅譜，7～9 Hz 之間持續存在條帶狀異常以及背景噪聲得到進一步壓制，并且貫穿全頻率的線狀異常也得到明顯壓制，使得3 次信號發射時的信號響應更加明顯。

A 測點通道2 的時間序列如圖3（a）所示，紅色曲線為原始電場數據，藍色曲線為一次奇異值分解濾波后的電場數據，青色曲線為二次奇異值分解濾波后的電場數據，可以看出一次奇異值分解濾波后噪聲大部分被壓制，二次奇異值分解濾波后噪聲被進一步壓制。圖3（b）為A 測點通道2 電場數據濾波前的時頻振幅譜，相比于通道1，通道2電場數據中3 次信號發射時的信號響應較小，通道2 與發射方向的夾角較大，有效接收極間距小于通道1。此外，通道2 中同樣在7～9 Hz 之間持續存在條帶狀異常、貫穿全頻率的線狀異常以及能量較強的背景噪聲等。圖3（c）為一次奇異值分解濾波后的時頻振幅譜，7～9 Hz 之間持續存在條帶狀異常以及背景噪聲得到壓制，但貫穿全頻率的線狀異常仍然存在。圖3（d）為二次奇異值分解濾波后的時頻振幅譜，7～9 Hz 之間持續存在條帶狀異常以及背景噪聲得到進一步壓制，并且貫穿全頻率的線狀異常也得到明顯壓制，使得3 次信號發射時的信號響應更加明顯。

為了定量地評價兩步奇異值分解濾波方法對電場數據的改善效果，引入信噪比（Signal-to-Noise Ratio，SNR）作為評價參數，其計算公式如下：

式中，Ps和Pn分別為信號和噪聲的功率。如果濾波方法有效，那么相比于濾波前，濾波后的數據應當具有更高的SNR。

這里，為了將SNR 這一評價參數應用到實際數據中，將式（6）中的Ps定義為某一頻率信號響應所在頻率的平均功率，將Pn定義為該頻率附近頻率的平均功率。實際上，這種方法計算得到的SNR 明顯大于嚴格按照定義計算得到的SNR，不過仍然可以作為噪聲壓制效果的評價參數。

截取3 次信號發射所在時間段3 個測點的電場數據，分別計算其濾波前后的SNR，如表2所示。經過兩步奇異值分解濾波之后，3 次信號發射時的信號響應所在時間段的SNR 均得到顯著提高。需要說明，B 測點通道2 與發射方向幾乎垂直，有效接收極間距很小，導致2 次8 Hz 發射信號均完全淹沒在噪聲中，已無法還原。這表明，在原始數據具有一定信噪比的前提下，兩步奇異值分解濾波方法可以有效壓制噪聲，從而提高數據的信噪比。

表2 3 個測點電場數據兩步奇異值分解前后的SNRTable 2 SNR of electric field data of 3 receivers before and after two-step SVD filtering

4 結束語

本文基于奇異值分解，提出了一種兩步奇異值分解濾波方法，用以改善極低頻電場透地通信數據的質量。某水庫3 個測點的實測數據處理結果表明：

在極低頻電場透地通信原始數據具有一定信噪比的情況下，兩步奇異值分解濾波方法在突出有效信號的同時可以有效壓制電場數據中的干擾噪聲，使得濾波后數據的信噪比明顯提升，從而極大地改善數據質量。