基于遞推最小二乘算法的光纖振動定位系統

羅義軍，尹 棋，李 勁

(1.武漢大學 電子信息學院，武漢 430072; 2.武漢紡織大學 電子與電氣工程學院，武漢 430200)

引 言

基于雙馬赫-曾德爾的分布式光纖振動定位系統對振動感知和應變十分敏感，具有功率損耗小、不受電磁干擾、造價低等優點，在油氣管道檢測報警以及周界安防中大量應用[1-4]，其中擾動定位算法以開環互相關算法和閉環最小均方(least-mean square,LMS)算法為主[5-9]。互相關以其開環結構對信號毛刺，電路噪聲敏感，擾動定位穩定性不足[5-10]。LMS的最優準則是使濾波器輸出與期望信號誤差平方的數學期望最小，根據輸入數據的長期統計特性得到的對一類數據的最佳濾波系數[11]。然而在判斷一次定位時，已知的僅是采集的一組振動數據，并且對于LMS，濾波器階數愈高，步長因子和輸入信號功率愈大，其失調系數也愈大[11]，會直接影響定位結果。

采用閉環遞推最小二乘(recursive least-squares,RLS)是對每一組已知數據尋求精確的最佳濾波器系數。因為RLS的最優準則是使濾波器輸出與期望信號誤差平方累計和最小，RLS收斂速度快于LMS[12-13]，避免了固定步長LMS收斂速度和穩態均方誤差相互制約的矛盾以及變步長的LMS待調參量多且調節困難的問題。由于濾波器收斂后，濾波器系數波形將近似于sinc函數，且兩路信號延遲點數在峰值位置獲得[7]。借鑒早遲門同步判決最佳采樣位置的原理，以鑒相曲線的近零點值位置作為定位的結果。現場可編程門陣列(field-programmable gate array,FPGA)可以支持高速多路并行運算，非常適合處理數據流，因此系統中定位算法由FPGA實現。

1 光纖定位系統原理闡述

1.1 雙馬赫-曾德爾定位原理

圖1中右側包含連續激光源；O1,O2,O3為3dB光耦合器；3根分別長為L1,L2,L3的鎧裝單模光纖；光電探測器D1,D2；環形器C1,C2。左側是對光電探測器接收信號的處理流程。在沒有振動時，經O1耦合的兩路相干光同時到達D1,D2。而在長為L2的振動臂上某一點有入侵行為，將存在時間差Δt，且入侵位置P可表示為[4-9]：

(1)

Fig.1 Optical path model and signal processing flow

設采樣率為f, 單個采樣點表示的定位距離精度r由相鄰兩個采樣點之間的時長決定，r可以通過下式求得[8,10]：

(2)

式中,c為光速,n0為光纖折射率，約為1.46。

1.2 RLS算法及橫式濾波器估計時延原理

RLS自適應濾波器是指能夠根據環境的改變，采用累計誤差平方和為代價函數的準則，自適應反饋調節濾波器系數使得代價函數最小[11]。圖2所示為M階RLS橫式濾波器結構。給定一個輸入采樣序列x(n),x(n-1)…x(n-M+1)，并用其沖激響應w0,w1,…,wM-1來表征該濾波器。在某一離散時刻n，濾波器的輸出為y(n)，先驗估計誤差eerr用期望響應d(n)與y(n)之差表示，通過迭代反饋使得下一刻的輸出值y(n+1)不斷地逼近期望信號d(n+1)，誤差將迫近為0。圖中，Z-1表示信號延遲1拍。

Fig.2 Structure of RLS adaptive filters

假定n時刻輸入M×1維矢量X(n)=[x(n),x(n-1),…,x(n-M+1)]T，濾波器系數M×1維矢量W(n)=[w0,w1,…,wM-1]T，T表示矢量轉置；M×M維迭代矩陣C(n)。振動信號為平穩信號，取遺忘因子λ=1；正則化參量δ為小正實數，與輸入信號x(n)的信噪比有關[12-13]。對RLS采用前加窗法，算法流程見下[11-13]。

(1)初始n=0時刻值：

(3)

式中，I是M×M維單位矩陣；迭代過程對n=1,2…，更新增益矢量g(n):

(4)

(2)更新濾波器系數矢量W(n)：

W(n)=W(n-1)+

g(n)[d(n)-[X(n)]TW(n-1)]

(5)

(3)更新迭代矩陣C(n)：

C(n)=λ-1[C(n-1)-

g(n)[X(n)]TC(n-1)]

(6)

由于濾波器的輸入信號x(n)會遍歷濾波器的0～M-1級延遲位置。由疊加定理知，相對于x(n)的延遲信號x(n-τ)可由濾波器各級沖激響應w0,w1,…,wM-1,表示為：

