含沙量監(jiān)測(cè)的wavelet—Kalman多尺度融合研究

付立彬 劉明堂 王麗 秦澤寧 楊陽(yáng)蕊

摘要：為解決黃河含沙量監(jiān)測(cè)時(shí)傳感器易受環(huán)境因素影響的問(wèn)題，簡(jiǎn)述了音頻共振法測(cè)量含沙量的原理，探討了音頻共振傳感器的諧振頻率和含沙量監(jiān)測(cè)之間的關(guān)系，提出了貫序式wavelet-Kalman多尺度融合模型，對(duì)音頻共振傳感器的諧振頻率進(jìn)行小波分解，把含沙量數(shù)據(jù)組成貫序式數(shù)據(jù)塊進(jìn)行多尺度分析，提取含沙量信號(hào)序列中的突變值，建立了Kalman融合方程，將溫度信息作為控制信號(hào)，消除了環(huán)境因素對(duì)含沙量監(jiān)測(cè)的影響，并進(jìn)行了含沙量測(cè)量的反演和誤差分析。結(jié)果表明：貫序式wavelet-Kalman多尺度融合模型能夠有效地消除環(huán)境影響，提高系統(tǒng)測(cè)量的精度和穩(wěn)定性，平均絕對(duì)誤差為3.95kg/m2，均方根誤差為3.13kg/m3，比其他反演模型的誤差小。

關(guān)鍵詞：音頻共振法；貫序式；小波多尺度分析；卡爾曼融合；含沙量監(jiān)測(cè)

中圖分類號(hào)：P335 文獻(xiàn)標(biāo)志碼：A doi：10.3969/j.issn.1000-1379.2018.09.005

懸浮泥沙使得高含沙河流具有特殊的產(chǎn)流、匯流、產(chǎn)沙和輸沙規(guī)律，會(huì)造成河床的劇烈變化。黃河是典型的懸浮泥沙含量較高的河流，也是世界上最復(fù)雜、最難治理的河流之一。目前黃河等多沙河流上的泥沙采樣設(shè)備大多仍采用傳統(tǒng)的采樣器，部分重要河段只能獲取低量程、低精度的含沙量信息。為了更好地研究黃河流域土壤侵蝕機(jī)理和防治水土流失，就必須建立黃河高懸浮含沙量在線實(shí)時(shí)監(jiān)測(cè)系統(tǒng)。

比較先進(jìn)的含沙量監(jiān)測(cè)方法有射線法、聲學(xué)法、光學(xué)法、電容差壓法、浮力法、濁度法、遙感光譜分析法等。這些方法在實(shí)際應(yīng)用中有各自的優(yōu)勢(shì)，但也有各自的局限性。監(jiān)測(cè)含沙量時(shí)，傳感器一般是直接插人水體，采集的信號(hào)會(huì)受到水溫、流速和深度等因素影響。這些影響因素會(huì)給監(jiān)測(cè)系統(tǒng)帶來(lái)一定的干擾，造成傳感器時(shí)間序列信號(hào)中存在相應(yīng)的突變信息，使得測(cè)量系統(tǒng)運(yùn)行失穩(wěn)。因此，研究含沙量傳感器的線性輸入輸出特性變化并進(jìn)行含沙量信息多尺度分析具有重要的意義。

小波（wavelet）分析方法是研究水文信息序列多尺度變化特性的有效工具。在含沙量時(shí)間序列分析中，同樣可以運(yùn)用小波的多尺度分析方法，提取含沙量傳感器突變的輸入輸出特性信息?？柭↘al-man）濾波是根據(jù)測(cè)量數(shù)據(jù)變化不斷融合的過(guò)程，可減少環(huán)境因素的干擾。因此，本文提出基于貫序式（sequential）小波多尺度分析和卡爾曼濾波融合的含沙量監(jiān)測(cè)模型，采用音頻共振法來(lái)監(jiān)測(cè)黃河含沙量，并同時(shí)監(jiān)測(cè)含沙水體的環(huán)境信息，進(jìn)行含沙量數(shù)據(jù)融合處理，以提高測(cè)量系統(tǒng)的精度和穩(wěn)定性。

1 硬件平臺(tái)及試驗(yàn)數(shù)據(jù)

1.1 硬件平臺(tái)搭建

圖1為含沙量監(jiān)測(cè)系統(tǒng)的整體結(jié)構(gòu)。該系統(tǒng)將含沙量、水溫等傳感器集成在一起，可完成多通道信息采集與處理工作。處理后的數(shù)據(jù)通過(guò)ZigBee技術(shù)無(wú)線傳送給信息處理中心，信息處理中心通過(guò)仿真界面實(shí)時(shí)顯示測(cè)量結(jié)果。

