光流算法在FY- 4A紅外圖像外推中的精度評價研究

石玉立，施聲偉

(南京信息工程大學 遙感與測繪工程學院， 南京 210044)

衛(wèi)星云圖為天氣預報提供云參數(shù)、大氣流場和各種大氣物理過程等重要的氣象信息，能夠直接明了地看出云系的發(fā)展和形態(tài)變化。彌補了海洋、高原和沙漠地區(qū)氣象觀測點少的不足，還能直觀的監(jiān)測到各種天氣系統(tǒng)的變化。很多氣象領域研究者利用衛(wèi)星產(chǎn)品數(shù)據(jù)研發(fā)定點、定量、定時的預報方法，尤其是使用衛(wèi)星資料和雷達產(chǎn)品相結合進行研究來提升預報水平，獲得了大量研究成果。然而當前氣象臺站對衛(wèi)星資料的應用水平還不高，衛(wèi)星云圖僅作為日常天氣分析和數(shù)值預報輔助工具，甚至被忽視。針對衛(wèi)星云圖圖像外推研究卻很少，云圖的圖像外推工作可提供氣象臺和氣象領域研究者所需要的下個時段云圖發(fā)展的資料，具有很高的利用價值。以往的數(shù)值天氣預報最新進展使我們能夠以很高的分辨率預報大氣動力學[1]，但對于更新周期頻繁的預報應用要求來說，計算成本通常很高，效率很低。因此我們需要進一步開發(fā)計算時間短，精度高的算法。光流法作為一種新的預報算法，外推時間較短、成本較低。光流法是假設云圖亮度強度保持恒定，從一系列衛(wèi)星云圖中跟蹤云圖特征點運動，基于該運動將特征點移動到未來即將到來的位置從而達到外推的目的。

根據(jù)對云圖特征的假設和文獻[2]提出的問題——外推方法擴展到氣象衛(wèi)星圖像的亮度溫度。將衛(wèi)星云圖外推大致分為兩大類，第一類是歐拉持續(xù)性，只需將當前觀測值等同于下時刻衛(wèi)星云圖數(shù)據(jù)，不用推斷系統(tǒng)在空間和時間上的演變。第二類是拉格朗日持續(xù)性，在假定云圖特征亮度強度和運動場是持續(xù)性的前提下對最近的云圖特征進行外推[2-3]。此外，我們還可以根據(jù)預測不確定性對現(xiàn)有的預測方法進行分類：與確定性方法相反，整體臨近預報試圖通過運動場的不同實現(xiàn)和云圖特征亮度強度本身的演變來說明預測的不確定性[4]。在本研究中，我們將模型開發(fā)的重點放在提供確定性的拉格朗日持續(xù)性模型組上。

拉格朗日方法包括2個步驟：跟蹤和預測(外推)[5]。在跟蹤步驟中，從一系列連續(xù)的衛(wèi)星云圖中以每個像素為基礎計算速度場[2-6]。在外推步驟中，使用該速度場平移最近時刻的云圖特征，即根據(jù)云圖的運動場將其移動到即將到來的位置，該步驟基于半拉格朗日方案[2]、插值算法[6]、基于網(wǎng)格的模型實現(xiàn)[4]。可以在每個步驟(跟蹤和預測)中使用不同的算法，從而實現(xiàn)云圖外推過程。

“光流”是很杰出的跟蹤算法，光流通常被理解為從連續(xù)圖像幀中推斷運動模式或速度場的一組方法[6-7]。對于速度場估計，假定亮度是恒定的，并且附加一組光流約束(OFC)條件。OFC的空間屬性標志著光流模型分為兩個主要類別：局部(微分)和全局(變分)[6-8]。局部模型嘗試僅在某些鄰域中設置OFC，而全局模型將OFC應用于整幅圖像。其中還有一組獨特的光譜方法，是將傅立葉變換應用于輸入，并在光譜(傅立葉)域中解析OFC[9]。Bowler等首次將局部光流算法應用到預報中，為模型的發(fā)展提供了新的方向[7]，Liu等提出對衛(wèi)星圖像的每個像素獨立使用局部Lucas-Kanade光流算法[6]，發(fā)現(xiàn)在紅外圖像外推中優(yōu)于全局的HS光流算法[10]，Yeung等使用了不同的全局光流算法在香港建立了業(yè)務預報的SWIRLS產(chǎn)品[8]。

