王倩 吳江雪 侯慶文 陳先中

在高爐煉鐵生產(chǎn)過程中,料面信息對安全生產(chǎn)、節(jié)能減排等目標(biāo)有著重要的指導(dǎo)作用,提取料面信息是控制爐況的關(guān)鍵.高爐內(nèi)部是高溫高壓、密閉的、多粉塵的惡劣環(huán)境,同時高爐料面是具有非均勻的流態(tài)化特征的礦物–煤氣–焦炭混雜的多元高溫熔固體共存的粗糙表面.在此基礎(chǔ)上,隨著雷達(dá)探測全壽命周期的晚期衰減,料面圖像特征會出現(xiàn)嚴(yán)重的模糊、缺失與淹沒,因此高爐料面研究成為難點與熱點.采用工業(yè)雷達(dá)調(diào)頻連續(xù)波(Frequency modulated continuous wave,FMCW)方法可以測量高爐內(nèi)部料面實時生產(chǎn)信息.通過可視化操作平臺輔助現(xiàn)場人員及時觀測和調(diào)整布料策略.

為了獲取有效的料面信息,近年來,陳先中等利用帶權(quán)值的雙算法綜合料面生成模型,采用最小二乘法和NURBS(Non-uniform rational B-splines)曲面全局插值算法進行曲面擬合,建立了基于圖像處理的在線檢測模型[1].Zankl 等基于單通道FMCW 信號模型,結(jié)合數(shù)字波束形成方法和進一步的先進信號處理算法,推導(dǎo)出一種多通道信號模型,采用CLEAN聚類算法[2].文獻[3]利用多回路充電結(jié)束時的立方曲線方程估計負(fù)荷剖面,解決了基于真實多雷達(dá)數(shù)據(jù)的多回路充電負(fù)荷分布實時估計問題[3].文獻[4]提出了改進的卡爾曼濾波和異常檢測模型,以協(xié)同跟蹤高爐爐料深度[4].而高爐內(nèi)部呈現(xiàn)高溫燃燒和強氣流沖擊,受到運動溜槽和固定十字測溫的多重物理干擾、微波近場天線在腐蝕和黏附等電磁干擾以及料面本身的瑞利雜波等問題,使得FMCW 雷達(dá)測得的料面回波信號信噪比低,狀態(tài)不穩(wěn)定[5].當(dāng)料面呈流態(tài)化甚至噴涌態(tài)時,回波會出現(xiàn)料線測點劇烈跳動現(xiàn)象,考慮Hawking-Penrose奇點定理條件[6],如果存在強能量條件時,則存在封閉陷獲面和緊致無邊非時序點集,依靠傳統(tǒng)尋峰法會導(dǎo)出錯誤料面形狀,嚴(yán)重誤導(dǎo)高爐布料操作.但是以上方法,均采用傳統(tǒng)信號處理方式處理料面回波,忽視料面流態(tài)化特征和顆粒物飛濺等問題,不適用于信噪比低和不穩(wěn)定的信號.桂衛(wèi)華院士的團隊提出了平行低光損背光高溫內(nèi)窺鏡的檢測方法來獲取高爐全料面圖像[7].這種設(shè)計實現(xiàn)了近距離取像,受粉塵干擾小,成像具有高分辨率和高清晰光學(xué)圖像特征,但定視角的爐料表面圖像不包含縱深信息,基于微波原理的料面點云圖像包含距離信息和料面流態(tài)化信息,可以同時實現(xiàn)料面形狀和氣流狀態(tài)的觀測.

以某鋼3號高爐為研究背景,借鑒目前遙感領(lǐng)域高精度地面SAR 成像技術(shù)原理,通過連續(xù)擺動雷達(dá)獲得料面徑向高密度連續(xù)點云回波信號,模擬SAR 成像方式.