1.2 試驗(yàn)材料和試驗(yàn)數(shù)據(jù)

試驗(yàn)材料來(lái)源于花園口水文站附近的黃土（含有細(xì)沙），同黃河泥沙具有一定的相似性，平均粒徑為176.42μm，中值粒徑為96.56μm。每次試驗(yàn)都記錄音頻共振傳感器的輸出值K和水溫T，同時(shí)用比重瓶獲取含沙量實(shí)測(cè)值。

從大量試驗(yàn)數(shù)據(jù)中經(jīng)過(guò)轉(zhuǎn)換、濾波、異常值剔除等環(huán)節(jié)后，本文選取92組數(shù)據(jù)塊（data block）進(jìn)行數(shù)據(jù)融合分析。其中第1組到第46組數(shù)據(jù)用來(lái)擬合，第47組到第92組數(shù)據(jù)用來(lái)驗(yàn)證。受試驗(yàn)條件限制，本文只進(jìn)行了水溫的環(huán)境變量監(jiān)測(cè)。

2 多尺度融合模型建立

2.1 諧振頻率的小波分解

音頻共振傳感器的一臂在共振激勵(lì)器驅(qū)動(dòng)下作簡(jiǎn)諧振動(dòng)，音頻共振傳感器的另一臂將產(chǎn)生共振，也作簡(jiǎn)諧振動(dòng)。當(dāng)音頻共振傳感器兩臂浸人到含沙水體時(shí)，兩臂間的介質(zhì)質(zhì)量發(fā)生變化，導(dǎo)致共振頻率發(fā)生變化，于是就可通過(guò)共振頻率的變化來(lái)間接實(shí)現(xiàn)含沙量的監(jiān)測(cè)。通過(guò)改變振子質(zhì)量m可改變音叉的共振頻率。設(shè)m0為音頻共振傳感器振子的等效質(zhì)量，mx為水體中含沙量，k為比例因子，則有

由諧振頻率f'（t）便可求出含沙量Mx。

實(shí)際應(yīng)用時(shí)，要用傳感器在清水中的固有頻率f0來(lái)標(biāo)定比例因子k。實(shí)測(cè)時(shí)mx不為零，于是通過(guò)諧振頻率f（t）就可以建立含沙量監(jiān)測(cè)的反演模型。

小波分析是將信號(hào)分解為不同頻率的信號(hào)，然后對(duì)這些信號(hào)進(jìn)行分析。本文將音頻共振傳感器一段時(shí)間內(nèi)獲取的含沙量數(shù)據(jù)組成貫序式數(shù)據(jù)塊進(jìn)行分析，實(shí)現(xiàn)每g個(gè)數(shù)據(jù)組成一個(gè)數(shù)據(jù)塊。任意含沙量時(shí)間序列數(shù)據(jù)塊信號(hào)f （t）的離散小波變換（discretewavelet transform，簡(jiǎn)稱DWT）系數(shù)為

WTf（j，k）==∫Rf（t）ψj，k（t）dt（2）式中：WTf（j，k ）為離散小波變換系數(shù)；j、k為小波基數(shù)；f為信號(hào)的頻率；ψj，k為小波變換核函數(shù)；R為實(shí)數(shù)域；q's.k（ t）為小波變換核函數(shù)的歸軛函數(shù)。

利用式（2）可以對(duì)f（t）進(jìn)行重構(gòu)，也就是小波變換的逆變換，其公式為式中：Cψ為對(duì)ψ（t）（小波變換核函數(shù)的實(shí)數(shù)域函數(shù)）提出的允許性條件；a、b分別為小波尺度；ψ（·）為小小波的位移與尺度伸縮。

2.2 卡爾曼和溫度融合模型

