孔德華, 張 東, 2, 張 卓, 宋志堯, 4

基于加速魯棒特征圖像匹配的云導(dǎo)風(fēng)計算方法

孔德華1, 張 東1, 2, 張 卓2, 3, 宋志堯2, 3, 4

(1. 南京師范大學(xué) 海洋科學(xué)與工程學(xué)院, 江蘇 南京 210023; 2. 江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心, 江蘇 南京 210023; 3. 南京師范大學(xué) 虛擬地理環(huán)境教育部重點實驗室, 江蘇 南京 210023; 4. 大規(guī)模復(fù)雜系統(tǒng)數(shù)值模擬江蘇省重點實驗室, 江蘇 南京 210023)

利用云導(dǎo)風(fēng)技術(shù)結(jié)合高分辨率氣象衛(wèi)星遙感數(shù)據(jù)獲取風(fēng)矢量, 在監(jiān)測臺風(fēng)等極端氣象災(zāi)害方面具有重要應(yīng)用。本文提出了一種基于加速魯棒特征(speeded up robust features, SURF)圖像匹配的云導(dǎo)風(fēng)計算方法, 利用SURF算法結(jié)合隨機(jī)抽樣一致算法(random sample consensus, RANSAC), 提取并匹配兩景連續(xù)時序云圖的特征點, 計算風(fēng)矢量, 并結(jié)合當(dāng)?shù)卮髿鉁囟壤€指定云高, 經(jīng)質(zhì)量控制得到云導(dǎo)風(fēng)矢量。運用該方法模擬了2018年臺風(fēng)“山竹”的云導(dǎo)風(fēng)矢量, 以美國威斯康星大學(xué)氣象衛(wèi)星研究合作所(CIMSS)的大氣運動矢量資料進(jìn)行驗證, 結(jié)果表明: (1) 風(fēng)速和風(fēng)向的相關(guān)系數(shù)分別為0.78和0.79, 均方根誤差分別為4.75 m·s–1和37.64°, 平均絕對百分比誤差分別為33.49%和22.55%, 整體具有良好的模擬精度; (2) 與CIMSS資料相比, 基于特征點匹配的SURF云導(dǎo)風(fēng)計算方法在反演密集云區(qū)的風(fēng)矢量有明顯優(yōu)勢, 可有效提高云區(qū)內(nèi)風(fēng)矢量的數(shù)量, 擴(kuò)大風(fēng)矢量的空間覆蓋范圍; (3) 圖像對比度增強(qiáng)處理對特征點的提取和風(fēng)矢量的空間分布有重要影響, 伽馬變換因子γ=5時, 能較好地平衡臺風(fēng)外圍螺旋云帶和中心附近云區(qū)的風(fēng)矢量數(shù)量, 反映臺風(fēng)風(fēng)場的整體特征。該方法作為基于尺度不變特征變換的云導(dǎo)風(fēng)計算方法的改進(jìn), 可為利用衛(wèi)星遙感影像數(shù)據(jù)進(jìn)行云導(dǎo)風(fēng)計算提供新的思路。

加速魯棒特征算法; 圖像匹配; 云導(dǎo)風(fēng); 風(fēng)矢量; 臺風(fēng)

云導(dǎo)風(fēng)技術(shù)是通過測算不同時相的氣象衛(wèi)星影像中云的運動來估計空間大范圍風(fēng)場的風(fēng)速和風(fēng)向的技術(shù), 所得結(jié)果稱為大氣運動矢量或風(fēng)矢量。風(fēng)矢量資料廣泛應(yīng)用于風(fēng)場數(shù)值模擬[1]、數(shù)值天氣預(yù)報[2]等研究領(lǐng)域, 在臺風(fēng)、暴雨等各種氣象災(zāi)害的監(jiān)測、分析和預(yù)測上發(fā)揮著重要的作用[3-4]。

云導(dǎo)風(fēng)技術(shù)最初是利用動畫技術(shù), 通過人工觀測衛(wèi)星云圖中云的移動, 獲得風(fēng)場數(shù)據(jù), 準(zhǔn)確性較低[5]。20世紀(jì)70年代, Leese等[6]提出了相關(guān)系數(shù)法, 利用連續(xù)衛(wèi)星云圖間圖像塊的相似性, 基于模板匹配技術(shù)追蹤示蹤云圖像塊的運動, 得出風(fēng)矢量。相關(guān)系數(shù)法作為云導(dǎo)風(fēng)技術(shù)的重要算法沿用至今, 但其運算量較大, 運算效率較低[7]。在探索不同于模板匹配的云導(dǎo)風(fēng)技術(shù)過程中, 王振會等[8]利用傅立葉相位分析法對示蹤云進(jìn)行頻域波譜分析, 通過諧波的相位變化計算波速, 得到風(fēng)矢量, 從而避免了相關(guān)系數(shù)法產(chǎn)生的“亞像元尺度位移”問題[9], 但是造成的相位重疊使得計算風(fēng)速偏小[10]。馬俠霖等[11]提出了基于尺度不變特征變換(scale invariant feature transform, SIFT)的云導(dǎo)風(fēng)計算方法, 通過在連續(xù)兩幅云圖中追蹤圖像特征點的位置變化得到風(fēng)矢量。該方法以像元而非圖像塊進(jìn)行匹配, 克服了示蹤云模板的制約, 因此不受云團(tuán)形變的影響, 得到的風(fēng)矢量精度較高。但在實際應(yīng)用中發(fā)現(xiàn), SIFT云導(dǎo)風(fēng)方法得到的風(fēng)矢量數(shù)量較少, 難以體現(xiàn)風(fēng)場的整體特征。

