楊江洪 顧宗山 盧志忠
(1.中國電子科技集團(tuán)公司第三十八研究所 合肥 230088;2.哈爾濱工程大學(xué) 哈爾濱 150001)
海浪是海洋動(dòng)力環(huán)境的重要因素,波高是衡量海水運(yùn)動(dòng)及海浪能量一個(gè)重要參量。實(shí)時(shí)提供精確的海浪監(jiān)測信息有助于人們的海洋作業(yè)安全,因此海浪監(jiān)測具有十分重要的科研意義和實(shí)用價(jià)值[1]。當(dāng)前海浪監(jiān)測方式有:人工觀測,依賴于觀測者經(jīng)驗(yàn)積累,準(zhǔn)確性較差;浮標(biāo)、流速壓力傳感器能直接獲得海浪參數(shù),測量精度高,但觀測范圍較??;聲學(xué)傳感器測量距離遠(yuǎn),然而自身噪聲較高,難以獲得高精度參數(shù)。以上的測量方法雖然可以直接獲得海浪參數(shù),但面臨著精度不夠高、花費(fèi)較大,觀測面小等問題。伴隨著雷達(dá)技術(shù)技術(shù)飛速發(fā)展,雷達(dá)遙感裝備觀測海浪成為當(dāng)前海洋測量技術(shù)的發(fā)展趨勢。利用導(dǎo)航雷達(dá)對海況進(jìn)行實(shí)時(shí)監(jiān)測兼顧導(dǎo)航和測浪的功能,最大化提升導(dǎo)航雷達(dá)用途。
目前,遙感雷達(dá)圖像反演海浪波高有兩種方法,分別是海浪譜反演法和陰影圖像反演法。海浪譜反演法是利用傅里葉變換獲得雷達(dá)圖像波數(shù)能量譜,計(jì)算信噪比,通過信噪比的平方根成線性關(guān)系來估測有效波高[2-5]。而海浪譜反演法難度較大,圖像譜中利用色散濾波器提取海浪波數(shù)能量譜,估流精度直接影響色散濾波器性能,實(shí)際中準(zhǔn)確估流較為困難[5],另外調(diào)制傳遞函數(shù)與雷達(dá)機(jī)制有關(guān),不同雷達(dá)調(diào)制傳遞函數(shù)不同,這對于推廣應(yīng)用帶來難度,操作難度提高通用性降低。而陰影圖像反演法避開這一復(fù)雜處理過程,它是基于海浪在雷達(dá)回波中各種物理調(diào)制特性映射到了海雜波散射模型,具備幾何光學(xué)粗糙面某些特性,通過圖像處理提取海浪參數(shù)[16]。1990年, Wetzel[6]利用理論分析海面假設(shè)為高斯面,分析低掠射角入射下的海面后向散射模型,利用散射強(qiáng)度來估測波高,奠定圖像法反演波高的理論基礎(chǔ)。1995年,Henschel[7]等人利用實(shí)測實(shí)驗(yàn)數(shù)據(jù)來驗(yàn)證Wetzel模型,并取得一定成果,然而由于缺乏大量數(shù)據(jù)與校驗(yàn),波高估測誤差較大。1998年,J.R Buckley等人[8]對Henschel估測方法的局限性原因進(jìn)一步分析,通過利用歸一化雷達(dá)目標(biāo)有效截面校驗(yàn)閾值,使得量測范圍變大、精度提高。2000年,R. Gangeskar[9]通過大量實(shí)驗(yàn)提取圖像灰度方差反演波高,相關(guān)系數(shù)為0.9,算法的精度得到很大提升。2001年,Buckley[10]通過仿真實(shí)驗(yàn)來揭示觀測值與Wetzel理論值間偏差原因,并利用不同位置上的回波強(qiáng)度來校驗(yàn)閾值,精確估計(jì)了不同海況下的有效波高;進(jìn)一步從理論上證實(shí)了閾值隨作用距離發(fā)生變化的。2012年 W. J. Plant[11]研究了極化方式對幾何陰影與局部陰影調(diào)制問題,實(shí)驗(yàn)結(jié)果表明局部陰影是海雜波陰影圖像中主要組成部分;W.J.Plant深度剖析了陰影成因的問題,為日后精確的檢測陰影提供支撐。2014年,Rune[12][16]通過分析W.J.Plant的實(shí)驗(yàn)發(fā)現(xiàn)工程應(yīng)用中,難以從圖像中辨別兩種陰影的影響的主次,通過大量數(shù)據(jù)分析,仿真試驗(yàn)得出且兩種陰影差異較小, Rune基于隨機(jī)粗糙面幾何陰影理論利用邊緣檢測提取圖像中陰影閾值,獲取有效波高,算法提升了估測精度,降低試驗(yàn)成本及簡化試驗(yàn)流程。2017年,盧志忠等[15-16]利用實(shí)測雷達(dá)圖像數(shù)據(jù)分析Rune結(jié)論,發(fā)現(xiàn)陰影閾值隨著海況變化。2017年,衛(wèi)延波博士[17]采用自適應(yīng)分塊以及淺水條件修正等改進(jìn)策略改進(jìn)幾何陰影法反演有效波高,提高了有效波高的估測精度。2019年,Giovanni Ludeno等[18]基于原始雷達(dá)圖像和相應(yīng)的非校準(zhǔn)波高圖像之間的相關(guān)性,利用比例因子改變波高圖像的幅度,建立雷達(dá)成像的數(shù)學(xué)模型。通過分析Rune的實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)陰影比例隨著海況而變化,而在掠射角小于2°區(qū)域的這種變化最明顯,對大量實(shí)驗(yàn)分析發(fā)現(xiàn)波高與陰影比例間存一定線性關(guān)系,還需要通過大量實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證。
本文在Rune的研究基礎(chǔ)上,結(jié)合圖像處理技術(shù)及統(tǒng)計(jì)學(xué)分析,研究不同波高下陰影比例的差異,用統(tǒng)計(jì)分析最小二乘擬合固定區(qū)域下陰影比例與波高間的線性模型,用擬合模型估測原始數(shù)據(jù)的波高,最終將估測的波高與參考波高作比較來檢驗(yàn)算法的可靠性與可行性。
海洋運(yùn)動(dòng)學(xué)上認(rèn)為海面是由多種隨機(jī)波動(dòng)疊加而成,在電磁散射模型中海面構(gòu)造為二維隨機(jī)粗糙面。將二維粗糙表面(例如海面、沙漠,戈壁等)簡化成一維粗糙表面模型[14],降低研究雜度,增加工程實(shí)用性。
低掠射下的海面由亮暗條紋組成,本文提出的波高反演策略基于粗糙面理論的延拓。如圖1所示為原始雷達(dá)圖像。
圖1 原始雷達(dá)圖像
估測陰影閾值,先提取圖像陰影邊緣。利用邊緣檢測找到陰影邊緣。
1.2.1 圖像陰影邊緣檢測
通過對雷達(dá)圖像單一方向上做像元卷積運(yùn)算式(1)獲取圖像梯度[12][16]。算子如式(2)依次對雷達(dá)圖像相鄰的8個(gè)方向上分別進(jìn)行卷積運(yùn)算。右側(cè)卷積為
IE1(r,θ)=I(r,θ)?H1(r,θ)
(1)
其中I(r,θ)是極坐標(biāo)表下原始雷達(dá)圖像,r是徑向斜距,θ是方位角,IE1(r,θ)是左側(cè)算子。像元差分算子給出
(2)
繪制原始圖像頻率曲線取前M%圖像的像素值作為閾值,遍歷原始雷達(dá)圖像,然后利用H1~i算子做卷積得到圖像梯度IE1~I(xiàn)E8。
以IE1中的第i幅邊緣圖像的閾值求取為例
(3)
N%為IE1映射到圖像中前N%的灰度值,之后每個(gè)方向的邊緣圖像疊加成一個(gè)合成的邊緣圖像IT,將各方向上邊緣圖像做和[12][16]
(4)
1.2.2 剔除孤立點(diǎn)及奇異值
對于孤立點(diǎn)去除利用邊緣比較法將每個(gè)點(diǎn)與閾值τF做比較
(5)
1.2.3 提取陰影邊緣區(qū)域
(6)
(7)
其中Ips為可疑陰影邊界,If為陰影邊界[12]。
通過對海雜波雷達(dá)圖像作個(gè)方向邊緣檢測得到If,繪制If的灰度直方圖。If映射到原始雷達(dá)圖像η中,并獲取統(tǒng)計(jì)直方圖FH(η),對直方圖取眾數(shù),即可得到陰影灰度閾值τS[12]為
τS=mode(FH(η))
(8)
將用閾值τS將原始圖像分割成陰影區(qū)與非陰影區(qū)并求取陰影比例Rs(γ)。
(9)
對圖像處理結(jié)果分析發(fā)現(xiàn),低掠射角下區(qū)域面積變化較大,通過觀察了大量圖像數(shù)據(jù)發(fā)現(xiàn)陰影比例大小與波高有關(guān)。因此選取掠射角小于10°,L×L個(gè)像素的區(qū)域求取陰影比例為
(10)
其中L為區(qū)域尺寸,SIs為所選區(qū)域內(nèi)陰影面積。
通過假設(shè)有效波高與陰影比例呈線性關(guān)系,目標(biāo)函數(shù)即為
Hs=A+B·Rs
(11)
最小二乘擬合有效波高關(guān)系式,目標(biāo)函數(shù)定義為
(12)
根據(jù)最小二乘法的計(jì)算原理[15],可以得到標(biāo)定系數(shù)A和B的估計(jì)值為
(13)
從2007年至今,哈爾濱工程大學(xué)在國家某裝備項(xiàng)目資助下積極開展航海雷達(dá)監(jiān)測海浪、海流、海面風(fēng)場、近岸水深、波面高度、溢油及海冰等方面的應(yīng)用研究,先后在北海、黃海、東海、南海等地開展科研試驗(yàn),積累了大量實(shí)測數(shù)據(jù)。實(shí)驗(yàn)航海雷達(dá)為X波段機(jī)掃導(dǎo)航雷達(dá)、天線為桿式端饋開槽波導(dǎo),天線轉(zhuǎn)速:26+2、-6/min,參數(shù)如表1所示,信號(hào)帶寬為20 MHz,分辨率為徑向7.5 m、角度向0.1°,雷達(dá)圖像回波強(qiáng)度是0~-2.5 v,A/D變換后映射為0~8191范圍內(nèi),圖像幀時(shí)間約2.5 s,岸基天線架設(shè)高度為45 m、船基為25 m。
表1 X波段導(dǎo)航雷達(dá)參數(shù)
本文選取數(shù)據(jù)源于2010~2013年期間東海海域?qū)崪y數(shù)據(jù)共210組,數(shù)據(jù)具有多種海況,遮擋小,圖像質(zhì)量好,具有一定代表性。所選岸基與海上數(shù)據(jù)的風(fēng)速在15 m/s左右,海上波高范圍為3~4 m,岸基為2 m左右。由于遠(yuǎn)區(qū)雜波競爭等因素,故本文選取的徑向圖像上限為2000 m,剔除遮擋因素角度向范圍為60°左右。原始數(shù)據(jù)采集由于雷達(dá)架設(shè)的高度不同,因此相同掠射角下所選擇的徑向距離點(diǎn)數(shù)不同。根據(jù)仿真要求選擇掠射角2~10°間的區(qū)域。所以岸基選1200~1700 m,艦載選600~1100 m作為數(shù)據(jù)處理區(qū)。
圖2是8個(gè)方向卷積之后疊加輸出結(jié)果。
圖2 經(jīng)過卷積后八幅疊加的梯度圖像(風(fēng)速為13.26 m/s,波高4.0 m)
將圖像前N%=10%的像素值作為梯度閾值檢測邊緣,用τF=6去除孤立點(diǎn)。通過灰度直方圖前H=30%判別獲得過渡區(qū)(陰影邊緣),最后建立邊緣的直方圖統(tǒng)計(jì)取眾數(shù)得到閾值τS,如圖3所示。
圖3 直方圖注:較高大的曲線為所選區(qū)域的直方圖曲線,較低的內(nèi)部曲線為邊界直方圖曲線,虛線是所得的閾值。
通過大量的數(shù)據(jù)試驗(yàn)重復(fù)上述方法求得陰影比例值,用Rune方法繪出波高與照度間的曲線。圖像表明在不同波高下掠射角越大陰影比例差異越小,掠射角越小在相同掠射角下陰影比例差異越大,特別是在有效掠射角范圍1~3°內(nèi)最為顯著,而小于1°的范圍內(nèi)背景噪聲較大故不作考慮。
圖4[16] 波高與照度曲線注:左邊在有效波高為2 m,風(fēng)速高于4 m/s,右邊有效波高4 m,風(fēng)速15 m/s,明顯看到在有效的掠射角范圍1~10°陰影比例隨波高變化最顯著是1~3°。
獲取閾值求取陰影比例后,將wavex測得的波高為外部參照。本文用了200組數(shù)據(jù)其中有10組海上數(shù)據(jù),海上數(shù)據(jù)選取波高值大都在3 m以上,風(fēng)速在15 m/s左右,剩余的都是岸基數(shù)據(jù),岸基海況上述已作介紹。利用其中一半數(shù)據(jù)來擬合直線,利用最小二乘法式(12)、式(13)估測出式(11)中A=0.693、B=2.83剩余的數(shù)據(jù)用來檢驗(yàn)。擬合效果如圖5所示。
圖5 利用最小二乘法擬合的直線注:橫坐標(biāo)表示陰影比例、縱坐標(biāo)表示有效波高。
利用最小二乘法得到線性關(guān)系式來估測有效波高,通過關(guān)系式反演的有效波高與外部參考波高作比較,評價(jià)算法可靠性,如圖6所示。
圖6 算法估測值Hm02與外部參考值比較
本文算法得到波高與外部參考波高求相關(guān)系數(shù)為0.89,均方根誤差為0.42,方差為0.21。通過比較統(tǒng)計(jì)參數(shù)得出本算法能準(zhǔn)確估測雷達(dá)圖像波高,而且參照值偏差較小,滿足估測精度要求。
目前對于有效波高的估測在工程上較為成熟的算法是譜分析法。圖像統(tǒng)計(jì)法大多停留在實(shí)驗(yàn)階段,本文算法得益于信噪比求波高啟發(fā),在圖像統(tǒng)計(jì)法中應(yīng)用線性關(guān)系結(jié)合外部參照來反演有效波高。
并利用實(shí)測雷達(dá)圖像數(shù)據(jù)來檢測算法精度,通過外部參考統(tǒng)計(jì)分析得到其相關(guān)系數(shù)為0.89,均方根誤差為0.42,方差為0.21。
試驗(yàn)結(jié)果表明:陰影比例與有效波高之間滿足線性關(guān),利用線性關(guān)系式能準(zhǔn)確估測有效波高,文獻(xiàn)[16]中提出假設(shè)得到驗(yàn)證,算法具有可行性與有效性。
然而由于數(shù)據(jù)有限算法仍有不足,算法沒有考慮電磁波散射模型對圖像陰影調(diào)制的影響。雷達(dá)極化方式、電磁波長、粗糙度等對陰影的作用。還需要進(jìn)一步考慮不同海域,低海況,低風(fēng)速等條件下算法適用性及可靠性。