基于梯度算子的高分三號SAR風(fēng)場反演方法

陳克海 張金蘭 解學(xué)通 鄧 丹

（1.廣東工貿(mào)職業(yè)技術(shù)學(xué)院 測繪遙感信息學(xué)院，廣東 廣州 510510；2.廣州大學(xué) 地理科學(xué)與遙感學(xué)院，廣東 廣州 510006；3.廣東國圖勘測地理信息有限公司，廣東 廣州 510180）

0 引言

目前海面風(fēng)場遙感監(jiān)測傳感器主要有散射計、高度計和合成孔徑雷達（synthetic aperture radar，SAR）。其中，高度計垂直向下觀測，能觀測星下線附近區(qū)域的風(fēng)速，卻無法觀測風(fēng)向值。散射計一般具有兩個旋轉(zhuǎn)的波束，這兩個波束的極化方式和入射角不同，在較大范圍內(nèi)對地重復(fù)觀測，在同一地面單元上可產(chǎn)生多個觀測值，通過反演算法反演出風(fēng)速、風(fēng)向值。但是，散射計一般在25 km×25 km 的范圍內(nèi)只有一個風(fēng)速、風(fēng)向反演值，其風(fēng)場數(shù)據(jù)僅適用于全球環(huán)境研究與應(yīng)用［1-2］。SAR 空間分辨率高，最高可達到米級，能觀測到海面上中小尺度風(fēng)場結(jié)構(gòu)，為區(qū)域環(huán)境研究和應(yīng)用提供數(shù)據(jù)，是目前海面上中小尺度風(fēng)場觀測的熱點方向之一［3］。

SAR 能捕捉到海面風(fēng)場引起的黑白相間的風(fēng)條紋，由于風(fēng)條紋走向基本平行于海面風(fēng)向，可以利用風(fēng)條紋信息從SAR 圖像中提取海面風(fēng)向［4-5］。目前SAR 圖像風(fēng)場反演一般分兩步走，第一步根據(jù)SAR 圖像風(fēng)條紋走向求解出風(fēng)向值，第二步在已知的后向散射模型函數(shù)條件下，根據(jù)圖像后向散射系數(shù)和已求解得到的風(fēng)向值推算出風(fēng)速值。可見，SAR 圖像風(fēng)場反演關(guān)鍵在于風(fēng)向值求解。目前，一般采用頻譜分析方法求解SAR 圖像風(fēng)向值。首先利用傅里葉變換等方法將SAR 圖像從空間域轉(zhuǎn)換到頻率域，然后計算功率譜，并基于功率譜上波峰走向垂直于風(fēng)條紋走向的假設(shè)，從波峰走向推算出海面風(fēng)向［6-8］。

頻譜分析方法具有較高的數(shù)學(xué)理論依據(jù)，但是存在一些缺點，比如頻譜分析法要求圖像中海洋區(qū)域必須完整，區(qū)域內(nèi)不能出現(xiàn)陸地、海島或者船只等其他物體，否則無法開展傅里葉變換及后續(xù)處理。另外，頻譜分析法要求圖像區(qū)域要足夠大，否則會降低精度。借鑒在圖像邊緣檢測技術(shù)中被普遍采用的灰度梯度算子在提取邊緣方向時對處理區(qū)域的形狀和大小沒有特別要求，在SAR 圖像風(fēng)場反演算法中引入灰度梯度算子，用于計算圖像各像素灰度值梯度方向。在統(tǒng)計意義上，風(fēng)條紋走向應(yīng)垂直于分布概率最大的梯度方向。據(jù)此，可以利用梯度算子推算出海面風(fēng)向。為了測試本文方法，對高分三號（GF-3）SAR圖像進行風(fēng)場反演實驗，反演得到的SAR 風(fēng)場與歐洲中尺度天氣預(yù)報中心（European Centre for Medium-Range Weather Forecasts，ECMWF）再分析風(fēng)場數(shù)據(jù)和美國國家浮標(biāo)中心（National Data Buoy Center，NDBC）浮標(biāo)風(fēng)矢量進行對比，分析SAR風(fēng)場質(zhì)量 。