風(fēng)矢量數(shù)量由圖像上匹配到的特征點對數(shù)量的多少決定, SIFT法造成風(fēng)矢量較少的原因與圖像特征點提取算法有關(guān)。在眾多圖像特征提取算法中, 由Hebert Bay等[12]提出的加速魯棒特征算法(speeded up robust features, SURF)是對SIFT算法的改進(jìn), 在保證了尺度和旋轉(zhuǎn)的魯棒性的同時, 大大減少了算法的運算量, 被廣泛應(yīng)用于圖像匹配[13]、配準(zhǔn)[14]、跟蹤[15]等方面, 取得了良好效果。因此, 本文擬從圖像灰度特征出發(fā), 運用SURF算法提取相鄰時序臺風(fēng)影像的特征點, 結(jié)合特征點匹配方法與隨機(jī)抽樣一致算法(random sample consensus, RANSAC)實現(xiàn)風(fēng)矢量的追蹤, 并在此基礎(chǔ)上構(gòu)建完整的云導(dǎo)風(fēng)計算方法, 探討其在高時空分辨率氣象衛(wèi)星遙感影像中的風(fēng)場提取能力, 為利用氣象衛(wèi)星遙感技術(shù)快速獲取臺風(fēng)風(fēng)場信息, 監(jiān)測和預(yù)測臺風(fēng)運動, 減少極端氣象災(zāi)害帶來的危害提供方法和數(shù)據(jù)支撐。

1 基礎(chǔ)數(shù)據(jù)

2018年第22號臺風(fēng)“山竹”于2018年9月7日在西北太平洋上生成, 11日晚發(fā)展為超強(qiáng)臺風(fēng), 15日上午轉(zhuǎn)為強(qiáng)臺風(fēng), 16日17時登陸我國廣東省江門市臺山海宴鎮(zhèn), 登陸時中心附近最大風(fēng)力14級(45 m·s–1), 是2018年登陸我國的最強(qiáng)臺風(fēng)[16]。本文以臺風(fēng)“山竹”為例, 收集了以下三類數(shù)據(jù), 實現(xiàn)基于云導(dǎo)風(fēng)方法的風(fēng)矢量計算試驗: 1) 高分辨率氣象衛(wèi)星影像數(shù)據(jù); 2) 大氣溫度垂直廓線數(shù)據(jù); 3) 風(fēng)矢量精度驗證數(shù)據(jù)。

1.1 衛(wèi)星影像數(shù)據(jù)

葵花-8號(Himawari-8)氣象衛(wèi)星是日本氣象廳在2014年10月發(fā)射的新一代地球同步靜止氣象衛(wèi)星[17], 全盤掃描觀測時間間隔10 min, 影像空間分辨率0.5～2 km, 共有16個觀測波段, B13紅外波段(波長10.4 μm)、B3可見光波段(波長0.64 μm)、B9水汽波段(波長7.0 μm)被用于業(yè)務(wù)化的云導(dǎo)風(fēng)矢量計算[18]。紅外波段具有全天時、全天候成像能力, 同時也是追蹤云層運動的有效波段[19]。因此, 本文選用Himawari-8氣象衛(wèi)星B13紅外波段的亮溫影像轉(zhuǎn)換得到的灰度影像來進(jìn)行風(fēng)矢量的反演。針對2018年第22號臺風(fēng)“山竹”的登陸過程, 選取了北京時間2018年9月15日8時、14時和20時; 9月16日2時、8時和14時的6個整點時刻及各個整點的后10 min時刻共計12個時刻的衛(wèi)星影像, 按整點時刻分為6組, 每組包含前后間隔10 min的兩幅影像。所有影像經(jīng)過輻射定標(biāo)、幾何校正后, 裁剪出臺風(fēng)“山竹”所在區(qū)域, 用于模擬臺風(fēng)登陸前后的云導(dǎo)風(fēng)風(fēng)場。

1.2 大氣溫度垂直廓線數(shù)據(jù)

