基于加速魯棒特征圖像匹配的云導風計算方法

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

基于加速魯棒特征圖像匹配的云導風計算方法

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

(1. 南京師范大學 海洋科學與工程學院, 江蘇 南京 210023; 2. 江蘇省地理信息資源開發與利用協同創新中心, 江蘇 南京 210023; 3. 南京師范大學 虛擬地理環境教育部重點實驗室, 江蘇 南京 210023; 4. 大規模復雜系統數值模擬江蘇省重點實驗室, 江蘇 南京 210023)

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

加速魯棒特征算法; 圖像匹配; 云導風; 風矢量; 臺風

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

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

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

1 基礎數據

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

1.1 衛星影像數據

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

1.2 大氣溫度垂直廓線數據

大氣溫度垂直廓線數據來自美國國家環境預報中心(NCEP)網站(https: //rda.ucar.edu)的全球再分析資料, 該資料包括地表溫度和大氣溫度廓線, 大氣垂直廓線分為26個高度層, 分別對應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°。本文選取覆蓋了臺風影像區域的六個整點時刻的全球再分析資料, 從中獲得大氣溫度垂直廓線數據, 用于風矢量的高度指定。

1.3 風矢量驗證數據

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

2 研究方法

本文提出基于SURF算法的云導風計算方法, 將圖像增強和特征點匹配技術相結合, 運用在時序衛星云圖上, 通過紅外波段云圖反演出風矢量, 具體技術路線如圖1所示。技術流程圖中的關鍵技術方法概要介紹如下。實驗中使用的SURF算法、匹配算法以及RANSAC算法通過Python語言調用OpenCV開源庫中的相應方法實現。

2.1 圖像增強

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

式中, 輸入灰度值Iin歸一化至[0, 1]區間; c為灰度縮放系數, 本文取c=255, 使輸出灰度值Iout的取值范圍在0～255之間。γ為伽馬因子, 當γ<1時, 圖像整體變亮; 反之圖像整體變暗。經試驗, 取γ=5可取得好的圖像增強效果, 具體對比如圖2a、2b所示。可以看到, 經伽馬變換后, 臺風外圍云區中的圖像細節被壓抑, 而密集云區中的圖像對比度增加, 細節得到增強。

a. 臺風云圖; b. 圖像增強后的臺風云圖

2.2 特征點提取

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

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

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

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

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

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

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

2.3 特征點匹配

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

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

2.4 誤匹配過濾

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

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

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

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

4) 將正確的變換矩陣下的內點保留, 外點去除, 實現對特征匹配結果的優化, 完成特征點誤匹配過濾。

2.5 風矢量高度指定

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

1) 根據風矢量的坐標, 在大氣溫度垂直廓線數據中提取該位置的26個高度層對應的大氣溫度集合emp = {1,2,3, …,26}, 其中1>2>3>…>26。

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

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

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

2.6 風矢量質量控制

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

式中,、為目標風矢量的風速和風向;v、d為鄰域風矢量的風速和風向;為鄰域內除目標風矢量以外的風矢量個數。

3 結果

3.1 臺風“山竹”的云導風矢量結果

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

3.2 風矢量精度驗證

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

注: a. 風速差絕對值; b. 風向差絕對值

3.3 風矢量空間分布對比

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

4 討論

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

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

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

4.2 特征點提取方法對云導風結果的影響

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

表2 SIFT云導風方法風矢量精度驗證

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

5 結論

本文提出了一種基于快速魯棒特征(SURF)圖像匹配的云導風計算方法, 在運用伽馬變換增強云區圖像對比度的基礎上, 利用SURF算法計算兩景相鄰時相臺風紅外波段影像的特征點, 以歐氏距離最小的原則完成相似特征點的匹配; 進一步引入RANSAC算法去除誤匹配點對, 根據每組匹配中兩點之間的位移計算出風速的大小和方向, 結合當地大氣溫度廓線指定云高, 經質量控制后得到最終的風矢量。以臺風“山竹”為試驗, 得出的研究結果表明:

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

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

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

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

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

[1] 何斌, 翟國慶, 高坤, 等. 云跡風資料同化對東亞海域風場數值模擬的影響[J]. 海洋學報(中文版), 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] 李華宏, 王曼, 薛紀善, 等. FY-2C云跡風資料在中尺度數值模式中的應用研究[J]. 氣象學報, 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云跡風資料同化應用對臺風預報的影響試驗研究[J]. 大氣科學, 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] 周兵, 徐海明, 吳國雄, 等. 云跡風資料同化對暴雨預報影響的數值模擬[J]. 氣象學報, 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] 許健民, 張其松. 衛星風推導和應用綜述[J]. 應用氣象學報, 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. 基于傅立葉相位分析的衛星云圖導風技術[J]. 氣象科學, 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] 許建明, 王振會. 傅立葉相位分析導風技術的適用范圍和相位重疊分析[J]. 氣象科學, 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] 馬俠霖, 羅鵬, 陳志斌, 等. 基于尺度不變特征變換的云導風計算方法[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] 黃春鳳, 劉守山, 別治峰, 等. 改進的SURF算法在圖像匹配中的應用[J]. 現代電子技術, 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] 唐穎復, 王忠靜, 張子雄. 基于改進SIFT和SURF算法的沙丘圖像配準[J]. 清華大學學報(自然科學版), 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] 權巍, 包鐵壯, 白寶興, 等. SURF與RANSAC算法結合的圖像跟蹤方法[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號臺風“山竹”引起浙江東北部大暴雨成因分析[J]. 海洋預報, 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衛星的逐時次海表溫度融合[J]. 海洋學報, 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]. 南京: 南京氣象學院, 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

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

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

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

(本文編輯: 康亦兼)