在過去幾十年里，光流算法在臨近預報中應用廣泛，很多研究者把研究重點放在降水、雷暴等研究上，而忽視了大范圍的氣象衛(wèi)星云圖外推研究。本研究中我們基于用于跟蹤步驟的光流公式Sparse(稀疏光流模型)和Dense(稠密光流模型)，以及圖像仿射變換和空間插值的簡化外推技術，通過調(diào)整其參數(shù)開發(fā)了衛(wèi)星云圖的外推程序，比較這2組(每組2種)光流模型在FY- 4A紅外圖像的外推精度，采用表現(xiàn)較好的光流組進行后續(xù)的預報工作。

1 研究數(shù)據(jù)

FY- 4A是我國新一代靜止氣象衛(wèi)星，裝載多種觀測儀器，包括多通道掃描成像輻射計(AGRI)、干涉式大氣垂直探測儀、閃電成像儀和空間環(huán)境監(jiān)測儀器等。其中AGRI包含14個通道。本文研究數(shù)據(jù)是2018年FY- 4A多通道掃描成像輻射計所探測的4 km分辨率全圓盤數(shù)據(jù)，時間間隔為1 h，讀取13通道的紅外圖像數(shù)據(jù)(地理范圍為：2°30′N～62°30′N，67°30′E～150°E)進行圖像外推工作。

衛(wèi)星觀測數(shù)據(jù)和地理信息數(shù)據(jù)是分別存放的，本文采用公式法將FY- 4A標稱上的經(jīng)緯度與行列號進行轉換。將定標級數(shù)與定標表中每一級觀測數(shù)據(jù)一一對應能讀取數(shù)據(jù)[11]。從指定地理范圍讀取相應通道數(shù)據(jù)存儲于1 201×1 651矩陣中，F(xiàn)Y- 4A紅外衛(wèi)星云圖示例如圖1所示。

圖1 FY- 4A紅外衛(wèi)星云圖

在實驗之前，對讀取的每組數(shù)據(jù)依次按照[T-23，T-22，…，T-1，T，T+1，…，T+12]重新命名編排，每組數(shù)據(jù)中選用前24個時間間隔為1 h的全圓盤數(shù)據(jù)作為紅外圖像外推的實驗數(shù)據(jù)，每組后12個數(shù)據(jù)做驗證數(shù)據(jù)(實測數(shù)據(jù))，外推出的圖像是12幅時間間隔為1 h的圖像。選取6組不同日期的實驗數(shù)據(jù)用于本次研究，每組數(shù)據(jù)是36個時間間隔為1 h FY- 4A數(shù)據(jù)，實驗數(shù)據(jù)具體如表1所示。

表1 FY- 4A基于光流模型圖像外推實驗組數(shù)據(jù)

2 光流模型

2.1 光流法原理介紹

2.1.1LK光流法

光流法本質(zhì)上是將真實世界的物體運動轉換成相機平面的二維圖像，并求解目標在圖像上的運動參數(shù)u、v的方法[12]。Lucas-Kanada稀疏光流法是一種2幀差分的光流估計算法，計算得到的一種稀疏光流場[13]。其假設相鄰幀之間圖像亮度不變；圖像中的運動隨時間變化較小；像素點與領域內(nèi)的點具有相似的運動，求解得到光流方程為：

Ixu+Iyv+It=0

(1)

式(1)中：I是整幅圖像；Ix、Iy、It分別是對x、y、t方向求導；u、v分別為x、y方向的速度。

