適用于Ku波段雷達高度計海冰干舷高度反演的積雪校正方法

張 晰，趙全芳,孟俊敏，劉眉潔

（1.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266590；2.自然資源部第一海洋研究所，山東 青島 266061;3.青島大學(xué) 物理學(xué)院，山東 青島 266071）

海冰直接影響著全球的氣候變化，并通過與大氣的相互作用反饋于全球環(huán)境系統(tǒng)。海冰厚度是重要的氣候環(huán)境變化表現(xiàn)因子。準確估算海冰厚度，對于氣候變化研究、極區(qū)航行保障具有重要意義。利用遙感手段準確估算大范圍海冰厚度始終是當前的國際難題。相比于人工實測、電磁感應(yīng)和仰視聲吶等手段，近年來發(fā)展的雷達高度計測高技術(shù)使得獲取連續(xù)大范圍的海冰厚度成為可能[1]。到目前為止，ERS-1/2、ENVISAT、HY-2、AltiKa、CryoSat-2和Sentinel-3A等雷達高度計相繼發(fā)射，為海冰厚度估算提供了多源的衛(wèi)星觀測手段。在眾多雷達高度計中，工作在Ku波段的CryoSat-2（CS-2）和Sentinel-3A（S3）合成孔徑雷達高度計相較于傳統(tǒng)雷達高度計，具有更高的空間分辨率（沿軌約0.3 km，CS-2/S3交軌約1.5/1.64 km），能夠提供較高分辨率的海冰厚度信息，是當前最先進的衛(wèi)星雷達高度計。

雷達高度計并不直接計算海冰厚度，是通過計算海冰的冰面出水高度（海冰干舷高度），并結(jié)合浮體定律實現(xiàn)海冰厚度估算[2-3]。因此，準確反演海冰干舷高度對于估算海冰厚度至關(guān)重要。由于衛(wèi)星高度計信號回波的實際時間跟蹤點與預(yù)設(shè)時間跟蹤點存在位置上的偏差，需要對回波波形進行重跟蹤處理，以計算實際時間跟蹤點，進而校正得到衛(wèi)星高度計質(zhì)心到地面點的真實距離。Ricker等[4]使用波形重跟蹤方法指出重跟蹤閾值的選擇對海冰干舷和厚度的估計值有顯著影響。在極區(qū)，海冰表面普遍覆蓋著積雪，當雷達高度計探測海冰表面高度時，Ku波段電磁波須穿過覆蓋于海冰之上的積雪層，因此在計算海冰干舷高度時，需考慮電磁波在積雪中穿透速度較慢的距離校正問題。在實際情況下，積雪層受海冰析鹽過程的影響，常會在積雪的上下層形成強鹽度梯度變化[5-7]；另外，積雪層內(nèi)粒徑也會與電磁波發(fā)生體散射[8]。這些因素使得微波雷達高度計的電磁波很難完全穿透積雪層直達海冰表面。對于短波長的Ku波段雷達，該現(xiàn)象更為明顯，這將導(dǎo)致雷達信號主散射面向上偏離雪-冰交界面，從而影響海冰干舷高度的估計。但大多數(shù)研究者在使用雷達高度計反演海冰干舷高度時，常假設(shè)Ku波段電磁波可以完全穿透積雪層，直接將雪-冰交界面當作雷達信號的主散射面[2-4,9-12]，利用海冰上的整體積雪深度進行距離校正。Armitage和Ridout[13]分析了CS-2雷達高度計電磁波在積雪中的雷達穿透系數(shù)，發(fā)現(xiàn)Ku波段電磁波在垂直入射時無法完全穿透積雪層，并指出若不考慮積雪層對電磁波的影響，海冰干舷高度將被高估，從而大大影響海冰厚度的估算。Ricker等[14]發(fā)現(xiàn)，未知雷達信號穿透到積雪層中的不確定性可能導(dǎo)致大約0.06～0.12 m的海冰干舷高度偏差。