x(n-τ)=w0x(n)+…+

wM-1x(n-M+1)

(7)

式中，τ為延遲點數。又根據卷積定理：

x(n-τ)=x(n)*δ(n-τ)

(8)

式中,δ(n)為采樣沖激響應函數。取：

(9)

式中，0≤i

(10)

即有限長橫向濾波器的系數wi與采樣函數sinc(i-τ)重合時，濾波器輸出近似為延遲信號x(n-τ),且系數最大值wτ的位置為延遲點數τ。

1.3 早遲門原理

尋找系數峰值的位置i=τ直接判別在噪聲存在的情況下易產生偏差[16]。假定不在峰值點對信號采樣，在i=τ-φ位置早采樣，在i=τ+φ位置遲采樣，早、遲采樣的絕對值都比峰值采樣值的絕對值小,φ表示偏離系數峰值位置的點數。再以遲采樣時刻值減去早采樣時刻值，得到鑒相曲線。由于收斂的系數相對于最佳采樣時刻i=τ是偶函數(sinc(i-τ))，早、遲采樣值的絕對值應該相等。因此，采樣峰值的位置就是在鑒相曲線為零的位置。

2 系統硬件平臺及定位算法實現

2.1 光纖定位系統硬件平臺

基于圖1左側電信號處理流程，搭建了圖3所示的硬件平臺。D1,D2是PIN激光管，將光信號轉換為電小信號，以高精度低噪聲的雙運放OPA2227放大5倍。經過截止頻率為100kHz的有源低通濾波器UAF42濾出有用振動信號。全差分運算放大ADA4930將單端信號轉換為差分信號，進一步地抑制噪聲影響，并放大2倍驅動后級高速模-數轉換器(analog-to-digital converter, ADC)。ADC選取ADI公司推出的低功率、16bit、雙通道、最大采樣率為125MHz的ADC9268，采樣時鐘定為10MHz。FPGA選取Altera公司的Cyclone Ⅳ系列高性價比低功耗、邏輯資源豐富的EP4CE115F23I7實現振動信號識別及定位算法。

Fig.3 Circuit block diagram of signal processing

2.2 實現定位的FPGA結構

圖4為處理ADC采集的D1,D2兩路信號FPGA程序框圖。兩路信號經過64階低通濾波后以兩路信號的方差判別振動起點，將RLS每次迭代獲取的濾波器系數經過滑動平均過程，降低噪聲對于后級早遲門判決最大系數位置的影響，最后將估計的延遲點數轉換為實際擾動位置。

Fig.4 FPGA block diagram

2.2.1 振動起點識別 光纖在靜止狀態下發生振動時會導致ADC采集的數據離散度迅速增大，可以取一滑動窗內的定點數據量的方差來判斷振動起點。方差獲取的simulink框圖(見圖5)，任意一路信號方差值連續多次超過閾值，則識別為振動的起點。之后所采集的數據均為振動數據，用于RLS定位算法。圖6為Quartus Ⅱ中計算方差的Verilog頂層例化模塊框圖。

Fig.5 Simulation block diagram of calculating the variance

Fig.6 FPGA block diagram of calculating the variance

2.2.2 RLS算法FPGA實現 在一段160m光路上進行定位測試，以10MHz(r=±10m精度)采樣后兩路信號最大延遲點數為16，取RLS橫向濾波器為M=20階實現定位。對于長距離的定位，可以相應地增加濾波器的階數，而這會影響算法收斂速度，也會成倍增加FPGA硬件資源消耗。可先對采集的數據進行下抽，降低定位精度，獲取擾動的大致距離區間。在此基礎上，以原來的采樣率進一步的細定位[17]。算法迭代流程中計算濾波器輸出y(n)，M×1維增益矢量g(n)，M×1維系數矢量W(n)及M×M維矩陣C(n)均是行矢量乘以列矢量的乘、加基本結構，因此FPGA中基本計算單元設計為1×M維行矢量乘以M×1維列矢量[18-19]。圖7為Quartus Ⅱ中基本計算單元的Verilog例化頂層模塊框圖。

Fig.7 FPGA block diagram of RLS basic computational unit

此外，迭代更新系數矢量W(n+1)及矩陣C(n+1)則是上一拍的舊值與校正值的迭代累加結構，其基本單元simulink結構見圖8。

Fig.8 Simulation block diagram of basic iterating unit

分析迭代流程可知，RLS的運算量主要集中于(4)式和(6)式中狀態量的計算和迭代更新[12]。RLS算法每次迭代需要3M2+3M次乘法，1次除法和2M2+2M次加減法，遞推的運算次數在O(M2)[12]。算法復雜度結構框圖如圖9所示。