在工業(yè)生產(chǎn)中,高爐爐頂壓力可通過模糊解耦得到精準(zhǔn)控制[8],電磁波發(fā)射信號在爐頂高溫高壓和粉塵煙霧環(huán)境下產(chǎn)生強烈電磁噪聲和干擾信號,但經(jīng)過特殊頻段和降噪處理,還是可以檢測到料面形狀,一般檢測范圍20 m,滿足工業(yè)SAR 雷達(dá)成像的范圍要求.本文將回波信號轉(zhuǎn)化成圖像進行料面有效特征提取.在高爐冶金過程中,物料狀態(tài)多樣,爐內(nèi)顆粒物分布復(fù)雜,呈現(xiàn)流態(tài)化特征.本文采用點云表現(xiàn)形式,體現(xiàn)了料面參數(shù)的概率分布值,反映了物料的反射特性和料面的流態(tài)化特征.根據(jù)高爐現(xiàn)場條件,物料顆粒物能量的高低變化復(fù)雜,故借助能量重心原則提取料面能量點.在雷達(dá)回波模型基礎(chǔ)上,對料面能量點云的距離頻率進行估計,引入克拉美羅下界(Cramer-Rao lower bound,CRLB)作為最優(yōu)估計性能.

本文對一維單點信號頻譜圖擴展到引入料面半徑方向空間信息的二維連續(xù)多點信號頻譜圖.將傳統(tǒng)信號處理方法與圖像處理方法結(jié)合,設(shè)計多種濾波方法處理高斯噪聲等干擾信號.以聚焦雷達(dá)回波高能量點,兼顧低能量點為原則,銳化料面峰脊,提出基于能量重心法的提取方法,進而提取料面主能量點生成3D點云料面模型并提取料線距離.通過料面信息的概率分布推導(dǎo)頻率估計的CRLB,同時對比Single-Tone算法,對現(xiàn)場數(shù)據(jù)的距離頻率估計誤差更接近CRLB,對實測高爐雷達(dá)回波數(shù)據(jù)進行成像驗證,在料面處于流態(tài)化和噴涌態(tài)時,本文所提方法能有效處理其產(chǎn)生的低信噪比回波信號,準(zhǔn)確提取料面信息.

1 雷達(dá)信號數(shù)據(jù)分析與預(yù)處理

1.1 雙擺動雷達(dá)探測模型

本文使用雙擺動FMCW 雷達(dá)對高爐內(nèi)部料面情況進行探測,圖1為單個雷達(dá)安裝示意圖.

圖1 雷達(dá)安裝及掃射范圍示意圖Fig.1 Radar installation and distance diagram

其中,雷達(dá)測距分辨率K為0.0915,平均每次擺動掃描收發(fā)50個回波信號[9].

1.2 粗糙料面回波信號分析

FMCW雷達(dá)高爐料面信號干擾主要分為隨機干擾與固定干擾.隨機干擾包含近場天線強物理干擾、瑞利雜波和焦炭粉塵等諸多干擾;固定干擾主要來自高爐內(nèi)部的十字測溫裝置、周期性運動的溜槽以及邊緣測量點的爐壁反射.除此之外,由于被測料面是非均勻流態(tài)化的粗糙表面,其粗糙度導(dǎo)致雷達(dá)回波能量分散衰弱.

實際生產(chǎn)中雷達(dá)所測兩種典型回波信號的頻譜圖如圖2所示.圖2(a)為理想信號,圖2(b)為典型低信噪比信號,此時干擾較大,復(fù)雜的散射特性使得傳統(tǒng)信號處理方法難以區(qū)分真實的料面回波與干擾信號.其中FFT的幅值為雷達(dá)回波能量值,即雷達(dá)散射截面(Radar cross section,RCS)的大小.

本文設(shè)計的雙擺動FMCW雷達(dá),將空間信息引入料面回波信號,經(jīng)過沿高爐料面半徑方向一次擺動掃描,即可獲取一組包含料面形狀與距離信息的雷達(dá)回波信號.由于料面距離有一定范圍,對原始時域信號后FFT變換后,前120個譜線號即包含料面所有信息,將其轉(zhuǎn)換成標(biāo)準(zhǔn)64階圖像進行分析.如圖3所示,可以看出高爐料面形狀是一條波動的具有典型特征的過渡帶狀曲線,與干擾噪聲有著明顯的視覺差異.