1 數(shù)據(jù)來源

1.1 GF-3 SAR圖像

我國GF-3 衛(wèi)星于2016 年8 月10 日成功發(fā)射［9］，其上搭載的SAR 工作在C 波段。它具有多個成像模式，本實驗采用其中的標(biāo)準(zhǔn)條帶成像模式，數(shù)據(jù)等級為L1A，文件保存為GeoTiff 格式。每幅圖像約覆蓋10 000 km2，空間分辨率約為5 m×7 m。在L1A 文件中，像素行號沿著衛(wèi)星前進方向遞增，而列號由近端向遠(yuǎn)端遞增。每個L1A 文件元數(shù)據(jù)文件包含定位、觀測參數(shù)等基礎(chǔ)數(shù)據(jù)，利用這些基礎(chǔ)數(shù)據(jù)可以進行輻射校正、空間定位和觀測幾何參數(shù)計算等預(yù)處理。

1.2 輔助數(shù)據(jù)

本文使用ECMWF 再分析風(fēng)場數(shù)據(jù)和NDBC浮標(biāo)風(fēng)矢量等輔助數(shù)據(jù)來評估風(fēng)場反演精度。

（1）ECMWF 再分析風(fēng)場數(shù)據(jù)是歐洲空間局事后整合全球衛(wèi)星數(shù)據(jù)和各種實測數(shù)據(jù)優(yōu)化得到的高精度數(shù)值模擬產(chǎn)品，特點是數(shù)據(jù)量多、精度較高，被廣泛用于海洋衛(wèi)星產(chǎn)品質(zhì)量評估［10-11］。選用的ECMWF 再分析風(fēng)場的空間、時間分辨率分別為0.125°×0.125°和6 h，通過在空間維度上采用二維線性插值和在時間維度上采用三次樣條函數(shù)插值，可得到任意時刻、位置的風(fēng)速和風(fēng)向值。

（2）NDBC 浮標(biāo)風(fēng)矢量數(shù)據(jù)由美國國家浮標(biāo)中心提供［12］。被選用的浮標(biāo)位于離岸50 km以上的海域，具有連續(xù)的觀測能力，每10 min 記錄一次風(fēng)速、風(fēng)向測量值，通過插值得到任意時刻的風(fēng)速、風(fēng)向值。

2 反演方法

2.1 風(fēng)向求解

在開始風(fēng)場反演之前，以100 m×100 m 的分辨率對SAR 圖像重新進行采樣，并將圖像劃分為25 km×25 km 的子圖［13］。風(fēng)場反演技術(shù)是在子圖上進行。

2.1.1 計算灰度梯度值

子圖各個像素在水平、垂直方向上的灰度梯度值Gx和Gy為

式中，σoij為窗口內(nèi)第i行、第j列的后向散射系數(shù)；Sxij和Syij為分別為在水平、垂直方向上的Sobel算子值，具體數(shù)值見圖1。在圖像空間上移動Sobel 算子，計算子圖各個像素在水平、垂直方向上的灰度梯度值Gx和Gy。

圖1 Sobel算子

2.1.2 計算灰度梯度方向值

根據(jù)梯度值Gx和Gy，計算子圖每個像素灰度梯度方向值D為

2.1.3 計算風(fēng)向可能解

根據(jù)各像素灰度梯度方向值，繪制灰度梯度方向概率密度曲線。由于噪聲的存在，該曲線存在較大起伏，使用均值濾波器進行平滑處理。對于風(fēng)條紋而言，風(fēng)條紋走向的灰度梯度方向必然對應(yīng)于概率密度曲線峰值，由此可以推算出海面風(fēng)向?為

式中，φ為平滑處理后概率密度曲線峰值對應(yīng)的灰度梯度方向值。由式（4）可知，由灰度梯度值計算得到的風(fēng)向不可避免地存在180°模糊，即存在兩個可能風(fēng)向，分別為?和?+180°。

