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

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

0 引言

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

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

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

1 數(shù)據(jù)來(lái)源

1.1 GF-3 SAR圖像

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

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

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

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

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

2 反演方法

2.1 風(fēng)向求解

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

2.1.1 計(jì)算灰度梯度值

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

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

圖1 Sobel算子

2.1.2 計(jì)算灰度梯度方向值

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

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

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

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

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

式中，Φ為雷達(dá)觀測(cè)方位角，正負(fù)號(hào)由雷達(dá)觀測(cè)方向而定，當(dāng)雷達(dá)觀測(cè)方向?yàn)樽笠晻r(shí)選擇正號(hào)，為右視時(shí)選擇負(fù)號(hào)。

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

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

2.2 風(fēng)速求解

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

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

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

圖2 擴(kuò)展后的CMOD5散射模型曲線圖

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

2.3 反演質(zhì)量評(píng)估

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

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

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果

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

現(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)向基本一致，這說(shuō)明該區(qū)域的ECMWF 風(fēng)向非常接近于實(shí)際風(fēng)向，可用于評(píng)估反演后的SAR 風(fēng)場(chǎng)質(zhì)量。從子圖反演結(jié)果來(lái)看，SAR 風(fēng)向非常接近于ECMWF 方向，可見(jiàn)該子圖風(fēng)場(chǎng)反演效果很好。

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

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

圖4（a）為整個(gè)SAR 圖像反演后風(fēng)場(chǎng)，圖中粗箭頭為各子圖反演得到的SAR 風(fēng)向值，細(xì)箭頭為插值后風(fēng)場(chǎng)，右上角最粗箭頭為浮標(biāo)風(fēng)向，背景灰度值表示不同的風(fēng)速值。圖4（b）為對(duì)應(yīng)的ECMWF風(fēng)場(chǎng)圖，細(xì)箭頭為ECMWF風(fēng)場(chǎng)。為了方便對(duì)比，將圖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)計(jì)，與ECMWF 風(fēng)場(chǎng)相比，SAR 風(fēng)速、風(fēng)向均方根誤差分別為1.79 m/s 和11.95°；與浮標(biāo)最近的子圖SAR 風(fēng)向跟浮標(biāo)相比，SAR 風(fēng)速、風(fēng)向絕對(duì)偏差分別為0.24 m/s和9.03°。

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

4 結(jié)束語(yǔ)

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