因此，選擇合適的冰間水道與海冰的重跟蹤閾值以及確定雷達信號在積雪中的穿透深度，對精確地反演海冰干舷高度起決定性作用?；诖耍疚姆謩e針對CS-2與S3兩種高度計衛(wèi)星開展了確定最優(yōu)波形重跟蹤閾值的分析，并詳細分析了積雪層對Ku波段CS-2和S3合成孔徑雷達高度計主散射面位置的影響，計算了Ku波段電磁波的積雪穿透系數(shù)，針對海冰干舷高度精確反演的需求，改進了積雪校正算法。

1 數(shù)據(jù)

1.1 CryoSat-2雷達高度計數(shù)據(jù)

CryoSat-2是歐洲空間局發(fā)射的合成孔徑雷達高度計，其工作在Ku波段，中心頻率為13.575 GHz，帶寬約為320 MHz，空間覆蓋范圍達到南北緯88°，運行重復(fù)子周期為30 d，其沿軌道分辨率可達0.3 km，交軌道分辨率優(yōu)于1.5 km，相比于足跡為10 km的傳統(tǒng)雷達高度計[15]，其測量精度顯著提升。該衛(wèi)星有LRM（低分辨率）、SAR和SARin 3種工作模式，本文采用的是SAR模式L1b級數(shù)據(jù)。在該模式下，CS-2的每一個波形采樣窗口大小為60 m，包含有256個距離門[16]。為了進行波形重跟蹤校正，得到海冰/水的表面高程，還需對L1b數(shù)據(jù)進行波形噪聲去除和地球物理校正（地球物理校正包含：干/濕對流層、電離層、逆氣壓、海洋潮汐、長周期平衡潮、海洋負荷潮、固體潮和地心極潮校正等處理）。本文使用的是北極區(qū)域、時間為2017年3月-4月和2018年4月的雷達數(shù)據(jù)。

1.2 Sentinel-3A雷達高度計數(shù)據(jù)

Sentinel-3A也是歐空局發(fā)射的合成孔徑雷達高度計，它是第一個100%以SAR模式覆蓋海洋、冰區(qū)和內(nèi)陸水域的雷達高度計。其工作的Ku波段中心頻率與CS-2相同，帶寬約為350 MHz，空間覆蓋范圍為南北緯81.35°，運行重復(fù)子周期為4 d，沿軌道分辨率約為0.3 km，交軌道分辨率大約為1.64 km[17]。本文采用的是SAR模式的L2 WAT Enhanced（Water Enhanced）數(shù)據(jù)。在該模式下波形采樣窗口大小為60 m，包含有128個距離門。與CS-2相同，使用的也是北極區(qū)域2017年3月-4月和2018年4月的數(shù)據(jù)，并對這些數(shù)據(jù)進行了與CS-2相同的地球物理校正處理。

1.3 Operation IceBridge數(shù)據(jù)

本文使用美國國家航空航天局（NASA）發(fā)布的OIB機載數(shù)據(jù)作為海冰干舷高度的驗證數(shù)據(jù)，其每年均會在北極、格陵蘭島和南極上空開展機載探測飛行。該數(shù)據(jù)由機載地形掃描測繪儀（Airborne Topographic Mapper,ATM）、數(shù)字測繪相機（Digital Mapping System Camera,DMS） 和雪雷達（Snow Radar）3種裝備采集得到。ATM用于測量表面高度，空間分辨率為0.4 m，垂直分辨率為0.03 m；DMS用于識別海冰、海水和冰間水道等，空間分辨率為0.1 m；雪雷達用于測量積雪厚度，空間分辨率為40 m，垂直分辨率為0.06 m[18]。綜合上述3種傳感器，機載OIB可提供海冰的激光干舷高度和積雪厚度等數(shù)據(jù)。圖1展示了2017-2018年春季機載OIB北極飛行路線及激光干舷高度分布圖。