Fig.9 Analysis of computational complexity for RLS algorithm

2.2.3 滑動平均處理及早遲門 RLS濾波器收斂后，其抽頭系數波形將近似于sinc函數主瓣峰，每次迭代更新后獲取的系數值也在濾波器理論最佳系數附近，波形不會發生本質性的變化。由于系數迭代從初始值開始至全部振動數據結束，每迭代一次便可以獲得一組濾波器系數。可以充分利用每次迭代獲取的系數值，讓其經過一滑動平均過程，在統計平均的意義上不改變波形本質且能抑制某一次迭代過程因為噪聲的隨機性而使系數波形隨機起伏[20]。圖10為滑動平均和早遲門的simulink結構圖。

Fig.10 Simulation block diagram of moving-average and early-late door

滑動平均后的系數取φ=1的早遲門相減得到鑒相曲線。圖11鑒相曲線計算零值位置的Verilog頂層例化模塊框圖。

Fig.11 FPGA block diagram of getting position by phase curve

2.3 定位結果處理

一次振動的時長在500ms以上，其中幅度較強的信號大約集中在最初200ms以內。當ADC采樣率為10MHz時，一次振動可以獲取幅度較強的信號點數約2×106個數據。RLS具有很快的收斂速率，約迭代1×105次(10ms采樣點數)后就可以得到良好的定位結果。因此FPGA設計時，將最初振動幅度較強的信號分為18段分別進行定位計算，記錄每一段的定位結果。而后剔除其中定位的最大值和最小值，剩余16段計算結果取平均，作為最終的定位結果。這樣，在不增加濾波器階數和ADC采樣率的情況下，提高定位精度同時也抑制噪聲而減少誤判概率。

3 實驗結果

基于圖1中雙馬赫-曾德爾光路框圖和圖3中的光纖振動硬件平臺，設計了一套光纖振動定位系統。在一段160m的光路上進行重復敲擊實驗，選用波長為1550nm，功率可調的穩定化臺式激光源，功率輸出范圍0.1mW～40mW。光路鋪設部分選用8芯單模鎧裝光纖(GYXTW-8B),并將光纖蛇形盤繞在小區圍欄上。

3.1 simulink中定位算法對比仿真

在simulink中對開環互相關，閉環最小均方LMS以及遞推最小二乘RLS進行定位對比仿真。選擇在90m附近敲擊光纖，上位機通過千兆網口接收ADC9268采集的振動數據作為仿真源。采樣率為10MHz(r=±10m精度)時，simulink中3種定位算法以50×104個采樣點計算振動位置。從圖12中的仿真結果可以看出，LMS收斂速度明顯慢于RLS；而互相關雖也得到振動位置，但隨著互相關累加點數增加，相關峰位置改變，穩定性差。比較之下，RLS在10×104個采樣點以內已經快速收斂出振動位置(9×10=90m),且定位收斂穩定。縱坐標為D1,D2兩路相關信號相差的延遲點數。

Fig.12 Results of positioning simulation

算法收斂后，RLS濾波器系數收斂波形及鑒相曲線如圖13所示，縱坐標分別為RLS濾波器系數的值和鑒相曲線的值。

Fig.13 Filter coefficient value and phase curve

3.2 FPGA定位結果

實驗中選擇在60m附近，以手握光纖、攀爬護欄、輕敲、重敲等模擬不同的實際擾動環境。圖14為某一次60m處攀爬護欄，用SignalTapⅡ抓取各模塊關鍵信號及ADC采集的振動波形。

在60m附近重復敲擊30次，記錄每一次上位機接收的定位結果，在MATLAB中繪制出圖15。從圖中可以大致看出，實際定位誤差在±6m。

Fig.14 Vibration wave in SignalTap II

Fig.15 Multiple tap positioning results near 60m

4 結 論

運用RLS、滑動平均和早遲門設計出一套光纖振動定位系統，并在FPGA硬件平臺上實現快速定位。對比于開環的互相關以及閉環LMS，RLS可以根據環境變換自適應快速的調整參量，針對每一組ADC采集數據得到不同的濾波器系數。并且RLS隨著迭代次數增加，單次定位結果穩定不會發生波動，定位收斂速度也是LMS的3倍左右，解決了開環互相關穩定性不足,避免了LMS收斂速度和穩態誤差相互制約的矛盾。由于RLS具有快速的收斂性，對振動數據分段定位處理，定位平均值作為最終定位結果，保證定位的穩定可靠，且系統精度保持在±6m。在實際應用中,包括刮風、下雨等大量非人為擾動因素存在的復雜環境，后續工作著重在于對不同振動信號的識別。