同時料面回波顯示,高爐料面是在一定范圍內(nèi)波動的模糊層.這是由于高爐料面具有典型的非穩(wěn)態(tài)粗糙流態(tài)界面特性,以及在上升煤氣強氣流和爐料滑動等影響下,會出現(xiàn)波動,這使傳統(tǒng)尋峰法難以判斷料面距離精確數(shù)值.為了更加準(zhǔn)確還原料面實際形狀與真實距離信息,需引入雷達(dá)擺動掃描時各點回波信號對應(yīng)的角度信息,將矩陣圖像A轉(zhuǎn)換成具有實際空間意義的扇形空間下料面回波信號2D原始頻譜圖,如圖4所示.

1.3 雷達(dá)信號去噪預(yù)處理

依據(jù)遙感SAR 旋轉(zhuǎn)雷達(dá)成像原理,把傳統(tǒng)的10～20點固定位置數(shù)據(jù)抽取成像,加密采樣成多個料面位置的點云數(shù)據(jù).

單個擺動周期內(nèi)連續(xù)采樣N個料面位置信息,經(jīng)FFT變換后,截取前120個譜線數(shù)據(jù),生成120×N的二維圖像矩陣A:

圖2 原始雷達(dá)信號頻譜對比Fig.2 Comparison of original radar signal spectrum

圖3 料面徑向回波信號2D頻譜圖對比Fig.3 Comparison of 2D spectra of radial echo signal of blast surface

圖4 扇形空間下料面回波信號2D原始頻譜圖對比Fig.4 Comparison of 2D original spectra of fan-shaped space blanking surface echo signal

其中,X1,X2,···,X N表示料面徑向上的采樣點的頻譜數(shù)據(jù).

由于高爐內(nèi)部環(huán)境惡劣,非穩(wěn)態(tài)粗糙流態(tài)界面的距離信息因為界面的粗糙、非穩(wěn)態(tài)特性,呈現(xiàn)出距離信息寬帶化;同時,又因為各種干擾的存在,在距離寬帶處又表現(xiàn)出距離信息與干擾信息共存的特征.針對此類非穩(wěn)態(tài)粗糙流態(tài)界面距離信息的提取,圖像矩陣A中包含大量干擾,因此采用多級濾波方法對其進行去噪處理.

1)中值濾波

中值濾波是一種基于統(tǒng)計的非線性濾波方法,顧及到一定鄰域內(nèi)噪聲點的異常值.其中選擇120點處理數(shù)據(jù),選擇合理的階次平滑數(shù)據(jù).

2)窗函數(shù)濾波

窗函數(shù)作用滿足式(2)[10]:

其中,h(n)是實際的濾波器的單位采樣響應(yīng),hd(n)是理想濾波器的單位響應(yīng),ω(n)為Kaiser窗函數(shù),滿足[11]:

由于料面信息集中在50～100 Hz頻段,故設(shè)定截止頻率以滿足帶通濾波的要求.其中過渡的頻譜段長3 Hz,阻帶截止頻譜號為40 Hz,則通帶的起始頻譜號是43 Hz;通帶截止頻譜號是90 Hz,阻帶起始頻譜號是93 Hz.

3)基于方差的帶寬提取法

根據(jù)單次掃描料面回波圖中反映的料面與干擾明顯的形狀特征區(qū)別,料面所在信號頻帶方差較大,因此本文引入基于方差的帶寬提取法濾除低頻干擾信號.

利用對二維矩陣A的每一行數(shù)據(jù),求方差排序,選取方差最大的前40行,濾除無用信息點,以達(dá)到分離頻帶的目的.

其中X i是頻譜號為i的觀測值,var(X i)為方差計算,V是閾值。

求取每個頻譜帶的方差均值,設(shè)經(jīng)過上步后寬帶數(shù)目為k,寬帶的頻譜號數(shù)目是m(i),i=1,2,···,k,則方差均值可由向量K V AR(i)表示:

只保留方差均值最大的頻譜寬帶,如式(6)所示,濾除較強的干擾.

