CryoSat－2衛(wèi)星海冰區(qū)域波形識別及海冰干舷高確定

王立偉 金濤勇 張勝軍 李大煒

1 武漢大學(xué)測繪學(xué)院，武漢市珞瑜路129號，430079

CryoSat－2衛(wèi)星是繼ERS－1、ERS－2、Envisat及ICESat后首個針對極地進(jìn)行觀測的衛(wèi)星，在高緯區(qū)域具有30d的子周期［1］。衛(wèi)星搭載Ku波段雷達(dá)高度計(jì)SIRAL，有3種測量模式：LRM（低分辨率測量模式）、SAR（合成孔徑雷達(dá)測量模式）、InSAR（干涉合成孔徑雷達(dá)測量模式）［2］。考慮到數(shù)據(jù)源及數(shù)據(jù)處理方式的一致性，本文選擇基于SAR 模式的觀測數(shù)據(jù)。為顧及數(shù)據(jù)的現(xiàn)勢性，使計(jì)算出的海冰干舷高能更好地反映某一時刻某一區(qū)域的特征，計(jì)算過程中采用沿軌內(nèi)插方法，而沒有采用常用的區(qū)域格網(wǎng)化方法。將兩種方法的計(jì)算結(jié)果進(jìn)行比較發(fā)現(xiàn)，二者在多年冰區(qū)域存在略微差異，但整體趨勢較為一致。

1 CryoSat－2衛(wèi)星波形重跟蹤

在近岸區(qū)域或海冰區(qū)域，衛(wèi)星測高波形通常會受到較為嚴(yán)重的污染，在數(shù)據(jù)處理之前，首先需要進(jìn)行波形重跟蹤，計(jì)算出星載預(yù)設(shè)門與重跟蹤門之間的距離，進(jìn)而對測量得到的衛(wèi)星到星下點(diǎn)之間的距離進(jìn)行改正。常見的波形重跟蹤方法包括重心偏移法（OCOG 法）、閾值法、改進(jìn)閾值法、神經(jīng)網(wǎng)絡(luò)算法、Beta－5參數(shù)算法、Beta－9參數(shù)算法以及海洋算法等［3－9］。

歐空局所發(fā)布的SAR 數(shù)據(jù)分為3 個級別：L1B數(shù)據(jù)，即波形數(shù)據(jù)；L2數(shù)據(jù)，直接給出相對于WGS84橢球的表面高；L2I數(shù)據(jù)，在L2數(shù)據(jù)基礎(chǔ)上增加海冰區(qū)域的表面識別信息，以便計(jì)算海冰干舷高（海冰露出水面的高度）。由L1B 數(shù)據(jù)得到的L2 數(shù)據(jù)和L2I數(shù)據(jù)的處理方式統(tǒng)一采用OCOG 波形重跟蹤算法，而OCOG 算法最初是針對海洋波形提出的，并不適合海冰的返回波形，而閾值法更適合于海冰區(qū)域的波形重跟蹤［1，8－11］。

閾值法的思路是：首先計(jì)算或者直接利用OCOG 方法獲得波形的振幅，然后對振幅乘以一個閾值因子，獲得重跟蹤點(diǎn)的幅值，最后線性內(nèi)插獲得重跟蹤點(diǎn)的采樣位置（圖1）。計(jì)算公式如下：

圖1 閾值法重跟蹤示意圖Fig.1 Sketch map of threshold retracking method

式中，p（i）為第i個采樣點(diǎn)所對應(yīng)的回波功率；Amax為回波功率的最大幅值，即波形振幅；DC 為利用前c個采樣點(diǎn)所計(jì)算的噪聲水平；TL為重跟蹤點(diǎn)所對應(yīng)的幅值；t為幅值大于TL值點(diǎn)所對應(yīng)的最近采樣門；LEG 為通過線性插值所獲得的波形前緣中點(diǎn)，即波形重跟蹤點(diǎn)所對應(yīng)的采樣門［12］。需要說明的是，很多學(xué)者在計(jì)算波形的噪聲水平時通常采用波形的前幾個采樣值，一般默認(rèn)從第一個采樣值開始。有的測高衛(wèi)星由于第一個采樣值或前幾個采樣值存在較大偏差，實(shí)際計(jì)算時是利用后面的一段采樣值對噪聲水平進(jìn)行估計(jì)［8］。本文在實(shí)際計(jì)算時利用每個波形的前5個采樣值來估計(jì)噪聲水平，即a＝0，c＝5。