應(yīng)用動(dòng)態(tài)多尺度系統(tǒng)理論建立基于Kalman濾波的多尺度狀態(tài)估計(jì)方程，能有效降低測(cè)量的均方誤差。但是，一般的Kalman濾波模型僅對(duì)狀態(tài)方程和測(cè)量方程進(jìn)行融合處理，并沒(méi)有考慮環(huán)境信息對(duì)Kalman濾波模型的融合作用。當(dāng)不考慮控制信號(hào)的作用時(shí)，Kalman濾波的估值效果往往不夠理想。此時(shí)，傳感器需要進(jìn)行環(huán)境信息的融合處理，以削弱環(huán)境因素的影響。一般來(lái)說(shuō)，溫度對(duì)傳感器測(cè)量過(guò)程的影響較大，需要進(jìn)行溫度補(bǔ)償。傳統(tǒng)的溫度補(bǔ)償方法是通過(guò)尋找試驗(yàn)數(shù)據(jù)的統(tǒng)計(jì)規(guī)律來(lái)進(jìn)行軟件補(bǔ)償。本文建立的wavelet-Kalman模型是在進(jìn)行Kalman估值時(shí)直接進(jìn)行溫度補(bǔ)償。當(dāng)把溫度信息作為控制信號(hào)時(shí)，Kalman濾波公式可變?yōu)?/p>

X（j，k+1）=F（j，k）WTf（J，k）+G（j，k）t（j，k）+

w（j，k）（k≥0）（4）

Z（j，k）=H（j，k）WTf（j，k）（k≥0，j=1，2，…，N）

（5）式中：X（j，k+1）為在某一尺度j上k+1時(shí)刻的系統(tǒng)狀態(tài)；WTf（J，k）∈Rn×1；F（j，k）為系統(tǒng)矩陣，F(xiàn)（j，k）∈Rn×n；G（j，k）為溫度控制參數(shù)；t（j，k）為j尺度下環(huán)境溫度控制變量；Z（j，k）為j尺度下k時(shí)刻的值；H（j，k）為小波變換的參數(shù)矩陣；w（j，k）為噪聲，且滿足式（6）。

E{w（j，k）}=0（k≥0）

E{w（j，k）w（j，l）T}=Q（j，k）δkl（l≥0）（6）式中；E（·）為數(shù)學(xué)期望值函數(shù)；Q（j，K）為非負(fù)定對(duì)稱陣；δkl為狄拉克函數(shù)。

卡爾曼濾波通過(guò)無(wú)偏估計(jì)來(lái)達(dá)到方差最小，實(shí)現(xiàn)卡爾曼最優(yōu)估計(jì)。整個(gè)計(jì)算流程分為兩部分，一部分是估值的計(jì)算流程，另一部分是估值增益的遞推計(jì)算流程。具體推導(dǎo)公式可參照卡爾曼最優(yōu)濾波的5條基本公式。

3 貫序式wavelet-Kalman多尺度融合分析

3.1 基于wavelet的單尺度分析

先用小波函數(shù)對(duì)音頻共振傳感器的輸出信號(hào)進(jìn)行單尺度分解，設(shè)、為音頻共振傳感器輸出的原始信號(hào)，cal為小波分解后的低頻系數(shù)，cdl為小波分解后的高頻系數(shù)?？梢钥闯鲇脝纬叨鹊皖l系數(shù)cat重構(gòu)的含沙量信息與傳感器輸出的原始信號(hào)波形基本上一致，且低頻系數(shù)cat和高頻系數(shù)cdl的信號(hào)幅值都是原始信號(hào)、的一半。

3.2 基于wavelet的多尺度分析

由多尺度分析結(jié)果可以看出，尺度3（小波分析的第3層）的低頻系數(shù)ca3占整個(gè)數(shù)據(jù)長(zhǎng)度的1/4，尺度2的低頻系數(shù)cat的長(zhǎng)度占整個(gè)數(shù)據(jù)長(zhǎng)度的1/4。剩下長(zhǎng)度的一半是尺度1的高頻系數(shù)cdl。由此可知，低頻系數(shù)反映了含沙量信號(hào)的輪廓和基本信息，并且各層的系數(shù)數(shù)量隨著分解層數(shù)的增大而減小，原始含沙量數(shù)據(jù)相應(yīng)的噪聲得到抑制。因此，可通過(guò)對(duì)數(shù)據(jù)進(jìn)行多尺度變換，提取小波低頻系數(shù)來(lái)表示含沙量的整體特征，同時(shí)利用最小二乘法進(jìn)行回歸分析，建立不同分解層的反演模型并進(jìn)行檢驗(yàn)，確定最佳的小波分解估算模型。

3.3 基于貫序式wavelet-Kalman的多尺度融合分析