對于NASA發(fā)布的機載OIB數(shù)據(jù)，其海冰激光干舷高度、積雪厚度和產(chǎn)品空間分辨率都統(tǒng)一為40 m[18]。在此需要指出的是：機載OIB數(shù)據(jù)攜帶的測高儀器為激光雷達，因此，得到的海冰激光干舷高度并非海冰干舷高度，而是海冰干舷高度與積雪厚度之和。為得到海冰干舷高度，需將OIB得到的激光干舷高度減去雪雷達得到的積雪厚度，才能得到真實的海冰干舷高度。根據(jù)文獻記載，機載OIB的海冰激光干舷高度探測精度為0.015±0.06 m[19]，積雪厚度的探測精度為0.01±0.05 m[20]。本文所用的OIB數(shù)據(jù)均與CS-2和S3雷達高度計數(shù)據(jù)時間、地點范圍相同（北極，2017年3月-4月和2018年4月）。由于CS-2、S3雷達高度計與OIB數(shù)據(jù)的空間分辨率不同，為進行逐一比較，將雷達高度計數(shù)據(jù)與OIB數(shù)據(jù)在25 km的網(wǎng)格中取平均，實現(xiàn)數(shù)據(jù)一一對應(yīng)。

圖1 春季機載OIB北極飛行路線及激光干舷高度分布圖

2 Ku波段雷達高度計波形重跟蹤閾值確定

計算海冰干舷高度可通過反演海冰面到雷達質(zhì)心的高程和冰間水道面到雷達質(zhì)心的高程，并計算二者的差值得到。因此，首先需要區(qū)分雷達高度計的回波波形，以識別海冰、開闊水域和冰間水道等不同地物類型，然后再反演各種地物類型表面的高程。

對于地物類型識別，主要是利用Ocean and Sea Ice SAF（OSI-SAF）發(fā)布的海冰密集度數(shù)據(jù)以及CS-2或S3的波形區(qū)分海冰、開闊水和冰間水道3種類型。常用后向散射系數(shù)（Radar Backscatter Coefficient,Sigma0）、波形前緣寬度（Leading Edge Width,LEW）和脈沖峰值（Pulse Peakiness,PP）3種回波波形特征進行類型識別[21]。Sigma0為雷達高度計接收到的地物的表面后向反射系數(shù)；LEW是最大峰值功率5%和95%點位之間的波形前緣寬度；PP是雷達波形最大峰值功率與同一采樣波形里的所有波形總功率的比率[4]，計算方法如下：

式中：NWF代表一個波形內(nèi)的距離門數(shù)；WF為雷達波形在第i個距離門處的功率。通常冰間水道的表面較平緩，雷達信號多為鏡面反射，PP值較高；對于開闊水域和海冰則以發(fā)生漫反射為主，PP值相對較低，LEW較寬（開闊水域的LEW相比海冰更大）。常用的分類參數(shù)值設(shè)置[22-23]見表1。

表1 CryoSat-2和Sentinel-3A的雷達波形分類參數(shù)設(shè)置

雷達高度計是通過計算信號發(fā)射與返回的時間差和光速來計算雷達到星下地物點之間的距離，因此常需要預(yù)先設(shè)置信號發(fā)射與返回時的時間跟蹤點。但在實際情況中，雷達接收信號的實際時間跟蹤點與預(yù)設(shè)時間跟蹤點之間存在偏差，需要根據(jù)實際情況進行波形重跟蹤校正[24]。本質(zhì)上，波形重跟蹤是將發(fā)生主散射的波形位置定為重跟蹤點，然后計算重跟蹤點和預(yù)設(shè)跟蹤點之間的偏移，以校正雷達高度計質(zhì)心到星下地物點之間的真實距離[25]。

TFMRA方法（Threshold First Maximum Retracker Algorithm）[26]是海冰干舷高度反演中常用的波形重跟蹤方法。該方法針對冰間水道和海冰等不同的地物類型，經(jīng)驗性地設(shè)置重跟蹤點閾值位置。表2總結(jié)了當前最為主要的一些重跟蹤閾值組合。為了更全面地比較不同重跟蹤閾值組合對干舷反演的影響，并得到最優(yōu)閾值組合，本文在表2的3種閾值組合的基礎(chǔ)上重新構(gòu)建了8種閾值組合，分別對應(yīng)為(40%，40%)，(50%，40%)，(50%，50%)，(60%，40%)，(60%，50%)，(70%，40%)，(70%，50%)，(70%，60%)（括號內(nèi)第一個閾值對應(yīng)于冰間水道，第二個閾值對應(yīng)于海冰）。以S3數(shù)據(jù)為例，這些閾值在雷達回波波形中的位置示意圖見圖2。同時，需要指出的是，冰間水道的閾值通常大于海冰的閾值，這是因為對于雷達高度計冰間水道的后向散射常高于海冰[27]。