本文采用閾值法對CryoSat－2衛(wèi)星在海冰區(qū)域的SAR 波形進(jìn)行重跟蹤。由于雪密度和雪濕度等的影響，實(shí)際中雷達(dá)波形并不能完全穿透雪深。一些研究者將閾值系數(shù)直接選為50%，是假設(shè)雷達(dá)波形能夠完全穿透雪深，使得計(jì)算出的冰層厚度存在較大誤差。本文將閾值系數(shù)選為40%，即假設(shè)雷達(dá)波形只能穿透一定深度的雪深［1］。

2 海冰干舷高反演

2.1 海冰干舷高計(jì)算原理

計(jì)算海冰干舷高的關(guān)鍵是要獲得高精度海冰表面高和可靠的瞬時海面高，通常采用格網(wǎng)化的方法，計(jì)算格網(wǎng)內(nèi)的平均海冰高度和瞬時海面高度［8］。高精度海冰表面高可通過重跟蹤方法獲得，平均瞬時海面高一般取為格網(wǎng)內(nèi)浮冰之間的海面高度。因此，對于海冰干舷高反演來說，最重要的就是進(jìn)行浮冰之間的海水面（Lead）識別，準(zhǔn)確測量出浮冰之間的局部瞬時海面高，然后內(nèi)插計(jì)算海冰處的局部瞬時海面高。

格網(wǎng)法的缺點(diǎn)是，一個格網(wǎng)內(nèi)的數(shù)據(jù)可能會含有多條軌跡，幾條軌跡間的時間跨度可能很長，通過取中值或平均值所得結(jié)果并不一定代表某一時刻的瞬時海面高，從而影響海冰干舷高的計(jì)算精度。

本文采用單條軌跡的沿軌內(nèi)插方法計(jì)算海冰位置處的瞬時海面高。先將DTU10平均海平面高內(nèi)插到衛(wèi)星的每一個星下點(diǎn)，然后通過波形特征參數(shù)識別出海冰的位置。以每一個海冰位置為節(jié)點(diǎn)，向軌跡的兩端同時搜索，識別出與該海冰節(jié)點(diǎn)距離最近的Lead。如果兩端均搜索到Lead，則同時計(jì)算兩個Lead處海平面異常，并采用距離加權(quán)方法將其內(nèi)插到海冰位置處，再加上海冰處的平均海平面高得到海冰處瞬時海平面高；如果僅一端搜索到Lead，則直接將該Lead處的海平面異常視為海冰位置處的海平面異常。事實(shí)上，僅一端能搜索到Lead的情況非常少見，可以忽略。

2.2 波形識別方法

為了進(jìn)行Lead和海冰的識別，將OSI SAF（海洋與海冰監(jiān)測）高緯處理中心發(fā)布的海冰濃度格網(wǎng)數(shù)據(jù)內(nèi)插到CryoSat－2數(shù)據(jù)中，刪除那些海冰濃度小于70%的數(shù)據(jù)（主要是海洋數(shù)據(jù)），然后再利用脈沖峰值、波束標(biāo)準(zhǔn)差、波峰陡峭度這3個波形特征參數(shù)對Lead和海冰進(jìn)行區(qū)分。各參數(shù)的門限值設(shè)置情況如表1。

表1 波形識別參數(shù)設(shè)置Tab.1 Set the waveform parameter identification

脈沖峰值是波形回波功率的函數(shù)［1］：

pp值越大，說明波形越尖銳。波束標(biāo)準(zhǔn)差用來描述波形后向反射角的變化［2］，而波峰陡峭度用來描述波形頂端的陡峭程度［13］。