小波的分解和重構(gòu)過(guò)程中，若小波分解層數(shù)偏小，則低頻子序列仍是含沙量多個(gè)時(shí)間尺度上信號(hào)的混疊，不能很好地消除噪聲；反之，當(dāng)小波分解層數(shù)過(guò)大時(shí)，會(huì)產(chǎn)生虛假含沙量信息，造成含沙量測(cè)量誤差累積等問(wèn)題。目前，小波分析時(shí)確定分解層數(shù)的常用方法是通過(guò)自相關(guān)分析對(duì)各層小波系數(shù)進(jìn)行白化檢驗(yàn)，但這種白化自相關(guān)分析方法存在分解層數(shù)選擇不合理的問(wèn)題。從表1可以看出，4種小波母函數(shù)中，Db8（Daubechies8）、Sym8（Symlets）和Bior6.8（bior-wavf）均是以小波分解的第5層系數(shù)構(gòu)建的模型的均方根誤差扭MSE）最?。?.717 4 _< RMSE <7.6654），穩(wěn)定性最好，且Sym8小波模型的均方根誤差小于Db8、Coifs（Coiflets）和Bior6.8的。

4 結(jié)果與討論

4.1 貫序式wavelet-Kalman多尺度融合結(jié)果

分別采用Db8、Sym8，Bior6.8、Coifs共4種小波進(jìn)行多層離散分解后得到小波低頻系數(shù)，并對(duì)另外一組數(shù)據(jù)進(jìn)行反演，選取均方根誤差作為指標(biāo)來(lái)衡量模型精度。RMSE越大反演模型的精度越低。

貫序式wavelet-Kalman多尺度融合模型選取Sym8小波的ca5層進(jìn)行小波分解。由于貫序式wavelet-Kalman多尺度融合模型考慮了環(huán)境因素的影響，因此輸出結(jié)果更接近實(shí)測(cè)值，反演效果更好，誤差更小，見(jiàn)圖2、圖3。

4.2 貫序式wavelet-Kalman融合誤差比較

為了衡量貫序式wavelet-Kalman多尺度融合模型的效果，給出了多個(gè)反演模型的均方根誤差、最大絕對(duì)誤差（MaxAE）、平均絕對(duì)誤差（MeanAE）和平均相對(duì)誤差（MeanRE），見(jiàn)表2。由表2可以看出，貫序式wavelet-Kalman多尺度融合模型的效果最好，其MeanAE為3.95kg/m3，RUSE為3.13kg/m3，MaxAE為14.15kg/m3，MeanRE為1.05kg/m3。表2中的wavelet-Kalman濾波模型是指用Kalman算法進(jìn)行濾波處理，沒(méi)有用Kalman算法進(jìn)行溫度融合處理，其處理效果較貫序式wavelet-Kalman多尺度融合模型的稍差一些，而Kalman和溫度融合模型的誤差比一元和多元線性回歸模型的小。

從表2還可以看出，經(jīng)過(guò)wavelet多尺度分析和Kalman融合后，音頻共振傳感器的測(cè)量誤差都比多元及一元線性回歸的小。整體而言，一元線性回歸的相對(duì)誤差較大，但平均相對(duì)誤差也在5%以內(nèi)，可滿足含沙量測(cè)量的需要。

一元線性回歸方程是采用前46組數(shù)據(jù)塊擬合得出的，其表達(dá)式為：

y=2.774x+356.3（7）式中：J為含沙量；x為加沙過(guò)程中含沙量信息通過(guò)音頻共振傳感器轉(zhuǎn)換的電流值。

多元線性回歸方程為

S=-477.9815+1.4152X1-0.2769X2（8）式中：X1為音頻共振傳感器的輸出信號(hào)；X2為溫度傳感器的輸出信息值。

一元線性回歸模型的MeanAE為16.22kg/m3，RUSE為20.18kg/m3，MaxAE為46.90kg/m3。多元融合后的誤差較一元線性回歸的小些，說(shuō)明考慮不同的環(huán)境因素影響，融合效果有提高，但還不是很理想。

5結(jié)語(yǔ)

（1）當(dāng)應(yīng)用多尺度分析方法提取含沙量信息的突變值并能融合環(huán)境因素時(shí)，監(jiān)測(cè)誤差相對(duì)小，監(jiān)測(cè)精度和穩(wěn)定性相對(duì)高。

（2）貫序式wavelet-Kalman多尺度融合模型能有效提取含沙量的突變信息，同時(shí)考慮了環(huán)境因素影響，提高了含沙量監(jiān)測(cè)的精度，增強(qiáng)了測(cè)量系統(tǒng)的穩(wěn)定性。

（3）音頻共振法在高懸浮含沙量監(jiān)測(cè)時(shí)誤差較小，監(jiān)測(cè)精度較高，表明音頻共振傳感器較適合于高含沙量河流中使用，特別適用于黃河流域的高懸浮含沙量監(jiān)測(cè)。