4)迭代閾值法

本文引入基于圖像分割的迭代閾值法[12].獲取迭代計算的自適應(yīng)閾值,將圖像分為前景和背景兩個區(qū)域[13],進行圖像的閾值濾波,去除隨機干擾,得到干凈清晰的料面距離帶寬.迭代閾值算法步驟如下.

步驟1.初始化:選取頻域最大值fmax,為頻域最小值fmin,初始閾值:

步驟2.根據(jù)將距離寬帶的頻域信息分割成兩個區(qū)域R1和R2,求區(qū)域均值:

其中,R1中包含N1個元素,R2中包含N2個元素.xi為區(qū)域R1中的第i個元素,yi Q為區(qū)域R2中的第i個元素.

步驟3.計算新閾值:

步驟4.求取閾值差:

步驟5.若?T

步驟6.輸出最佳閾值.

步驟7.結(jié)束.

2 基于CRLB的料面點云提取模型

2.1 克拉美羅下界分析

針對雷達(dá)回波信號,分析其克拉美羅下界(Cramer-Rao lower bound,CRLB),可計算無偏估計中能獲取的最佳估計精度,更可評估不同參數(shù)估計方法的性能[14].

由高爐料面特性可知,料面圖像中的噪聲主要來自高斯白噪聲和泊松噪聲,故采用高斯函數(shù)作為點擴散函數(shù)(Point spread function,PSF).高爐料面距離信息是料面成像的重要參數(shù).根據(jù)FMCW雷達(dá)測距原理,距離R與頻率f0關(guān)系表達(dá)式為:

其中,B為雷達(dá)調(diào)制帶寬,f s為采樣頻率,c為光速.因此對料面距離的估計也就是對回波信號的頻率估計.通過推導(dǎo)頻率估計的CRLB可衡量本文方法對料面距離的估計精度[15].雷達(dá)料面回波信號模型為:

計算其對數(shù)似然函數(shù)可以得到:

對對數(shù)似然函數(shù)求參數(shù)二階導(dǎo)數(shù)并求期望得到費雪信息(Fisher)[16]:

針對該信號模型,料面距離頻率估計的CRLB為[17]:

2.2 峰脊銳化

為了更加科學(xué)估計料面距離,根據(jù)邊緣能量聚焦和峰脊銳化等數(shù)理統(tǒng)計和圖像處理方法利用加權(quán)采樣法篩選與料面相關(guān)度更高的標(biāo)志點.

將二維單次掃描雷達(dá)料面回波信號頻譜圖增加一維,引入笛卡爾坐標(biāo)系的z軸,代表料點的幅值大小,反映料線峰脊的高低,故將幅值作為峰脊銳化、識別料面的重要特征.圖5多視角展示了理想信號(a)與低信噪比信號(b)在3維空間下峰脊分布情況.其物理意義是雷達(dá)回波的強度,間接代表其料點的數(shù)據(jù)質(zhì)量,故將提取的數(shù)據(jù)點分為重要標(biāo)志點和次要標(biāo)志點兩大類[19].z軸幅值較高的數(shù)據(jù)點亮度較高,即雷達(dá)回波強度較強,屬于重要標(biāo)志點(高能量點);反之,則為次能量點.

本文根據(jù)數(shù)理統(tǒng)計方法,采用加權(quán)采樣的方式統(tǒng)計標(biāo)志點的幅值,從以下兩方面達(dá)到峰脊銳化的目的:1)聚焦重要標(biāo)志點,應(yīng)用小步長采樣,達(dá)到高能量點密集點云的效果,反應(yīng)料面流態(tài)化情況;2)兼顧次能量點,采用較大的采樣步長,濾除邊緣雜波點,達(dá)到次能量點稀疏點云的效果.具體算法步驟如下:

步驟1.采樣點線性分區(qū):根據(jù)采樣點線性劃分N個子區(qū)域:V1,V2,···,VN,識別出子區(qū)域中的料面位置信息并進行空間對齊.

步驟2.逐區(qū)間求矩陣峰值:統(tǒng)計子區(qū)域內(nèi)矩陣峰值,記為R=