表2 文獻中常用的重跟蹤閾值組合

圖2 Sentinel-3A雷達高度計海冰回波波形中第30～60距離門處歸一化回波功率閾值位置示意圖

在完成重跟蹤閾值設(shè)置后，利用波形重跟蹤校正后的海冰表面高程減去相應(yīng)的平均海表面高度與海表面高程異常值，得到海冰雷達干舷高度FR（Radar Freeboard）。由于電磁波在穿過積雪層時的速度與真空中不同，所以還需要利用積雪厚度數(shù)據(jù)對雷達海冰干舷高度FR進行距離校正，以得到真正的海冰干舷高度F。因此，當雷達信號完全穿透積雪層時，進行積雪校正的海冰干舷高度計算公式[4]如下：

式中：cs為積雪層中雷達信號的傳播速度；c為雷達信號在空氣中的傳播速度；hs(1-cs/c)≈0.22hs為積雪校正值。需要注意的是，公式（2）的使用條件為：電磁波能夠完全穿透雪層直達雪-冰交界面。若電磁波無法穿透雪層，應(yīng)用上式計算海冰干舷高度將會不可避免地引入誤差。

圖3 2017年3月、4月與2018年4月北極地區(qū)不同閾值組合下海冰干舷高度與OIB海冰干舷高度散點圖

表3 CryoSat-2衛(wèi)星數(shù)據(jù)不同閾值組合計算的海冰干舷高度與OIB海冰干舷高度對比

將2017年3月、4月與2018年4月CS-2與S3在8種閾值組合方案下得到的海冰干舷高度分別與同期的OIB海冰干舷實測數(shù)據(jù)進行對比，從而確定兩種Ku波段合成孔徑雷達高度計針對海冰干舷高度反演的最優(yōu)波形重跟蹤閾值組合。為保證精度，在本文的處理中，積雪厚度數(shù)據(jù)采用的機載OIB數(shù)據(jù)提供的積雪厚度測量值。圖3和表3給出了不同閾值組合得到的海冰干舷高度與機載OIB得到的海冰干舷高度對比結(jié)果。

從表3中可以看出，對于CS-2衛(wèi)星數(shù)據(jù)，最優(yōu)的閾值組合為70%（冰間水道）、60%（海冰）。在該組合下得到的海冰干舷高度與機載OIB數(shù)據(jù)的海冰干舷高度非常接近。觀測表3可知，其平均絕對值差值和均方根誤差均為最小值，分別為0.060 7 m和0.077 6 m。

對于S3衛(wèi)星數(shù)據(jù)，由表4可知，最優(yōu)的閾值組合為50%（冰間水道）、50%（海冰）。在該組合下的平均絕對值差值和均方根誤差分別為0.075 0 m和0.098 2 m。

同時，由表3、表4中的數(shù)據(jù)可以發(fā)現(xiàn)，CS-2衛(wèi)星數(shù)據(jù)在最優(yōu)閾值組合下得到的海冰干舷高度，無論是反演精度還是相關(guān)性，均高于S3衛(wèi)星數(shù)據(jù)。

表4 Sentinel-3A衛(wèi)星數(shù)據(jù)不同閾值組合計算的海冰干舷高度與OIB海冰干舷高度對比