以上海面風(fēng)向是在圖像坐標(biāo)系上計算得到的，需轉(zhuǎn)化到地理坐標(biāo)系上。對于GF-3 SAR，觀測方向垂直于衛(wèi)星方向，列號沿著觀測方向遞增，可見圖像坐標(biāo)系x軸正方向為觀測方向；行號沿著衛(wèi)星前進方向遞增，可見y軸正方向為衛(wèi)星前進方向。根據(jù)觀測幾何關(guān)系，SAR 風(fēng)向?可從圖像坐標(biāo)系轉(zhuǎn)換到地理坐標(biāo)系，公式為

式中，Φ為雷達觀測方位角，正負(fù)號由雷達觀測方向而定，當(dāng)雷達觀測方向為左視時選擇正號，為右視時選擇負(fù)號。

2.1.4 去除風(fēng)向模糊

根據(jù)ECMWF 再分析風(fēng)場數(shù)據(jù)，從兩個模糊解中挑選出與ECMWF 風(fēng)向最接近的一個作為最終風(fēng)向解。

2.2 風(fēng)速求解

GF-3 SAR 工作在C 波段。CMOD5 為目前比較成熟的C 波段散射模型，但是只有垂直（VV）極化模型，缺少水平（HH）極化模型［14］。對于重置（HH）極化模型值，可使用極化比系數(shù)由VV 極化模型值轉(zhuǎn)化得到［15］。

式中，σH0H和σV0V分別為HH、VV 的后向散射系數(shù)模型值；u和?分別為風(fēng)速、風(fēng)向值；θ為雷達入射角；α為極化比系數(shù)，這里取值1。

圖2為擴展后的CMOD5散射模型曲線圖（入射角為40°），其中圖2（a）、2（b）分別對應(yīng)垂直（VV）和HH 極化，圖中各曲線從下到上分別表示風(fēng)速為5、10、15、20、25、30 m/s 時不同相對風(fēng)向的后向散射系數(shù)模型值。

圖2 擴展后的CMOD5散射模型曲線圖

從子圖中計算出風(fēng)向之后，利用擴展后的CMOD5散射模型，搜索出與該區(qū)域后向散射系數(shù)測量值均值最接近的風(fēng)速值，作為該子圖的風(fēng)速值。

2.3 反演質(zhì)量評估

使用ECMWF 再風(fēng)場數(shù)據(jù)和NDBC 浮標(biāo)風(fēng)矢量數(shù)據(jù)，對反演得到的SAR 風(fēng)場質(zhì)量進行評估。對于ECMWF 再風(fēng)場數(shù)據(jù)，以風(fēng)速、風(fēng)向均方根誤差為評估參數(shù)；對于NDBC浮標(biāo)風(fēng)矢量數(shù)據(jù)，以風(fēng)速、風(fēng)向絕對偏差為評估參數(shù)。

設(shè)在第i個匹配數(shù)據(jù)中，SAR 風(fēng)速、風(fēng)向分別為si和di，ECMWF或NDBC風(fēng)速、風(fēng)向分別為si0和di0，則以上評估參數(shù)可表示為：

（1）風(fēng)速絕對偏差Esi

（2）風(fēng)向絕對偏差Edi

（3）風(fēng)速均方根誤差Rs

（4）風(fēng)向均方根誤差Rd

以上風(fēng)速、風(fēng)向絕對偏差反映SAR 風(fēng)速、風(fēng)向值與參考值的絕對偏差值；風(fēng)向、風(fēng)速均方根誤差綜合反映SAR 風(fēng)速、風(fēng)向值與參考值之間的偏差，是衡量反演精度的主要參考指標(biāo)之一。

3 實驗結(jié)果