大氣溫度垂直廓線數(shù)據(jù)來自美國國家環(huán)境預(yù)報中心(NCEP)網(wǎng)站(https: //rda.ucar.edu)的全球再分析資料, 該資料包括地表溫度和大氣溫度廓線, 大氣垂直廓線分為26個高度層, 分別對應(yīng)1000、975、950、925、900、850、800、750、700、650、600、550、500、450、400、350、300、250、200、150、100、70、50、30、20和10 hpa氣壓高度[20], 空間分辨率為1°×1°。本文選取覆蓋了臺風(fēng)影像區(qū)域的六個整點時刻的全球再分析資料, 從中獲得大氣溫度垂直廓線數(shù)據(jù), 用于風(fēng)矢量的高度指定。

1.3 風(fēng)矢量驗證數(shù)據(jù)

用于對比和驗證的云導(dǎo)風(fēng)資料來自美國威斯康星大學(xué)氣象衛(wèi)星研究合作所(CIMSS)網(wǎng)站(http:// tropic.ssec.wisc.edu)的大氣運動矢量(Atmospheric Motion Vectors, AMVs)數(shù)據(jù)集。該資料包含有風(fēng)矢量的經(jīng)緯度、風(fēng)速、風(fēng)向、氣壓高度等信息。搜集了與實驗結(jié)果相對應(yīng)的6個整點時刻的紅外波段大氣運動矢量資料, 用于驗證和評價本實驗得到的風(fēng)矢量精度。

2 研究方法

本文提出基于SURF算法的云導(dǎo)風(fēng)計算方法, 將圖像增強(qiáng)和特征點匹配技術(shù)相結(jié)合, 運用在時序衛(wèi)星云圖上, 通過紅外波段云圖反演出風(fēng)矢量, 具體技術(shù)路線如圖1所示。技術(shù)流程圖中的關(guān)鍵技術(shù)方法概要介紹如下。實驗中使用的SURF算法、匹配算法以及RANSAC算法通過Python語言調(diào)用OpenCV開源庫中的相應(yīng)方法實現(xiàn)。

2.1 圖像增強(qiáng)

對比度強(qiáng)、細(xì)節(jié)特征明顯的圖像有利于特征點提取。從收集到的Himawari-8衛(wèi)星紅外波段影像可以看出(圖2a), 影像中云區(qū)像元的灰度值較高, 無云區(qū)像元的灰度值較低。特別是對于臺風(fēng)區(qū)域, 云區(qū)的覆蓋范圍大且密集, 云區(qū)像元的灰度值接近, 使得圖像整體偏亮且對比度不高。為提高圖像對比度, 突出云區(qū)細(xì)節(jié)特征, 采用伽馬變換對圖像進(jìn)行增強(qiáng)處理。伽馬變換是對輸入灰度值in進(jìn)行的非線性操作, 變換公式如下:

式中, 輸入灰度值Iin歸一化至[0, 1]區(qū)間; c為灰度縮放系數(shù), 本文取c=255, 使輸出灰度值Iout的取值范圍在0～255之間。γ為伽馬因子, 當(dāng)γ<1時, 圖像整體變亮; 反之圖像整體變暗。經(jīng)試驗, 取γ=5可取得好的圖像增強(qiáng)效果, 具體對比如圖2a、2b所示?？梢钥吹? 經(jīng)伽馬變換后, 臺風(fēng)外圍云區(qū)中的圖像細(xì)節(jié)被壓抑, 而密集云區(qū)中的圖像對比度增加, 細(xì)節(jié)得到增強(qiáng)。

a. 臺風(fēng)云圖; b. 圖像增強(qiáng)后的臺風(fēng)云圖

2.2 特征點提取

根據(jù)特征點在兩幅相鄰時刻臺風(fēng)云圖中的位置變化, 可以計算得到其所在像元上的風(fēng)矢量, 因此特征點提取算法對最終的云導(dǎo)風(fēng)計算結(jié)果至關(guān)重要。SURF算法[21]基于不同尺度下的近似Hessian矩陣檢測圖像特征點, 尺度空間分成4組, 每組包含4層圖像。圖中某點(,)的近似Hessian矩陣surf定義為:

式中,D,D,D分別是高斯濾波模板簡化后得到的盒子濾波器與積分圖像函數(shù)在該點處的卷積;為尺度因子,為以像元個數(shù)計數(shù)的盒子濾波器的邊長。由此可得近似Hessian矩陣的行列式如下:

det(surf) = DD– (0.9D)2, (4)

對尺度空間中的每層圖像計算近似Hessian行列式得到局部極值點, 以相鄰三層圖像的中間層的每個局部極值點為中心, 在當(dāng)前層和上、下層中分別選取該點周圍3×3鄰域內(nèi)的像元, 構(gòu)成3×3×3的立體鄰域。若該極值大于立體鄰域其他26個像元的近似Hessian行列式值, 則確定該局部極值點為特征點。

得到圖像特征點后需要確定特征點的主方向和特征向量, 作為特征點匹配時的依據(jù)。以特征點為圓心, 對半徑為6的圓形區(qū)域內(nèi)的像素進(jìn)行Haar小波變換, Haar小波濾波器模板如圖3所示, 模板邊長為4, 其中黑色區(qū)域權(quán)重為–1, 白色區(qū)域權(quán)重為1。統(tǒng)計所有特征點在和方向上的Haar小波響應(yīng)值, 以60°的扇形窗口遍歷整個圓形區(qū)域, 計算窗口內(nèi)的特征點的Harr小波響應(yīng)值相加形成的局部方向向量, 最長向量方向則為該特征點主方向。如圖4所示, 紅色箭頭方向即為該特征點的主方向。

注: a. Haar小波在方向上的濾波器模板; b. Haar小波在方向上的濾波器模板

特征向量的計算過程如圖5所示, 以每個特征點為中心, 構(gòu)建邊長為20的正方形窗口, 并旋轉(zhuǎn)到特征點的主方向(圖5中紅色箭頭方向), 將其分為16個邊長為5的正方形子區(qū)域, 計算每個子區(qū)域的像元在水平方向和垂直方向上的Haar小波響應(yīng)值之和Σd、Σd及其絕對值之和Σ|d|與Σ|d|, 得到一個4維特征向量= {Σd, Σd, Σ|d|, Σ|d|}, 串接16子區(qū)域的特征向量, 得到該特征點的64維特征向量。

2.3 特征點匹配

分別提取出兩幅時序云圖中的所有特征點后, 基于特征向量的相似度, 對前一時刻影像的每個特征點, 在后一時刻找到與之最相似的特征點進(jìn)行匹配。選擇歐氏距離Dis作為特征向量相似度的評判依據(jù), 歐氏距離最小時認(rèn)定兩點的特征向量最相似, 兩個特征點相互匹配。Dis計算公式如下:

式中,V表示前一時刻影像中第個特征點對應(yīng)的特征向量的第個元素;V表示后一時刻影像中第個特征點對應(yīng)的特征向量的第個元素。

2.4 誤匹配過濾

利用歐氏距離進(jìn)行特征點匹配的過程中, 會存在少量特征向量最相似的特征點對并非合理匹配的情況。為提高特征點匹配的準(zhǔn)確度, 減少后續(xù)風(fēng)矢量計算時出現(xiàn)異常值的概率, 本文采用RANSAC算法來減少特征點誤匹配情況[22]。RANSAC算法的具體步驟如下:

1) 從全部組匹配點對中隨機(jī)提取4組匹配點對, 求解對應(yīng)的單應(yīng)性變換矩陣, 其表達(dá)式為:

2) 根據(jù)該變換矩陣, 對前一時刻特征點集中的其余-4個特征點進(jìn)行變換, 計算變換后的每個點的坐標(biāo)與原匹配點的坐標(biāo)之間的距離。若距離小于閾值, 則認(rèn)為該特征點為此變換下的內(nèi)點, 否則為外點, 記錄內(nèi)點數(shù)量, 閾值設(shè)定為5個像元的長度。

3) 重復(fù)1-2步驟若干次, 內(nèi)點數(shù)量最多時的變換矩陣′即為正確的變換矩陣。

4) 將正確的變換矩陣下的內(nèi)點保留, 外點去除, 實現(xiàn)對特征匹配結(jié)果的優(yōu)化, 完成特征點誤匹配過濾。

2.5 風(fēng)矢量高度指定

得到正確匹配的特征點對結(jié)果后, 根據(jù)已知時間間隔內(nèi)每組匹配點對中的兩個特征點間的位置變化, 計算出風(fēng)矢量的大小和方向, 得到初始風(fēng)矢量。由于風(fēng)矢量在云頂高度附近, 因此推算云高即可得出風(fēng)矢量的高度估計值。云頂高度通過結(jié)合大氣溫度垂直廓線數(shù)據(jù)和紅外云圖云頂亮溫值來推算[23], 其中大氣溫度垂直廓線中的溫度采用經(jīng)轉(zhuǎn)換后的亮度溫度。設(shè)風(fēng)矢量所在位置的云頂亮溫值為(K), 對應(yīng)的氣壓高度為(hpa), 風(fēng)矢量高度計算的步驟如下:

1) 根據(jù)風(fēng)矢量的坐標(biāo), 在大氣溫度垂直廓線數(shù)據(jù)中提取該位置的26個高度層對應(yīng)的大氣溫度集合emp = {1,2,3, …,26}, 其中1>2>3>…>26。

2) 對于風(fēng)矢量所在像元的云頂亮溫值, 若能在集合emp中找到與相等的亮溫值, 則直接查得所對應(yīng)的氣壓高度值;

3) 若不能在集合emp中直接找到與相等的亮溫值, 則通過對數(shù)線性插值法計算對應(yīng)的氣壓高度。對數(shù)線性插值公式如下:

式中,1、2分別為亮溫值1、2對應(yīng)的大氣溫度垂直廓線高度層的氣壓高度。1、2的取值分三種情況: 若>1, 則令1=2,2=1; 若<26, 則令1=26,2=25; 若26<<1, 則將集合中與最相近的兩個亮溫值分別記為1、2, 且2<<1。

2.6 風(fēng)矢量質(zhì)量控制

為確保每個風(fēng)矢量與周圍風(fēng)矢量的屬性差異在合理范圍內(nèi), 分別通過計算每個風(fēng)矢量與其鄰域內(nèi)風(fēng)矢量的風(fēng)速均方根誤差和風(fēng)向均方根誤差, 設(shè)置閾值剔除誤差較大的風(fēng)矢量, 實現(xiàn)對初始風(fēng)矢量的質(zhì)量控制[10]。Holmlund認(rèn)為在進(jìn)行風(fēng)矢量質(zhì)量控制時, 以目標(biāo)矢量為中心的正方形鄰域的邊長應(yīng)設(shè)定為100～200 km[24]。由于實驗采用影像的空間分辨率為4 km, 故設(shè)定目標(biāo)風(fēng)矢量正方形鄰域搜索范圍的邊長為50個像元。Endlic和McLean發(fā)現(xiàn)對于高空急流, 在1°的緯度范圍內(nèi), 風(fēng)速的衰減約為15%, 且氣流越強(qiáng)風(fēng)速的衰減越大[25], 因此風(fēng)速均方根誤差的閾值設(shè)定為正方形鄰域內(nèi)最大風(fēng)速的40%。風(fēng)向均方根誤差的閾值參考馬俠霖等的取值, 設(shè)定為100°[10]。風(fēng)速均方根誤差R和風(fēng)向均方根誤差R的計算公式如下:

式中,、為目標(biāo)風(fēng)矢量的風(fēng)速和風(fēng)向;v、d為鄰域風(fēng)矢量的風(fēng)速和風(fēng)向;為鄰域內(nèi)除目標(biāo)風(fēng)矢量以外的風(fēng)矢量個數(shù)。

3 結(jié)果

3.1 臺風(fēng)“山竹”的云導(dǎo)風(fēng)矢量結(jié)果

利用收集到的臺風(fēng)“山竹”12個時刻的Himawari-8衛(wèi)星遙感影像數(shù)據(jù), 采用SURF云導(dǎo)風(fēng)方法計算得到6個整點時刻的風(fēng)矢量結(jié)果如圖6所示。對于這6個時刻臺風(fēng)“山竹”的云導(dǎo)風(fēng)矢量結(jié)果, 從風(fēng)速分布來看, 北部云區(qū)的風(fēng)速普遍大于南部云區(qū), 風(fēng)速超過30 m·s–1的風(fēng)矢量集中分布在臺風(fēng)中心附近以北的云區(qū), 且最大風(fēng)速均小于40 m·s–1。由這些區(qū)域向外, 風(fēng)速逐漸減小, 風(fēng)速小于5 m·s–1的風(fēng)矢量集中分布在南部的外圍螺旋云系。從風(fēng)速變化來看, 圖6(a)～(c)中, “山竹”西南側(cè)云系的風(fēng)速增大, 增強(qiáng)的氣流為中心云系輸送了能量, 風(fēng)眼逐漸出現(xiàn), 臺風(fēng)強(qiáng)度有增強(qiáng)趨勢, 眼區(qū)附近風(fēng)速逐漸增大; 圖6(d)中臺風(fēng)眼最為清晰, 圍繞眼區(qū)分布的風(fēng)矢量的最為密集; 圖6(e)～(f)中風(fēng)速較大的風(fēng)矢量分散到外圍云區(qū), 風(fēng)眼逐漸消散, 來自東南側(cè)的急流云系維持著向中心云區(qū)的能量輸送。計算得到的風(fēng)矢量主要分布在400 hpa氣壓高度以上, 150 hpa氣壓高度以下, 圖7所示為3個整點時刻中氣壓高度在150～200 hpa之間的風(fēng)矢量?？梢钥闯? 紅外波段臺風(fēng)影像中, 亮白云團(tuán)的高度較高, 臺風(fēng)中心附近的亮白云團(tuán)呈順時針向外輻散的運動趨勢?？傮w上, SURF云導(dǎo)風(fēng)方法計算得到的風(fēng)矢量能有效顯示臺風(fēng)發(fā)展變化過程中云團(tuán)運動速度的變化情況。

3.2 風(fēng)矢量精度驗證