2.1.2DIS光流法

Dense Inverse Search(DIS)光流法關注時間復雜度，提出了一個非常低的時間復雜度和具有競爭力的精度解決方案用于計算稠密光流[17]，它包括2個部分：

1) 快速逆向搜索。在此認為W(X;u)是在像素X上由X=(x，y)T參數(shù)化決定的圖像變換，使得W(X;u)=(x+u，y+u)，逆搜索的目標函數(shù)則為

∑x[It+1(W(X;u))-T(X)]2

(2)

It+1是參考圖像，T(X)是給定的圖像塊。通過迭代最小化找到圖像變換參數(shù)，將圖像變換參數(shù)更新為u←u+Δu得到:

∑x[It+1(W(X;u+Δu))-T(X)]2

(3)

使用高斯-牛頓梯度下降法，將式(3)一階泰勒展開得到：

(4)

對于參數(shù)更新為Δu，則存在一最小二乘法的閉合形式解。將式(3)的偏導數(shù)設為0得到：

(5)

Δu=H-1∑xST·[T(X)-It+1(W(X;u))]

(6)

其中H=∑xSTS是一個n×n大小近似于海森(Hessian)矩陣。

由于S依賴于圖像梯度，因此在每次迭代時必須重新計算位移u、S、海森矩陣H，所以逆向搜索算法是將參考圖像與給定圖像塊的角色顛倒一下，從而避免每次迭代中重新計算Jacobiao和Hessian，目標函數(shù)如下：

∑x[T(W(X;Δu))-It+1(W(X;u))]2

(7)

式(7)一階泰勒展開得到：

(8)

式(8)中：W(X;0)是初始的圖像變換，▽T是W(X;0)在的圖像梯度。使用最小二乘法得到：

(9)

(10)

歸一化后得到：

(11)

2.2 稀疏光流模型組

這組方法的核心思想是識別FY- 4A紅外云圖中用于跟蹤的典型特征，這些特征是衛(wèi)星圖像中云特征強度梯度變化顯著的點，俗稱”角點”。該方法不依賴于尺度大小，比將風暴看作連續(xù)對象進行跟蹤經(jīng)典算法更具有普遍性[14]，因為無需指定任意和尺度相關的云圖特征，而角點的識別僅取決于單元領域內(nèi)梯度的變化。在這一組中，使用了2個在跟蹤和外推方面略有不同的模型。

2.2.1SparseSD光流模型

SparseSD(Sparse Single Delta)讀取兩幅連續(xù)時刻的FY- 4A紅外云圖進行識別、跟蹤和外推圖像特征，假設T時刻是FY- 4A圖像起始時間，模型參數(shù)如表2所示，具體步驟如下：

1) 使用Shi-Tomasi角點檢測器[15]識別T-1時刻衛(wèi)星云圖中的特征。該檢測器基于角質(zhì)量度量公式min(λ1，λ2)來確定圖像中強角點。原理是Shi-Tomasi在Harris檢測算法上進行了改進，其自相關矩陣的2個特征值中較小的一個大于設定的閾值，λ1、λ2是自相關矩陣的特征值，則可以得到強角點。

2) 使用局部Lucas-Kanade光流算法[16]跟蹤T時刻衛(wèi)星云圖上特征點(T-1時刻識別出的)并計算出大氣運動矢量(AMV)。該算法使用最小二乘法求解局部特征鄰域中的一組光流方程，從而在T時刻識別我們先前在T-1時刻FY- 4A紅外圖像中識別出的特征的位置。

3) 使用線性外推法外推特征點運動，以便預測出在每個步長時間為n的特征點位置。

4) 使用最小二乘法，根據(jù)T到T+n時刻紅外圖像所有已識別特征點的位置，計算每個時間步長n的仿射變換矩陣[17]。