2.3 實(shí)驗(yàn)結(jié)果及比較

德國AWI利用40%閾值法分析2011－01～2013－03共27個月的數(shù)據(jù)，利用格網(wǎng)化方法計(jì)算得到海冰干舷高。選取其最后3個月的數(shù)據(jù)與本文結(jié)果進(jìn)行比較，如圖2。可以看出，二者得到的北冰洋海冰干舷高整體分布趨勢較為一致，AWI在多年冰區(qū)域的海冰干舷高要更高一些。圖3給出2013－01～03 二者的差值直方圖。可以看出，這3個月二者的差值直方圖非常接近于正態(tài)分布。最后利用OSI SAF發(fā)布的每天多年冰和一年冰格網(wǎng)數(shù)據(jù)，分別計(jì)算AWI和本文數(shù)據(jù)的多年冰區(qū)域和一年冰區(qū)域的平均海冰干舷高（圖4）。

從圖4看出，本文計(jì)算的2011～2013年北冰洋多年冰和一年冰區(qū)域的平均海冰干舷高結(jié)果與AWI較為一致，海冰干舷高呈現(xiàn)明顯的周期性變化。由于夏季海冰融化，使得波形識別受到影響，利用返回波形特征識別Lead和海冰的方法不再可用，所以沒有對6～8月的海冰進(jìn)行處理。AWI利用波形識別方法得到的6～8月北冰洋海冰干舷高很不理想。對于夏季海冰的識別，還需采用其他更為有效的方法。

圖5為2011～2013年北冰洋各年同期數(shù)據(jù)對比圖，無論是本文結(jié)果還是AWI結(jié)果，各年同期海冰干舷高結(jié)果均較為接近。但同時也要注意到，在2～5月，本文得到的2013年數(shù)據(jù)無論是在多年冰區(qū)域還是在一年冰區(qū)域都較同期其他數(shù)據(jù)偏低，而在9～12月略微偏高，尤其是2013－09的數(shù)據(jù) 較2011－09和2012－09呈現(xiàn)明顯的增長趨勢。由于AWI數(shù)據(jù)只公布到2013－03，所以為了驗(yàn)證本文結(jié)果，還需要采用實(shí)地測量數(shù)據(jù)或其他衛(wèi)星測量數(shù)據(jù)進(jìn)行更好的解釋和說明。

圖2 本文海冰干舷高計(jì)算結(jié)果與AWI計(jì)算結(jié)果的比較Fig.2 Comparison between the results of sea ice freeboard calculation in this paper and the results of AWI

圖3 2013－01～03AWI與本文海冰干舷高差值直方圖Fig.3 Histogram of the sea ice freeboard height difference（2013－01～03）

3 結(jié) 語

本文采用沿軌內(nèi)插方法計(jì)算北冰洋區(qū)域的平均海冰干舷高，與常用的格網(wǎng)化方法計(jì)算結(jié)果整體趨勢較為一致，說明沿軌內(nèi)插方法完全可以用來進(jìn)行海冰干舷高的反演。在計(jì)算中發(fā)現(xiàn)，2013－09以后北冰洋區(qū)域無論是多年冰區(qū)域還是一年冰區(qū)域，平均海冰干舷高均呈現(xiàn)明顯上升趨勢，這還需要通過其他數(shù)據(jù)進(jìn)行驗(yàn)證。為了對本文結(jié)果和AWI結(jié)果進(jìn)行更好的比較，下一步將采用實(shí)測數(shù)據(jù)和機(jī)載數(shù)據(jù)，對二者的計(jì)算結(jié)果進(jìn)行更詳細(xì)的比較，進(jìn)而評價(jià)格網(wǎng)化方法和沿軌內(nèi)插方法在計(jì)算海冰干舷高中的優(yōu)劣，最終將兩種方法合理組合，更好地反演海冰干舷高。

圖4 2011～2013年北冰洋平均海冰干舷高變化趨勢Fig.4 Variation trend of the average sea ice freeboard height from 2011to 2013year