由于探空測站主要位于陸地, 難以獲取臺風(fēng)在海上階段的實測風(fēng)速數(shù)據(jù), 因此采用CIMSS的大氣運動矢量數(shù)據(jù)集來對SURF云導(dǎo)風(fēng)方法反演的風(fēng)矢量進(jìn)行精度驗證。由于得到的風(fēng)矢量具有不同的高度, 因此精度驗證時, 以本方法所得風(fēng)矢量為基準(zhǔn), 在數(shù)據(jù)集中選擇垂直方向氣壓高度相差100 hpa以內(nèi)、水平方向經(jīng)緯度相距0.1°范圍內(nèi)的最鄰近的矢量, 采用平均絕對誤差(Mean Absolute Error, MAE)、平均絕對百分比誤差(Mean Absolute Percentage Error, MAPE)、均方根誤差(Root Mean Square Error, RMSE)和相關(guān)系數(shù)()四個指標(biāo)分別從風(fēng)速和風(fēng)向兩方面對風(fēng)矢量的精度進(jìn)行驗證, 結(jié)果如表1所示。六個時刻的風(fēng)速和風(fēng)向平均相關(guān)系數(shù)分別為0.78和0.79, 均方根誤差分別為4.75 m·s–1和37.64°, 平均絕對百分比誤差分別為33.49%和22.55%。風(fēng)速與風(fēng)向的絕對值誤差分布如圖8所示。風(fēng)速絕對值誤差小于6 m·s–1的風(fēng)矢量占總數(shù)量的82.6%, 風(fēng)向絕對值誤差小于40°的風(fēng)矢量占總數(shù)量的83.7%, 可見本方法所得風(fēng)矢量的風(fēng)速、風(fēng)向與CIMSS大氣運動矢量資料吻合良好。

注: a. 9月15日08時; b. 9月15日14時; c. 9月15日20時; d. 9月16日02時; e. 9月16日08時; f. 9月16日14時

注: a. 9月15日14時; b. 9月15日20時; c. 9月16日02時

表1 SURF云導(dǎo)風(fēng)方法風(fēng)矢量精度驗證

注: a. 風(fēng)速差絕對值; b. 風(fēng)向差絕對值

3.3 風(fēng)矢量空間分布對比

圖9顯示了SURF云導(dǎo)風(fēng)方法風(fēng)矢量與CIMSS風(fēng)矢量在所有高度上的空間分布對比?？梢钥闯? SURF云導(dǎo)風(fēng)方法所得風(fēng)矢量與CIMSS資料風(fēng)矢量在空間上具有不同的分布特征。從Himawari-8衛(wèi)星紅外波段臺風(fēng)影像上看, 在臺風(fēng)外圍區(qū)域, 云層較薄, 像元灰度值小, 圖像細(xì)節(jié)特征較少, 不利于特征點的提取和匹配, 因此基于圖像特征點追蹤的SURF云導(dǎo)風(fēng)方法得到的風(fēng)矢量明顯少于CIMSS資料風(fēng)矢量。而在臺風(fēng)云區(qū), 由于臺風(fēng)氣旋的螺旋型形態(tài), 導(dǎo)致不同云高的像元灰度值之間存在差異, 圖像對比度增大, 細(xì)節(jié)特征增強(qiáng), 因此在有CIMSS風(fēng)矢量的區(qū)域, SURF云導(dǎo)風(fēng)方法基本上均能反演得到風(fēng)矢量, 而在臺風(fēng)眼下方區(qū)域, SURF云導(dǎo)風(fēng)方法反演的風(fēng)矢量明顯多于CIMSS數(shù)據(jù)集在該區(qū)域的風(fēng)矢量。雖然這部分風(fēng)矢量缺少實測數(shù)據(jù)可以直接驗證其精度, 但是定性來看, 風(fēng)矢量的風(fēng)速大小與周圍的風(fēng)矢量相近, 風(fēng)矢量的方向也呈現(xiàn)明顯的逆時針螺旋型分布趨勢, 因此可以認(rèn)為其能夠代表當(dāng)前的臺風(fēng)風(fēng)場?？梢娤啾扔贑IMSS風(fēng)矢量資料, SURF云導(dǎo)風(fēng)方法在反演密集云區(qū)的風(fēng)矢量時具有較大優(yōu)勢。

4 討論

4.1 伽馬因子對特征點提取結(jié)果的影響

由于SURF算法基于圖像尺度空間的極值計算特征點, 圖像增強(qiáng)處理的結(jié)果直接影響圖像特征點的提取, 進(jìn)而影響風(fēng)矢量的計算結(jié)果以及空間分布。為探究圖像增強(qiáng)處理時伽馬因子的取值對后續(xù)圖像特征點計算的影響, 以確定合適的值, 分別對各個時刻影像取= 1, 2, 3, …, 12進(jìn)行伽馬變換和基于SURF的特征點提取。由于臺風(fēng)“山竹”云系龐大, 直徑范圍達(dá)1 000 km, 因此以臺風(fēng)中心為圓心、500 km半徑確定圓形區(qū)域, 將臺風(fēng)“山竹”的云圖分為內(nèi)外兩個云區(qū), 內(nèi)云區(qū)位于圓形區(qū)域內(nèi), 包含臺風(fēng)眼區(qū)、云墻區(qū)和部分螺旋云帶; 外云區(qū)位于圓形區(qū)域外, 主要為外圍螺旋云帶。圖10顯示了不同γ取值下內(nèi)、外云區(qū)中的特征點數(shù)量對比。

可以看出, 伽馬因子的取值大于1時, 隨的增大, 外云區(qū)的特征點數(shù)量不斷減少, 而內(nèi)云區(qū)的特征點數(shù)量則波動增加, 達(dá)到最大值后出現(xiàn)逐漸減少的趨勢。在多數(shù)時刻下, 大致在=5時內(nèi)云區(qū)的特征點數(shù)量會出現(xiàn)一個峰值。若繼續(xù)增大, 內(nèi)云區(qū)特征點的數(shù)量不再增加或增幅較小, 而外云區(qū)特征點的數(shù)量則持續(xù)減少。因此, 若從選取單一值的角度, 取=5, 可以有效突出內(nèi)云區(qū)臺風(fēng)眼和云墻區(qū)的風(fēng)矢量, 同時較好地維持臺風(fēng)風(fēng)場總的風(fēng)矢量數(shù)量, 體現(xiàn)整體臺風(fēng)云區(qū)的螺旋形風(fēng)場特征。

4.2 特征點提取方法對云導(dǎo)風(fēng)結(jié)果的影響