5) 通過仿射變換外推T時刻處的FY- 4A紅外圖像，具體過程是將T時刻的FY- 4A紅外云圖每個像素點的位置通過仿射變換矩陣轉換為將來T+1～T+n時刻的紅外云圖像素點位置。對預測圖像中的剩余像素值進行線性插值，獲得與T時刻相對應的FY- 4A紅外圖像云圖特征強度圖[18]。

2.2.2Sparse光流模型

Sparse模型讀取24幅時間間隔一小時的FY- 4A紅外圖像，僅考慮在整個時間段(24個時間步長)中保持不變的特征。模型參數(shù)如表2所示，具體步驟如下：

1) 使用Shi-Tomasi角點檢測器識別T-23時刻紅外圖像上的特征點[15]。

2) 使用局部Lucas-Kanade光流算法跟蹤T-22到T時刻紅外圖像上這些特征點(T-23時刻識別出的)[16]。

3) 依據(jù)時間(T-23到T)為每個成功跟蹤的特征獨立地建立線性回歸模型。

4) 同SparseSD外推的3-5步驟。

2.3 稠密光流模型組

稠密模型組使用Dense Inverse Search(DIS)全局光流算法[17]。該算法通過對兩幅連續(xù)的FY- 4A紅外圖像的分析估計每個圖像像素點的移動速度。DIS算法與其他全局光流算法(如DeepFlow、PCAFlow)相比有更高的準確性和更高的計算效率[18-19]。

該組中的兩個模型僅在外推(或平流)步驟方有所不同。第1個模型(Dense)使用恒定矢量平流方案[7]，而第2個模型(DenseRotation)使用半拉格朗日平流方案[2]。2種方法之間的主要區(qū)別在于，假定運動場本身是持續(xù)的(如圖2所示)，恒定矢量方案不允許旋轉運動。半拉格朗日方案允許表示大尺度旋轉運動。

圖2 4種平流方案的位移向量圖，P是一個地塊的網(wǎng)格點

2種平流方案的實施方式有2種可能的選擇：時間上向前(在空間的下游)或時間上向后(在空間的上游)(如圖2所示)。在向前方案中，從網(wǎng)格點P開始，將地塊向下游平流直至點Q。因此，“向前”表示在時間上向前，在空間上在下游。而在向后方案中，向上游移動并確定地塊的原點O，該原點將在網(wǎng)格點P處結束。通常，O和Q都不完全與網(wǎng)格點重合。在向前方案中，平流值被重新分配到相鄰的網(wǎng)格點，而在向后方案中，需要插值方法來確定點O處的值。目前尚不清楚哪種方案可以被視為基于FY- 4A數(shù)據(jù)圖像外推的最合適和通用的解決方案，一方面是關于質(zhì)量守恒，另一方面是小尺度的能量損失。因此，進行了實驗，將向前與向后，恒定矢量與半拉格朗日平流進行了任何可能的組合。將向后方案用作Dense和DenseRotation模型的參數(shù)。

Dense和DenseRotation模型都使用線性插值方法，以便將預測位置的衛(wèi)星云圖像素值插值到對應的FY- 4A衛(wèi)星云圖位置。稠密組模型的模型參數(shù)如表2所示。

表2 FY- 4A基于光流模型圖像外推參數(shù)

具體實現(xiàn)步驟如下：

1) 基于T-1和T時刻的FY- 4A紅外圖像，使用全局DIS光流算法[17]計算速度場。

2) 對于T+n時刻來說根據(jù)速度場使用向后恒定向量[7]或向后半拉格朗日方案[2]逐步外推每個像素值，對于半拉格朗日方案，我們通過將速度場線性插值到該時間步的像素位置來更新每個預測時間步n處位移像素的速度。

3) 平流步驟的結果，獲得了不規(guī)則速度場，由從原始位置偏移的原圖像像素點組成的。使用在T+n時刻預測位置處每個位移像素的強度，對原始FY- 4A紅外圖像每個網(wǎng)格點處的像素強度進行反距離加權插值[6]。在插值之前對時間步長內(nèi)每個像素進行平流處理從而將數(shù)值擴散降到最低(如Germann和Zawadzki在2002年提出的”一次插值”方法)[2]。這樣，以避免插值效應使空間中的圖像特征變得平滑。