進一步分析圖3，可以發(fā)現(xiàn)，無論是哪種閾值組合（即便采用本文得到的最優(yōu)閾值組合），均普遍存在雷達高度計探測的海冰干舷高度大于機載OIB探測的海冰干舷高度的情況。這說明對于Ku波段雷達，積雪層不能視為透明，需要考慮其無法穿透積雪層的情況。當積雪覆蓋在一年冰上時，海冰在析鹽過程中形成的薄鹽水層向上滲入到積雪層中，從而改變了積雪的介電特性，降低了雷達信號的積雪穿透能力[31]。對于覆蓋在多年冰上的積雪，不僅積雪表面的粗糙度更高[32]，而且經(jīng)過多年的消融凍結(jié)循環(huán)過程，會形成包含多個不同密度層的積雪層[33]。這些因素綜合在一起，一方面會增加積雪層的后向散射強度，另一方面還會使雷達信號的主散射面由雪-冰交界面上移至積雪層中。這將導(dǎo)致海冰干舷高度被高估，從而引起海冰厚度的高估[14]。

3 積雪層對海冰干舷高度反演的影響分析

為詳細評估積雪層對雷達信號的影響，計算了CS-2與S3的雷達穿透系數(shù)f[13]：

式中：d為積雪厚度hs與式（2）中的積雪校正值0.22hs之和；FR_sat為雷達高度計測的雷達海冰干舷高度；FR_OIB為OIB得到的雷達海冰干舷高度；OIB雷達海冰干舷高度由OIB的海冰干舷高度減去積雪校正值0.22hs得到。式（3）中，當雷達穿透系數(shù)f=0時，表示雷達信號沒有穿入雪層，雷達的主散射面位于空氣-雪交界面；f=1表示雷達信號完全穿透雪層，雷達信號的主散射面位于雪-冰交界面。0＜f＜1表示雷達信號未完全穿透雪層，雷達信號的主散射面位于積雪層內(nèi)部。

圖4和表5給出了2017年3月-4月和2018年4月北極區(qū)域CS-2與S3各自的雷達穿透系數(shù)分布與統(tǒng)計結(jié)果。在圖4中還利用高斯分布對雷達穿透系數(shù)進行了擬合，并將雷達穿透系數(shù)的平均值作為擬合曲線的中心值。

圖4 雷達穿透系數(shù)分布直方圖

表5 CryoSat-2和Sentinel-3A的雷達穿透系數(shù)平均值

需要注意的是，圖4存在雷達穿透系數(shù)f＜0或者f＞1的情況，對于該問題Armitage和Ridout[13]已給出過分析和評估。這主要是由于：（1）在計算海冰干舷時，為了有足夠的覆蓋范圍與機載OIB數(shù)據(jù)進行比較，衛(wèi)星數(shù)據(jù)的平均時間是遠長于機載OIB數(shù)據(jù)的；（2）衛(wèi)星高度計和機載OIB系統(tǒng)的足印大小存在較大差異。由于衛(wèi)星海冰干舷與機載OIB海冰干舷的兩種測量差異，使得衛(wèi)星高度計的海冰干舷分布范圍大于機載OIB的海冰干舷分布范圍，從而出現(xiàn)雷達穿透系數(shù)f不總在[0,1]范圍區(qū)間的情況。

從圖4和表5可知，Ku波段雷達高度計對一年冰上積雪的穿透能力普遍強于多年冰。例如：CS-2對于一年冰上積雪的雷達穿透系數(shù)f為0.950，對于多年冰為0.889；S3對于一年冰上積雪的雷達穿透系數(shù)f為0.873，對于多年冰為0.856；對于全部海冰類型CS-2的雷達穿透系數(shù)f為0.912，S3的雷達穿透系數(shù)f為0.867。總體而言，CS-2的雷達穿透系數(shù)f大于 S3。

所以普遍來說，對于Ku波段雷達高度計，其主散射面均位于雪-冰交界面上方的積雪層中，不能再利用公式（2）進行積雪校正，必須發(fā)展新的校正方法。

4 積雪校正新方法與海冰干舷高度反演

由雷達穿透系數(shù)f的定義可知，當0＜f＜1時，雷達高度計測得的海冰干舷高度會比實際的海冰干舷高度高(1-f)hs，所以需要補償這一距離差。同時，還需對雷達信號穿透的那部分積雪厚度進行距離校正0.22fhs。因此，改正的積雪校正模型可表述如下[13]：