為研究特征點提取方法對云導(dǎo)風(fēng)結(jié)果的影響, 以SIFT云導(dǎo)風(fēng)方法對臺風(fēng)“山竹”進(jìn)行了風(fēng)矢量反演, 除了特征點提取方法不同以外, 其余數(shù)據(jù)處理步驟都與SURF云導(dǎo)風(fēng)方法一致。同樣以CIMSS大氣運動矢量資料對風(fēng)矢量反演結(jié)果進(jìn)行驗證, 結(jié)果如表2所示。SIFT云導(dǎo)風(fēng)方法得到的風(fēng)矢量的風(fēng)速和風(fēng)向平均相關(guān)系數(shù)分別為0.81和0.85, 均方根誤差分別為4.11 m·s–1和31.22°, 平均絕對百分比誤差分別為25.25%和21.75%。與SURF云導(dǎo)風(fēng)方法相比, 兩種方法的誤差比較接近, 總體上SIFT云導(dǎo)風(fēng)方法得到的風(fēng)矢量的精度上略高。

表2 SIFT云導(dǎo)風(fēng)方法風(fēng)矢量精度驗證

圖11顯示了基于SURF算法和SIFT算法得到的風(fēng)矢量反演結(jié)果的空間分布對比?？梢钥吹? SURF方法得到的風(fēng)矢量數(shù)量明顯多于SIFT方法的結(jié)果。SIFT云導(dǎo)風(fēng)方法得到的風(fēng)矢量集中分布在云區(qū)邊緣和臺風(fēng)中心附近區(qū)域, 而SURF云導(dǎo)風(fēng)方法得到的風(fēng)矢量在臺風(fēng)云區(qū)有更廣泛的分布, 特別是在臺風(fēng)眼附近及下方區(qū)域, 因此臺風(fēng)風(fēng)場的螺旋形形態(tài)特征也更明顯。造成這種差異的原因是SURF算法在圖像尺度空間構(gòu)造、特征點檢測的過程、特征向量方向的確定和特征描述符的生成四個方面優(yōu)化和改進(jìn)了SIFT算法, 在圖像的更多位置提取到了特征點, 最終得到了更多風(fēng)矢量。不同時刻下兩種方法得到的臺風(fēng)“山竹”的風(fēng)矢量數(shù)量對比如表3所示, 基于SURF的云導(dǎo)風(fēng)方法得到的風(fēng)矢量在數(shù)量上平均增加了49.6 %。因此, 在臺風(fēng)風(fēng)矢量模擬精度相近的情況下, SURF云導(dǎo)風(fēng)方法比SIFT云導(dǎo)風(fēng)方法對臺風(fēng)風(fēng)場有更好的代表性, 也能更好地描述臺風(fēng)風(fēng)場的空間形態(tài)。

5 結(jié)論

本文提出了一種基于快速魯棒特征(SURF)圖像匹配的云導(dǎo)風(fēng)計算方法, 在運用伽馬變換增強(qiáng)云區(qū)圖像對比度的基礎(chǔ)上, 利用SURF算法計算兩景相鄰時相臺風(fēng)紅外波段影像的特征點, 以歐氏距離最小的原則完成相似特征點的匹配; 進(jìn)一步引入RANSAC算法去除誤匹配點對, 根據(jù)每組匹配中兩點之間的位移計算出風(fēng)速的大小和方向, 結(jié)合當(dāng)?shù)卮髿鉁囟壤€指定云高, 經(jīng)質(zhì)量控制后得到最終的風(fēng)矢量。以臺風(fēng)“山竹”為試驗, 得出的研究結(jié)果表明:

(1) 以CIMSS數(shù)據(jù)集風(fēng)矢量為驗證, SURF云導(dǎo)風(fēng)計算方法得到的風(fēng)矢量的風(fēng)速和風(fēng)向的相關(guān)系數(shù)分別為0.78和0.79, 均方根誤差分別為4.75 m·s–1和37.64°, 平均絕對百分比誤差分別為33.49%和22.55%, 具有良好的模擬精度, 可有效用于臺風(fēng)風(fēng)場的模擬, 為極端氣象災(zāi)害的風(fēng)場數(shù)值模擬提供參考數(shù)據(jù)。

(2) 對比CIMSS數(shù)據(jù)集的風(fēng)矢量空間分布, SURF云導(dǎo)風(fēng)計算方法在反演密集云區(qū)的風(fēng)矢量時, 可以利用密集云區(qū)的特征點匹配, 得到更多的風(fēng)矢量。因此該方法可以彌補基于模板匹配的云導(dǎo)風(fēng)方法由于云團(tuán)形變原因等無法找到合適的示蹤云塊, 從而無法得到密集云區(qū)的風(fēng)矢量的缺陷, 進(jìn)一步擴(kuò)大云區(qū)內(nèi)風(fēng)矢量的空間覆蓋范圍。

(3) 對比基于尺度不變特征變換(SIFT)的云導(dǎo)風(fēng)計算方法, 基于快速魯棒特征(SURF)的云導(dǎo)風(fēng)計算方法得到的風(fēng)矢量的風(fēng)速和風(fēng)向模擬精度相近, 但是風(fēng)矢量的數(shù)量平均增加49.6%, 顯著擴(kuò)大了風(fēng)矢量的空間分布, 能夠更好地描述臺風(fēng)風(fēng)場的空間形態(tài)。

(4) 圖像增強(qiáng)處理對特征點的提取數(shù)量具有重要影響。當(dāng)伽馬因子≥1時,取值較小對于反演臺風(fēng)外圍螺旋云帶的風(fēng)矢量較為有利;取值較大適用于反演臺風(fēng)中心附近云區(qū)的風(fēng)矢量。從選取單一值的角度, 發(fā)現(xiàn)當(dāng)=5時, 能夠適中平衡臺風(fēng)外圍螺旋云帶和臺風(fēng)中心附近云區(qū)的風(fēng)矢量數(shù)量, 突出SURF云導(dǎo)風(fēng)計算方法在反演內(nèi)云區(qū)風(fēng)矢量時的優(yōu)勢, 同時兼顧外云區(qū)風(fēng)矢量, 較好地反映臺風(fēng)風(fēng)場的整體特征。此外, 若針對內(nèi)云區(qū)和外云區(qū)采用不同的, 風(fēng)矢量的數(shù)量可能進(jìn)一步增加, 從這一角度出發(fā),的選取值得進(jìn)一步研究。