步驟3.峰值排序.

步驟4.確定采樣步長:根據(jù)排序后的矩陣R賦予子區(qū)域權(quán)重,決策其采樣步長.

步驟5.雙向采樣非零數(shù)據(jù)點:同時從y軸和x軸雙向采樣以聚焦邊緣能量.

2.3 基于能量重心的料線提取

圖5 三維空間料面回波信號強度峰脊分布Fig.5 Peak ridge distribution of echo signal strength in three dimensional space

在料面上部存在顆粒物噴涌,且受到高斯白噪聲和散粒噪聲的影響下,即使經(jīng)過濾波與銳化峰脊處理,代表料面實際位置的回波信號附近也彌散著部分干擾回波信號,但這些干擾回波信號能量較弱.同時,料面的流態(tài)化特性,以料面回波信號能量為一條具有寬度、能量變化的波動曲帶的形式反映在2D料面回波信號頻譜圖上.因此,為了能提取出料線信息,本文根據(jù)能量重心法(The energy centrobaric correction)對峰脊銳化后料面點云的進行料線提取,濾除多余點云數(shù)據(jù)[20].

能量重心法適用于頻譜校正,n=3時效果佳,即通過三點卷積校正幅值的方法校正頻率[21].在高斯白噪聲背景下,功率譜滿足以下分布:

高斯白噪聲影響下,能量重心校正的歸一化頻率滿足:

其中,kr:功率譜最大值譜線;Pk r+i:第kr條功率譜線.

故主瓣中心坐標(biāo)k0(k1為最大值譜線):

通過能量重心法對單點料面點云數(shù)據(jù)進行主能量點提取后,對其進行傅里葉擬合,提取出一條波動光滑料線.根據(jù)能量重心理論剔除野值后,點云數(shù)據(jù)需擬合得到理想的光滑料線.采用5階傅里葉級數(shù)的線性傅里葉擬合(Fourier linear combiner,FLC)進行曲線重構(gòu)[22].公式:

其中,an是正弦項系數(shù),bn是余弦項系數(shù),C是常數(shù),T是周期.

綜上所述,將高爐料面點云銳化成像算法可表示為圖6所示流程圖.

圖6 高爐料面點云銳化成像算法流程圖Fig.6 Algorithm flow chart of burden surface point cloud sharpen image processing

3 高爐料面點云成像仿真驗證

本文以某鋼3#高爐2號雷達(dá)2018年4月實測回波信號為實例,驗證分析本文所設(shè)計的基于CRLB的料面點云銳化成像算法性能.為了證明新方法能夠有效處理低信噪比且不穩(wěn)定信號,把理想信號與非理想信號處理結(jié)果進行對比.

3.1 多級濾波效果驗證

本文設(shè)計了三種不同的多級濾波模型,MWATF法(Mid-value and window and artificial threshold filtering)、VITF法(Variance and iterative threshold filtering)、MWITF法(Mid-value window and iterative threshold filtering).每種濾波器由不同濾波方式串聯(lián)組成,具體組成方式見表1.

表1 三種濾波模型Table 1 Three filtering models

通過處理多組低信噪比料面回波實驗數(shù)據(jù)進行對比,MWITF法效果最好.MWATF法由于人工設(shè)定閾值,適用普遍性差,濾波后圖中仍存在明顯的干擾帶;而VITF法與MWITF法相比,存在濾波過渡的問題,MWITF保留了更豐富的料點,使圖像更真實;因此選擇MWITF法為本成像系統(tǒng)濾波方法.

對于料面回波信號,三種濾波效果如圖7.

圖7 三種濾波模型效果Fig.7 Three filtering model effects

某鋼雷達(dá)的現(xiàn)場安裝圖如圖8所示,在高溫高壓黏附和強粉塵環(huán)境下,連續(xù)運行3個月獲得的清晰料面電磁散射圖,粉塵幾乎被完全壓制,料線特征十分明顯,如圖9所示.

圖8 高爐雷達(dá)的現(xiàn)場安裝圖Fig.8 Field installation of the blast furnace radar