以產(chǎn)品號為2131355 的GF-3 SAR 圖像為實驗數(shù)據(jù)，先將圖像劃分為25 km×25 km 的子圖，然后使用本文方法進行風(fēng)場反演實驗，反演得到的SAR 風(fēng)場與ECMWF 再分析風(fēng)場和NDBC 浮標(biāo)矢量數(shù)據(jù)進行對比分析。

現(xiàn)在以該SAR 圖像中第2 行第3 列的子圖反演情況為例。圖3（a）為該子圖反演得到的風(fēng)向和相應(yīng)ECMWF 風(fēng)向。圖中灰度圖為SAR 圖像，粗、細(xì)箭頭分別表示SAR 風(fēng)向和ECMWF 風(fēng)向。顯然，該子圖具有明顯的明暗相間的風(fēng)條紋，風(fēng)條紋走向與ECMWF 風(fēng)向基本一致，這說明該區(qū)域的ECMWF 風(fēng)向非常接近于實際風(fēng)向，可用于評估反演后的SAR 風(fēng)場質(zhì)量。從子圖反演結(jié)果來看，SAR 風(fēng)向非常接近于ECMWF 方向，可見該子圖風(fēng)場反演效果很好。

圖3 GF-3 SAR子圖反演結(jié)果示例

為了說明本文算法，這里給出該子圖風(fēng)向計算的中間結(jié)果。圖3（b）為該子圖對應(yīng)的灰度梯度方向概率密度曲線，圖中細(xì)曲線為原始的灰度梯度方向概率密度曲線，抖動較大，經(jīng)平滑處理后為粗曲線。平滑后概率密度曲線峰值對應(yīng)的灰度梯度方向垂直于海面風(fēng)向。

圖4（a）為整個SAR 圖像反演后風(fēng)場，圖中粗箭頭為各子圖反演得到的SAR 風(fēng)向值，細(xì)箭頭為插值后風(fēng)場，右上角最粗箭頭為浮標(biāo)風(fēng)向，背景灰度值表示不同的風(fēng)速值。圖4（b）為對應(yīng)的ECMWF風(fēng)場圖，細(xì)箭頭為ECMWF風(fēng)場。為了方便對比，將圖4（a）SAR 風(fēng)向和浮標(biāo)風(fēng)向也顯示在圖4（b）上。從圖4（b）可以看出，各子圖SAR風(fēng)向非常接近于ECMWF 風(fēng)向和浮標(biāo)風(fēng)向。經(jīng)統(tǒng)計，與ECMWF 風(fēng)場相比，SAR 風(fēng)速、風(fēng)向均方根誤差分別為1.79 m/s 和11.95°；與浮標(biāo)最近的子圖SAR 風(fēng)向跟浮標(biāo)相比，SAR 風(fēng)速、風(fēng)向絕對偏差分別為0.24 m/s和9.03°。

圖4 GF-3 SAR圖像風(fēng)場反演結(jié)果

4 結(jié)束語

本文針對SAR 圖像提出了一種基于灰度梯度算子的風(fēng)場反演方法，并使用GF-3 SAR 圖像進行實驗，反演得到的SAR 風(fēng)場與ECMWF 再分析風(fēng)場對比，風(fēng)速、風(fēng)向均方根誤差分別為1.79 m/s 和11.95°；與NDBC 浮標(biāo)風(fēng)矢量對比，風(fēng)速、風(fēng)向絕對偏差分別為0.24 m/s 和9.03°，可見該方法可以有效從SAR 圖像中反演出較高精度的風(fēng)場。由于缺乏足夠的GF-3 SAR 圖像作為實驗數(shù)據(jù)，本文僅作為SAR 海面風(fēng)場反演的一種探索性研究。另外，由于風(fēng)條紋的產(chǎn)生與海洋大氣邊界層的穩(wěn)定性等因素有關(guān)，并非所有風(fēng)場都能產(chǎn)生風(fēng)條紋，這意味著并非所有SAR 圖像都存在風(fēng)條紋。對于沒有風(fēng)條紋的SAR 圖像，需采用其他方法進行風(fēng)場反演。