致謝: 感謝河海大學(xué)海洋學(xué)院葵花衛(wèi)星地面接收站提供的4 km分辨率的臺風(fēng)“山竹”Himawari-8氣象衛(wèi)星紅外波段影像數(shù)據(jù); 感謝美國國家環(huán)境預(yù)報中心(NCEP)提供的全球再分析資料和美國威斯康星大學(xué)氣象衛(wèi)星研究合作所(CIMSS)的大氣運動矢量數(shù)據(jù)。

[1] 何斌, 翟國慶, 高坤, 等. 云跡風(fēng)資料同化對東亞海域風(fēng)場數(shù)值模擬的影響[J]. 海洋學(xué)報(中文版), 2007, 29(6): 23-32.

HE Bin, ZHAI Guoqing, GAO Kun, et al. The impact of cloud motion wind data on numerical simulation of wind field in the sea area of Eastern Asia[J]. Acta Ocea-no-lo-gica Sinica, 2007, 29(6): 23-32.

[2] 李華宏, 王曼, 薛紀(jì)善, 等. FY-2C云跡風(fēng)資料在中尺度數(shù)值模式中的應(yīng)用研究[J]. 氣象學(xué)報, 2008, 66(1): 50-58.

LI Huahong, WANG Man, XUE Jishan, et al. A study on the application of FY-2C cloud drift wind in a meso-scale numerical model[J]. Acta Meteorologica Sinica, 2008, 66(1): 50-58.

[3] 劉瑞, 翟國慶, 王彰貴, 等. FY-2C云跡風(fēng)資料同化應(yīng)用對臺風(fēng)預(yù)報的影響試驗研究[J]. 大氣科學(xué), 2012, 36(2): 350-360.

LIU Rui, ZHAI Guoqing, WANG Zhanggui, et al. Impact of application of cloud motion wind data from FY-2C satellite on simulation of typhoon cases[J]. Chi-nese Journal of Atmospheric Sciences, 2012, 36(2): 350-360.

[4] 周兵, 徐海明, 吳國雄, 等. 云跡風(fēng)資料同化對暴雨預(yù)報影響的數(shù)值模擬[J]. 氣象學(xué)報, 2002, 60(3): 309-317.

ZHOU Bing, XU Haiming, WU Guoxiong, et al. Numerical simulation of CMWDA with it’s impacting on torrential rain forecast[J]. Acta Meteorologica Sinica, 2002, 60(3): 309-317.

[5] IZAWA T, FUJITA T. Relationship between observed winds and cloud velocities determined from pictures obtained by the ESSA III, ESSA V and ATS-I satellites[J]. Space Research IX, Amsterdam, North-Holland Publ. Co, 1969: 571-579.

[6] LEESE J A, NOVAK C S, CLARK B B. An automated technique for obtaining cloud motion from geosynchronous satellite data using cross correlation[J]. Journal of Applied Meteorology and Climatology, 1971, 10(1): 118-132.

[7] 許健民, 張其松. 衛(wèi)星風(fēng)推導(dǎo)和應(yīng)用綜述[J]. 應(yīng)用氣象學(xué)報, 2006, 17(5): 574-582.

XU Jianmin, ZHANG Qisong. Status review on atmospheric motion vectors-derivation and application[J]. Journal of Applied Meteorological Science, 2006, 17(5): 574-582.

[8] 王振會, 許建明, G.KELLY. 基于傅立葉相位分析的衛(wèi)星云圖導(dǎo)風(fēng)技術(shù)[J]. 氣象科學(xué), 2004, 24(1): 9-15.

WANG Zhenhui, XU Jianming, G.KELLY. Deriving cloud motion vectors from high temporal resolution ima-ges based on Fourier phase analysis technique[J]. Journal of the Meteorological Sciences, 2006, 17(5): 574-582.

[9] PURDOM J F. W. Detailed cloud motions from satellite imagery taken at thirty second, one and three minute intervals[C]//Proceeding to the 3rd international wind workshop in Ascona, Switzerland. 1996: 10-12.

[10] 許建明, 王振會. 傅立葉相位分析導(dǎo)風(fēng)技術(shù)的適用范圍和相位重疊分析[J]. 氣象科學(xué), 2004, 24(3): 309-313.

XU Jianming, WANG Zhenhui. The application condition of Fourier analysis technique and the analysis of phrase overlap[J]. Journal of the Meteorological Sciences, 2004, 24(3): 309-313.

[11] 馬俠霖, 羅鵬, 陳志斌, 等. 基于尺度不變特征變換的云導(dǎo)風(fēng)計算方法[J]. 氣象科技, 2014, 42(3): 391-396.

MA Xialin, LUO Peng, CHEN Zhibin, et al. A method for cloud motion wind vector prediction based on Scale- Invariant Feature Transform[J]. Meteorological Science and Technology, 2014, 42(3): 391-396.

[12] BAY H, TUYTELAARS T, VAN GOOL L. Surf: Spee-ded up robust features[C]//European conference on com-puter vision. Springer, Berlin, Heidelberg, 2006: 404-417.

[13] 黃春鳳, 劉守山, 別治峰, 等. 改進(jìn)的SURF算法在圖像匹配中的應(yīng)用[J]. 現(xiàn)代電子技術(shù), 2020, 43(10): 111-115.

HUANG Chunfeng, LIU Shoushan, BIE Zhifeng, et al. Application of improved SURF algorithm in image ma-tching[J]. Modern Electronics Technique, 2020, 43(10): 111-115.

[14] 唐穎復(fù), 王忠靜, 張子雄. 基于改進(jìn)SIFT和SURF算法的沙丘圖像配準(zhǔn)[J]. 清華大學(xué)學(xué)報(自然科學(xué)版), 2021, 61(2): 161-169.

TANG Yingfu, WANG Zhongjing, ZHANG Zixiong. Re-gistration of sand dune images using an improved SIFT and SURF algorithm[J]. Journal of Tsinghua University (Science and Technology), 2021, 61(2): 161-169.

[15] 權(quán)巍, 包鐵壯, 白寶興, 等. SURF與RANSAC算法結(jié)合的圖像跟蹤方法[J]. 計算機(jī)仿真, 2016, 33(9): 268-272.