式中：雷達穿透系數(shù)f可取為表5中的統(tǒng)計值，所以式（4）可根據(jù)雷達高度計和海冰類型的不同，分別針對CS-2與S3，得到本文新改進的積雪校正公式。

對于CS-2一年冰：

對于CS-2多年冰：

對于S3一年冰：

對于S3多年冰：

為了評估新改進的積雪校正方法的海冰干舷高度反演性能，本文以機載OIB實測海冰干舷高度數(shù)據(jù)為基礎(chǔ)，將其與改進的積雪校正法的海冰干舷反演結(jié)果、通用積雪校正方法（公式（2））的反演結(jié)果和ESA發(fā)布的海冰干舷高度產(chǎn)品進行了對比分析。ESA海冰干舷高度產(chǎn)品是歐空局利用CS-2和S3數(shù)據(jù)制作的海冰干舷高度產(chǎn)品（L2級），該數(shù)據(jù)可公開下載（ftp://science-pds.cryosat.esa.int）。

表6 通用積雪校正方法、提出的積雪校正方法、ESA海冰干舷高度產(chǎn)品與機載OIB數(shù)據(jù)的海冰干舷高度對比

表6和圖5給出的是3種海冰干舷反演方法的結(jié)果與機載OIB海冰干舷高度的差異統(tǒng)計表和分布圖。根據(jù)圖與表的統(tǒng)計結(jié)果可知，對于ESA發(fā)布的海冰干舷高度產(chǎn)品，CS-2發(fā)布的海冰干舷高度產(chǎn)品值明顯偏大，而S3發(fā)布的海冰干舷高度產(chǎn)品值明顯偏小，但總體上S3海冰干舷高度產(chǎn)品的誤差低于CS-2。對于通用積雪校正方法，CS-2和S3反演的海冰干舷高度均偏高。而本文改進的積雪校正方法顯著校正了通用積雪校正方法存在的高估問題，CS-2和S3海冰干舷高度的反演精度無論是平均絕對值差值還是均方根誤差均小于前兩種方法。另外需要指出的是，本文所改進的積雪校正方法對于一年冰和多年冰，其精度均有所提高，呈現(xiàn)了良好的穩(wěn)定性。

圖5 海冰干舷高度差異分布圖

5 結(jié)論與討論

工作在Ku波段的CryoSat-2和Sentinel-3A合成孔徑雷達高度計是當前最先進的高度計。本文針對雷達回波重跟蹤閾值的位置，以及積雪對雷達信號穿透能力的干擾影響，發(fā)展了改進的積雪校正算法，提高了海冰干舷高度反演精度。主要結(jié)論有：

（1）對于CS-2雷達高度計最優(yōu)的波形重跟蹤閾值組合為70%（冰間水道）、60%（海冰）；對于S3雷達高度計，其最優(yōu)的閾值組合為50%（冰間水道），50%（海冰）。

（2）Ku波段雷達信號未完全穿透一年冰與多年冰上的積雪層，并且CS-2的雷達穿透系數(shù)大于S3。CS-2對于一年冰上積雪層的雷達穿透系數(shù)為0.95，對于多年冰為0.889；S3對于一年冰上積雪層的雷達穿透系數(shù)為0.873，對于多年冰為0.856；對于全部海冰類型而言，CS-2的雷達穿透系數(shù)為0.912，S3的雷達穿透系數(shù)為0.867。

（3）相比于通用積雪校正法和ESA產(chǎn)品法，在本文的改進積雪校正法中得到的海冰干舷高度與OIB數(shù)據(jù)之間的誤差最小，精度最高。

本文的研究僅使用了北極春季的數(shù)據(jù)，為提高算法的普適性，未來將會利用北極其他季節(jié)的數(shù)據(jù)開展雷達穿透系數(shù)分析，從而進一步優(yōu)化算法。在下一步的研究工作中，準備利用本文改進的方法開展南極地區(qū)的海冰干舷高度反演與校正。

致謝：感謝歐洲空間局提供的CryoSat-2和Sentinel-3A雷達高度計數(shù)據(jù)及海冰干舷高度數(shù)據(jù)；感謝美國冰雪中心提供Operation IceBridge數(shù)據(jù)。