2.4 持續(xù)性(Persistence)基準模型

歐拉持續(xù)性基準模型是一種啟發(fā)式方法，假設云圖運動場在時間和空間上是不變的，即對于任何時間步長n，云圖運動場與T時刻運動場相同。這種預測技術比較粗糙，但很簡單，對于很短的步長時間卻是一個強大的模型，這使得歐拉持續(xù)性模型在現(xiàn)實世界的應用程序中是一個易于使用的選項。在目前的分析背景下，歐拉持續(xù)性型是一個基線，用來評估增加的復雜性水平是否可以提高預測技能。同時，它可以很好地衡量不同事件在時間上的去相關性。

3 光流模型評估

為了進行驗證模型效果，使用4種類型的評價方式：平均絕對誤差(MAE)用以評價外推圖像與實際圖像之間的差異，值越小差異越小。均方根誤差(RMSE)用以評價預測值與真實值之間的差異，值越小模型精度越高。平均絕對百分比誤差(MAPE)用以評價預測值與真實值逼近程度，值越小預測值與真實值越接近，模型精度越高。與歐拉基準模型(Persistence)比較，光流模型外推效果整體上得優(yōu)于基準模型。

3.1 評價指標

使用Mean Absolute Error(MAE)用以評價預測值與真實值之間的差異。公式如下：

綜上所述,高血壓患者采取綜合護理干預,臨床價值較高,可有效控制血壓,促使其配合度的提高,利于預后效果改善。

(14)

使用Root Mean Square Error(RMSE)用來衡量預測值與真實值之間的偏差。公式如下：

(15)

使用Mean Absolute Percentage Error(MAPE)用以計算預測值與真實值的逼近程度。公式如下：

(16)

其中simi和obsi分別指外推的衛(wèi)星云圖第i個像素值和實際云圖的第i個像素值，n是像素數(shù)。

3.2 結果分析

研究時間內(nèi)，所有模型(Sparse，SparseSD，Dense，DenseRotation，Persistence)均用程序計算1～12 h(時間間隔為1 h)的12幅外推圖像。選擇時間步長為3、6、9和12 h的外推圖像示例如圖3所示。

從圖3來看，2組模型(每組2種)在外推結果上總體較好，對鋒區(qū)及臺風重疊區(qū)域的位置均能較為準確的推演，隨著時間步長增加，外推出的圖像與實際相比整體圖像亮度有損失，外推云圖與實際相比部分輪廓不清晰，整體上不影響實際工作中對鋒區(qū)和臺風影響范圍的預報。由于衛(wèi)星云圖邊緣特征點不明顯，Sparse模型組在使用Shi-Tomasi角點檢測器難檢測出云圖左上角邊緣特征，并且局部檢測出的特征無法代表整體云圖的運動模式，隨著外推步長時間增加，Sparse模型組外推的圖像邊緣有像素點缺失并隨時間步長增加缺失數(shù)目增多。而Dense模型組外推圖像整體完整，外推云圖里視覺上與實際差距不大，內(nèi)部輪廓不清晰。

圖3 2018年7月20日09∶00 Sparse和Dense模型組圖像外推圖像

圖4表示外推圖像像素值與實際圖像的差，就空間上來看外推效果，Sparse模型組和Dense模型組在圖像y軸600～800區(qū)間和x軸0～500區(qū)間區(qū)域及圖像右下角外推較不好，隨著時間步長增加外推效果下降。由于實際圖像中云團云系類型較多較雜導致兩組模型都外推效果均不佳。Sparse和Dense模型組在步長時間為3時外推效果整體較好。

圖4 2018年7月20日09∶00模型組圖像外推圖與實際圖差