QUAN Wei, BAO Tiezhuang, BAI Baoxing, et al. A method of real-time picture tracking and mapping ba-sed on SURF algorithm and RANSAC algorithm[J]. Computer Simulation, 2016, 33(9): 268-272.

[16] 方艷瑩, 錢燕珍, 申華羽, 等. 1822號臺風(fēng)“山竹”引起浙江東北部大暴雨成因分析[J]. 海洋預(yù)報, 2020, 37(4): 86-96.

FANG Yanying, QIAN Yanzhen, SHEN Huayu, et al. Causes analysis on the heavy rainfall in northeastern Zhejiang related to the typhoon “Mangkhut” (1822)[J]. Marine Forecasts, 2020, 37(4): 86-96.

[17] 周旋, 葉小敏, 周江濤, 等. 基于Himawari-8衛(wèi)星的逐時次海表溫度融合[J]. 海洋學(xué)報, 2021, 43(1): 137-146.

ZHOU Xuan, YE Xiaomin, ZHOU Jiangtao, et al. Hourly sea surface temperature fusion based on Himawari-8 sa-tel-lite[J]. Acta Oceanologica Sinica, 2021, 43(1): 137- 146.

[18] SHIMOJI K. Introduction to the Himawari-8 atmospheric motion vector algorithm[J]. Meteorological Satellite Center Technical Note, 2017, 62: 73-77.

[19] OYAMA R, SAWADA M, SHIMOJI K. Diagnosis of tropical cyclone intensity and structure using upper tropospheric Atmospheric Motion Vectors[J]. Journal of the Meteorological Society of Japan. Ser. II, 2018, 96: 3-26.

[20] 王靜, 趙朝方. 最小方差算法反演NOAA/AMSU海面大氣溫度垂直廓線[J]. 海洋湖沼通報, 2008(2): 16-23.

WANG Jing, ZHAO Chaofang. Retrieval of atmosphe-ric vertical temperature profile over sea from NOAA/AMSU data with minimum variance algorithm[J]. Transactions of Oceanology and Limnology, 2008(2): 16-23.

[21] BAY H, ESS A, TUYTELAARS T, et al. Speeded-up robust features (SURF)[J]. Computer vision and image understanding, 2008, 110(3): 346-359.

[22] FISCHLER M A, BOLLES R C. Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography [J]. Com-munications of the ACM, 1981, 24(6): 381-395.

[23] 許建明. 一維傅立葉相位導(dǎo)風(fēng)技術(shù)的研究[D]. 南京: 南京氣象學(xué)院, 2003.

XU Jianming. Research on one-dimensional Fourier phase wind guide technology[D]. Nanjing: Nanjing me-teorological College, 2003.

[24] HOLMLUND K. The utilization of statistical properties of satellite-derived atmospheric motion vectors to derive quality indicators [J]. Weather and Forecasting, 1998, 13(4): 1093-1104.

[25] ENDLICH R M, MCLEAN G S. The structure of the jet stream core [J]. Journal of Atmospheric Sciences, 1957, 14(6): 543-552.

A cloud-motion-based wind retrieval method based on the SURF algorithm

KONG De-hua1, ZHANG Dong1, 2, ZHANG Zhuo2, 3, SONG Zhi-yao2, 3, 4

(1. College of Marine Science and Engineering, Nanjing Normal University, Nanjing 210023, China; 2. Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application, Nanjing 210023, China; 3. Key Lab of Virtual Geographic Environment under the Ministry of Education, Nanjing Normal University, Nanjing 210023, China; 4. Jiangsu Provincial Key Laboratory for Numerical Simulation of Large Scale Complex System, Nanjing 210023, China)

Cloud-motion-based wind retrieval technology combined with high-resolution meteorological satellite remote sensing data to obtain wind vectors has crucial applications in monitoring extreme meteorological disasters such as typhoons. In this study, a cloud-motion-based wind retrieval method based on the speeded-up robust features (SURF) image matching algorithm is proposed. The SURF and random sample consensus algorithms are used to extract and match the feature points of two consecutive time-series cloud images, calculate the wind vectors, and specify cloud height in combination with the local atmospheric temperature profile. The cloud-motion-based wind retrieval vectors are obtained through quality control. The proposed method is used to simulate the cloud- motion-based wind retrieval vectors of Typhoon “Mangkhut” in 2018, which is verified using the Cooperative Institute for Meteorological Satellite Studies (CIMSS) atmospheric motion vector data. The results indicate that: (1) The correlation coefficients of wind speed and wind direction are 0.78 and 0.79, respectively, the root mean square errors are 4.75 m·s?1and 37.64°, respectively, and the average absolute percentage errors are 33.49% and 22.55%, respectively, which has good simulation accuracy overall. (2) Compared with the CIMSS data, the cloud-motion- based wind retrieval method based on the SURF feature matching algorithm has obvious advantages in retrieving wind vectors in dense cloud areas, which can effectively improve the number of wind vectors in cloud areas and expand the spatial coverage of wind vectors. (3) Image contrast enhancement highly impacts the extraction of feature points and the spatial distribution of wind vector= 5, which can better balance the number of wind vectors in the spiral cloud belt and the area near the center of the typhoon, reflecting the overall characteristics of the typhoon wind field. As an improvement of the cloud-motion-based wind retrieval method based on scale- invariant feature transformation, the proposed method can be a new approach for wind vector calculation using satellite remote sensing image data.

SURF algorithm; image matching; cloud-motion-based wind retrieval; wind vectors; typhoon

Nov. 13, 2021

TP751

A

1000-3096(2022)09-0064-13

10.11759/hykx20211113001

2021-11-13;

2022-03-12

國家重點研發(fā)計劃項目(2018YFB0505500, 2018YFB0505502)

[National Key R & D Program of China, Nos. 2018YFB0505500, 2018YFB0505502]

孔德華(1996—), 男, 江蘇南京人, 碩士研究生, 主要從事風(fēng)矢量遙感反演研究, E-mail: 1290304214@qq.com; 張東(1975—), 男,通信作者, 江蘇南通人, 博士, 副教授, 研究方向為海洋信息技術(shù)與海岸帶資源開發(fā)管理, E-mail: zhangdong@njnu.edu.cn

(本文編輯: 康亦兼)