3.2 距離頻率估計CRLB驗證

通過加權(quán)采樣法將峰脊銳化,提取料面標(biāo)志點數(shù)據(jù).為了去除部分野值點影響,采用能量重心法提取能量點,使目標(biāo)點盡可能逼近主能量點的情況下,又考慮到分散在料面周圍的次能量點能量幅值大小的影響,根據(jù)最終得到的能量點算出距離,提取料線.效果如圖10.

圖9 高爐料面電磁散射圖Fig.9 Electromagnetic scattering image of burden surface

圖10 能量點提取與料線擬合(Data1為峰脊銳化后的標(biāo)志點,Data2為提取的能量點)Fig.10 Energy point extraction and material line fitting(Data1 is the mark point after peak ridge sharpening,and Data2 is the energy point extracted)

為了驗證本文算法對實測信號頻率估計的有效性,選取Single-Tone算法進行比較,并與CRLB限進行對比.Single-Tone 算法是基于DFT從極大似然估計量MLE和加權(quán)平均估計量WAE 入手的頻率估計方法[23].選取?5 d B到5 d B的高爐雷達(dá)料面回波信號,在每一個信噪比下用本文方法與Single-Tone方法各對隨機選取的100個實測信號進行相對均方跟誤差計算.圖11為不同信噪比下兩個算法頻率估計相對均方誤差(RMSE)比較.由圖11可知,從估計精度而言,本文估計方法更接近CRLB,估計準(zhǔn)確度較好.從算法穩(wěn)定性分析,本文算法頻率估計RMSE隨著信噪比增加逐步減小,逼近CRLB,而Single-Tone算法波動較大,因此本文算法提取料線性能更加真實可靠,運行更加穩(wěn)定,魯棒性強,抗噪性好.

圖11 不同信噪比下各算法頻率估計性能Fig.11 The performance of each algorithm under different SNR is estimated

3.3 3D料面點云成像驗證

在得到離散的料面能量點云數(shù)據(jù)與料線后,假設(shè)高爐料面是對稱的曲面,將點云數(shù)據(jù)與料線進行旋轉(zhuǎn),重構(gòu)出3D料面點云圖像與光滑料面曲面圖像,高爐加焦的料面模型如圖12(a)所示,平臺型料面如圖12(b)所示.隨機選取某鋼3號高爐2號雷達(dá)2018年4月和8月1 000組掃描信號,與原信號頻譜圖對比,符合預(yù)期料面特征的信號816組,料面提取準(zhǔn)確率達(dá)到81.6%.相比傳統(tǒng)信號處理方法,基本達(dá)到工業(yè)實用的要求.

4 結(jié)論

結(jié)合SAR 成像原理,本文從圖像處理角度,對實際料面回波頻譜圖像進行多級濾波處理,濾除低頻段強干擾;通過峰脊銳化與能量重心法提取離散的代表料面特征的能量點數(shù)據(jù),根據(jù)這些點云數(shù)據(jù)進行3D料面點云曲面與光滑料面重構(gòu),同時基于CRLB分析本文算法對實測信號的頻率估計精度,為高爐雷達(dá)測距準(zhǔn)確度提供了科學(xué)參考,驗證了本文算法頻率估計精確度高.并且相較于傳統(tǒng)單點采樣擬合料面法,本文保留了大量的料面形狀細(xì)節(jié),同時直觀反映了料面流態(tài)化特征,得到的料面模型更加貼近實際生產(chǎn)中的高爐料面,料面提取準(zhǔn)確率達(dá)到81.6%.工業(yè)雷達(dá)在高爐環(huán)境下,其天線的探測全壽命約1年,達(dá)到壽命晚期的時候,天線污染嚴(yán)重,出現(xiàn)大量低信噪比信號,此時圖像特征的模糊、缺失與淹沒現(xiàn)象加劇,目前更加智能化的圖像提取方法和信號處理算法還不完善,故該類算法是今后研究工作的重點,有待繼續(xù)研究.

圖12 3D料面點云成像效果Fig.12 Imaging effect of 3D surface point cloud