3.3 模型比較

圖5 Sparse模型組和Dense模型組平均絕對誤差(MAE)曲線

圖6 Sparse模型組和Dense模型組均方根誤差(RMSE)曲線

圖7 Sparse模型組和Dense模型組平均絕對百分比誤差(MAPE)曲線

圖5和圖6分別表示平均絕對誤差(MAE)，均方根誤差(RMSE)。其表示外推衛(wèi)星云圖圖像像素值與實際的差異，MAE和RMSE越小說明與實際圖像差異越小，從而外推精度越高。從圖6中看出：Dense模型組與歐拉持續(xù)基準對比明顯優(yōu)于Sparse模型組。

圖7是MAPE，從圖7中可以看出2組光流模型總體精度在11%之內(nèi)，整體精度較高，Dense組2種模型精度優(yōu)于Sparse組2種模型。

從誤差圖7中可以看出Dense模型組在歐拉持續(xù)性方面明顯優(yōu)于Sparse模型組，Sparse模型組所使用的Shi-Tomasi角點角檢測器識別出的局部特征(角點)可能無法代表云圖的整體運動：該檢測器集中于云圖特征有明顯的強度和梯度變化的特征點，其運動可能不代表占主導的中尺度運動模式。Dense模型組將圖像中每個像素與位移關聯(lián)能夠更好的外推出云圖圖像，Dense模型和DenseRotation之間的差異(換句話說是恒定向量和半拉格朗日方案之間的差異)可以忽略不計。半拉格朗日方案的理論優(yōu)越性可能會通過持續(xù)旋轉運動在其他事件(臺風、氣旋等預測)中體現(xiàn)出來。因此，應在以后的研究中進行更全面的分析。Dense和DenseRotation具有相同的插值過程，插值包括圖像變形的后處理(稀疏模型)和網(wǎng)格化臨近預報的計算中(Dense模型的一部分)。通常，這樣的插值步驟會導致數(shù)值擴散，從而導致小尺度特征強度的退化或損失[2]。然而，通過在特定時間步長對預測都執(zhí)行插值，我們就能夠?qū)ο∈枘Ｐ秃统砻苣Ｐ徒M都包含此類影響。正如文獻[2]所說，對于較長的步長時間，這些影響可能很明顯，取決于所使用的外推技術。

對于圖像外推而言，模型計算性能可能是滿足頻繁更新的用戶需求的一個重要標準。在具有Intel(R)Core(TM)i7-9750H CPU(16 G內(nèi)核，2.6 GHz)的筆記本上使用Jupyter notebook運行的圖像外推模型程序，得到稀疏組生成12幅外推圖像(1 h的時間步長)的平均時間是3.93～7.28 s，稠密組為13.6～19.1 s。因此Dense模型組的計算成本高。

4 結論

對Sparse和Dense光流模型組每組兩種模型在FY- 4A紅外圖像外推中應用與歐拉持續(xù)性基準模型對比進行了精度評價。結果顯示光流模型比以往的數(shù)值預報成本低，計算效率高，Dense模型組采用的DIS光流算法和向后平流方案外推效果優(yōu)于Sparse模型組采用的Shi-Tomasi角點檢測和LK光流算法方案。Dense模型組通過將圖像中的每個像素點與前后兩幀的圖像的位移關聯(lián)起來，從而外推效果較好，Dense模型組比歐拉基準模型效果更好。更能夠為缺少氣象站點的區(qū)域提供下個時段云圖發(fā)展資料。

從平均絕對誤差、均方根誤差、平均絕對百分比誤差評價指標上看，Dense模型組明顯優(yōu)于Sparse模型組，兩組模型各有其優(yōu)點，Sparse模型組在局部外推較好，而邊緣有信息缺失，計算時間短。Dense模型組外推較為完整，計算時間長。因此綜合考慮在對FY- 4A紅外圖像進行圖像外推工作時，優(yōu)先考慮Dense模型組。