圖5 2011～2013年北冰洋平均海冰干舷高同期數(shù)據(jù)比較Fig.5 Arctic sea ice freeboard period average data comparison from 2011to 2013

［1］Wingham D J，F(xiàn)rancis C R，Baker S，et al.CryoSat：A Mission to Determine the Fluctuations in Earth’s Land and Marine Ice Fields［J］.Adv Space Res，2006，37（4）：841－871

［2］Ricker R，Hendricks S，Helm V，et al.Sensitivity of CryoSat－2Arctic Sea－Ice Volume Trends on Radar－Waveform，Interpretation［J］.The Cryosphere Discuss，2014，8（2）：1 831－1 871

［3］高永剛，郭金運(yùn)，黃金維，等.湖泊范圍內(nèi)TOPEX／Poseidon衛(wèi)星測高波形分析［J］.大地測量與地球動力學(xué)，2008，28（1）：45－49（Gao Yonggang，Guo Jinyun，Huang Jinwei，et al.Waveform Analysis of Topex／Poseidon Satellite Altimeter over Lake［J］.Journal of Geodesy and Geodynamics，2008，28（1）：45－49）

［4］褚永海，李建成，張燕，等.ENVISAT 測高數(shù)據(jù)波形重跟蹤分析研究［J］.大地測量與地球動力學(xué)，2005，25（1）：76－80（Chu Yonghai，Li Jiancheng，Zhang Yan，et al.Analysis and Investigation of Waveform Retracking Data of Envisat［J］.Journal of Geodesy and Geodynamics，2005，25（1）：76－80）

［5］褚永海.衛(wèi)星測高波形處理理論研究及應(yīng)用［D］.武漢：武漢大學(xué)，2007（Chu Yonghai.The Research of Satellite Altimetry Waveform Theory and Application［D］.Wuhan：Wuhan University，2007）

［6］常曉濤，李建成，郭金運(yùn)，等.一種多前緣多閾值的波形重構(gòu)算法［J］.地球物理學(xué)報(bào)，2006，49（6）：1 629－1 634（Chang Xiaotao，Li Jiancheng，Guo Jinyun，et al.A Multileading Edge and Multi－threshold Waveform Retrackef［J］.Chinese J Geophys，2006，49（6）：1 629－1 634）

［7］European Space Research Institute（ESRI N）.European Space Agency and Mullard Space Science Laboratory－University College London［EB／OL］.http：／／earth.esa.int／web／guest／missions／cryosat／ipf－baseline，2013

［8］Rose S K.Measurements of Sea Ice by Satellite and Airborne Altimetry［D］.Denmark：DTU Space，National Space Institute，2013

［9］Davis C H.A Robust Threshold Retracking Algorithm for Measuring Ice－Sheet Surface Elevation Change from Satellite Radar Altimeters［J］.IEEE Trans Geosci Remote Sens，1997，35（4）：974－979

［10］Helm V，Humbert A，Miller H.Elevation Change of Greenland and Antarctic Derived from CryoSat－2［J］.The Cryosphere Discuss，2014，（8）：1 673－1 721

［11］Stenseng L，Andersen O B.Preliminary Gravity Recovery from CryoSat－2Data in the Baffin Bay［J］.Adv Space Res，2012，50（8）：1 158－1 163

［12］阮海林.衛(wèi)星測高波形重跟蹤算法在臺灣周邊海域的應(yīng)用研究［D］.廈門：國家海洋局第三海洋研究所，2010（Ruan Hailin.Research and Application of Algorithm in the Waters Surrounding Taiwan Satellite Altimetry Waveform Retracking［D］.Xiamen：The Third Oceanography Research Institute，2010）

［13］Brenner A C，Zwally H J，Bentley C R，et al.Derivation of Range and Range Distributions from Laser Pulse Waveform Analysis for Surface Elevations，Roughness，Slope，and Vegetation Heights［EB／OL］.http／／glas.wff.nasa.govl，2000