基于貝塞爾曲線的波形重跟蹤改進(jìn)算法

吳 斌 李 興

(山東科技大學(xué) 測(cè)繪與空間信息學(xué)院, 山東 青島 266590)

0 引言

內(nèi)陸湖泊作為出海水以及冰川以外的重要水資源,在維護(hù)生態(tài)平衡、減輕洪澇災(zāi)害以及滿足日常生活需求等方面發(fā)揮著不可替代的作用[1]。湖泊范圍的擴(kuò)張與縮減、水位的上升與下降以及帶來(lái)生態(tài)環(huán)境的變化是全球以及局部氣候或者構(gòu)造共同的結(jié)果。湖泊的變化與人的生活息息相關(guān),它可以反映一個(gè)人地區(qū)的氣候變化、環(huán)境變化、濕度變化,降水量變化,因此有效監(jiān)測(cè)湖泊對(duì)局部地區(qū)的變化有重要的意義[2]。

中國(guó)有眾多的內(nèi)陸水資源,面對(duì)復(fù)雜的地理環(huán)境以及氣候變化,如何在雨季和旱季很好地監(jiān)測(cè)湖泊尤為重要。反映湖泊變化最直接的指標(biāo)就是水位的變化,監(jiān)測(cè)水位變化可以通過(guò)建設(shè)水文站,但水文站存在一定程度上的缺陷。一方面不是所有的河流都適合建立水文站,例如，青藏高原的許多湖泊,另一方面水文站是呈點(diǎn)狀分布的,對(duì)連續(xù)水位的變化監(jiān)測(cè)并沒(méi)有那么有效。因此采用衛(wèi)星測(cè)高,即雷達(dá)高度計(jì)可以有效解決以上問(wèn)題,國(guó)內(nèi)外眾多學(xué)者利用雷達(dá)高度計(jì)對(duì)內(nèi)陸湖泊進(jìn)行監(jiān)測(cè)[3-4]。雷達(dá)高度計(jì)測(cè)高過(guò)程中存在誤差,需要進(jìn)行距離改正,改正的過(guò)程稱為波形重跟蹤。大部分重跟蹤方法更多的應(yīng)用在海冰以及遠(yuǎn)海,但近年來(lái)越來(lái)越多學(xué)者將算法應(yīng)用于近海以及內(nèi)陸湖泊[5-6]。Hwang[7]提出一種改進(jìn)的閾值算法對(duì)臺(tái)灣的近海波形進(jìn)行重跟蹤,最終得到的海面高程精度好于先前的海面高程結(jié)果。Tseng[8]引入一種統(tǒng)計(jì)方法,通過(guò)引入深海海洋波形倆修復(fù)近海破損的波形進(jìn)而保存更多波形,最后利用閾值法等算法來(lái)反演出精度更高的海平面高程。高永剛[9]在研究?jī)?nèi)陸湖泊時(shí)提出一種改進(jìn)波形重構(gòu)算法,利用該算法對(duì)國(guó)內(nèi)湖泊研究過(guò)程中均取得較好的結(jié)果。以上算法大多是基于物理參數(shù)或者是傳統(tǒng)的閾值法和傳統(tǒng)重心偏移(the Offset centre of gravity,OCOG)等[10-11],在Lee[12]和Deng[13]有關(guān)于物理參數(shù)法重跟蹤的表述,物理參數(shù)法存在著解算復(fù)雜度高,傳統(tǒng)閾值法的結(jié)果精度相對(duì)較差,因此越來(lái)越多的學(xué)者提出一種波形擬合重跟蹤算法用于海冰以及海水的參數(shù)反演。本文的算法是在貝塞爾曲線的基礎(chǔ)上,提出的一種基于主波峰提取的貝塞爾曲線擬合算法,曲線擬合法波形重跟蹤原始的貝塞爾曲線對(duì)于波形復(fù)雜程度較低的海冰波形具有較好的擬合效果,在Shen[14]的文章中有所體現(xiàn),而湖泊存在比較復(fù)雜的波形,需要首先提取出主波峰,然后擬合出準(zhǔn)確的波形前緣部分可以更好地用于波形重跟蹤。

1 波形重跟蹤方法

1.1 水位反演原理

星載雷達(dá)高度計(jì)主動(dòng)發(fā)射脈沖經(jīng)地表反射后由接收機(jī)接受后,測(cè)量得到發(fā)射脈沖和接受脈沖連個(gè)時(shí)刻的時(shí)間差Δt,從而確定衛(wèi)星質(zhì)心到星下點(diǎn)的距離,最后得處相應(yīng)的高度,距離計(jì)算如式(1)所示。

(1)

式中,R是衛(wèi)星與地面點(diǎn)的距離;c為光速(299 792 458 m/s)。

回波波形的上升部分對(duì)應(yīng)為前緣,前緣的半功率點(diǎn)對(duì)應(yīng)的采樣時(shí)間對(duì)應(yīng)時(shí)間延遲,借此可以計(jì)算出水面高度信息。雷達(dá)高度計(jì)在經(jīng)過(guò)海面或者湖面以外的地物時(shí)波形會(huì)發(fā)生變形,從而導(dǎo)致原本的重跟蹤點(diǎn)與實(shí)際跟蹤點(diǎn)有差別,進(jìn)而會(huì)影響星地距離的計(jì)算。因此需要采用合適的方法對(duì)波形進(jìn)行重跟蹤,進(jìn)而獲得改正值以對(duì)高度計(jì)到地物的距離進(jìn)行改正,改正計(jì)算如式(2)、(3)所示。

式中,Crt、Cnt分別為改正后的重跟蹤點(diǎn)和原定的重跟蹤點(diǎn);t為脈沖寬度;Rc為改正后距離;R為原始測(cè)得的衛(wèi)星到星下點(diǎn)的距離。Cnt在低分辨率模式(low resolution mode,LRM)下一般取63,t一般取3.125 ns。

1.2 改進(jìn)的貝塞爾曲線波形重跟蹤法

本文采用的波形重跟蹤方法是改進(jìn)的貝塞爾曲線擬合法波形重跟蹤,改進(jìn)的貝塞爾曲線擬合法波形重跟蹤增加了一個(gè)主波峰提取的流程,主波峰提取[12]如式(4)所示。

(8)

式中,B是模擬波形的功率值;Pi代表參數(shù),需要擬合產(chǎn)生;t為坐標(biāo)的橫坐標(biāo)值。

經(jīng)過(guò)波形重跟蹤后可以得到距離的改正值,而在衛(wèi)星測(cè)高中還會(huì)受到對(duì)流層、潮汐以及電離層等的影響,因此下一步需要對(duì)距離進(jìn)行干濕對(duì)流層校正、潮汐改正、電離層改正以及海況偏差校正等,相關(guān)的改正參數(shù)在L1b數(shù)據(jù)中均可以得到,這里將所有的改正值得和記為ΔR。最終水位反演的計(jì)算如式(9)所示。

(9)

式中,H代表最終的水位;Halt為測(cè)高衛(wèi)星的高度;Hgeoid為大地水準(zhǔn)面的高度。

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)研究區(qū)及實(shí)驗(yàn)數(shù)據(jù)

本文的研究區(qū)域巢湖在安徽省中部,是位于長(zhǎng)江中下游的五大淡水湖之一,東西長(zhǎng)55 km、南北寬21 km,湖岸線周長(zhǎng)176 km,平均水深2.89 m,面積約780 km2。研究區(qū)域在梅雨季節(jié)降雨嚴(yán)重的月份水位會(huì)大漲,本文利用過(guò)境的CryoSat-2數(shù)據(jù)對(duì)湖水水位進(jìn)行監(jiān)測(cè)。

本文采用的CryoSat-2雷達(dá)高度計(jì)由歐洲空間局(European Space Agency,ESA)于2010年發(fā)射,該衛(wèi)星是Ku波段的干涉雷達(dá)高度計(jì),主要用于極地海冰監(jiān)測(cè),進(jìn)而用于分析全球氣候變化與極地海冰變化的聯(lián)系等等。衛(wèi)星總共有三種工作模式,分別是低分辨率指向星下點(diǎn)的高度計(jì)測(cè)量模式(LRM),合成孔徑雷達(dá)(synthetic aperture radar,SAR)測(cè)量模式以及干涉合成孔徑雷達(dá)模式(SAR interferometry,SARIn)。本文采用經(jīng)過(guò)研究區(qū)域的LRM模式下的L1b(level 1 B)數(shù)據(jù),該數(shù)據(jù)是D基線數(shù)據(jù)以新的格式存儲(chǔ)。實(shí)驗(yàn)數(shù)據(jù)采用的是8月9日足跡經(jīng)過(guò)巢湖的數(shù)據(jù),研究區(qū)域和足跡如圖1所示。

圖1 研究區(qū)域

2.2 數(shù)據(jù)預(yù)處理

對(duì)于CryoSat-2雷達(dá)高度計(jì)的LRM數(shù)據(jù)下載完成之后,依據(jù)巢湖的流域范圍對(duì)數(shù)據(jù)進(jìn)行篩選,同時(shí)借助L1b數(shù)據(jù)中提供的地物類別幫助篩選出湖水波形。將篩選得到的數(shù)據(jù)進(jìn)行進(jìn)一步剔除其中受到噪聲影響較為嚴(yán)重的波形,篩選方法是將返回波形功率值波峰出現(xiàn)在前20個(gè)跟蹤門(mén)或者后10個(gè)跟蹤門(mén)的剔除。

2.3 基于貝塞爾曲線的擬合波形

利用式(4)～式(5)提取出主波峰后,在向波形前緣和后緣分別擴(kuò)展兩個(gè)點(diǎn)后,為保證擬合的精度,將前緣向前延伸5到10個(gè)點(diǎn),后緣延伸的點(diǎn)為主波峰峰值的80%、60%、40%、20%點(diǎn)處。度延伸后的波形進(jìn)行分段,分段點(diǎn)位于波峰前緣部分峰值的20%處和80%處,波形后緣部分峰值的80%處以及20%處,對(duì)每一段進(jìn)行單獨(dú)的曲線擬合,擬合后的結(jié)果如圖2所示。從圖2可以看出,該曲線對(duì)波形前緣有很好的擬合效果,因?yàn)楸疚牡乃惴P(guān)鍵在于波形前緣的擬合,所以對(duì)于波形后緣的擬合效果并不好,但不影響水位的反演。

圖2 貝塞爾曲線擬合結(jié)果圖

2.4 實(shí)驗(yàn)結(jié)果分析

首先對(duì)本文算法得到的水位結(jié)果用3σ準(zhǔn)則剔除異常值,去除異常值的目的是將明顯存在測(cè)量誤差的數(shù)據(jù)剔除。將本文改進(jìn)波形重跟蹤的算法得到的湖水水位與CryoSat-2 LRM的2級(jí)數(shù)據(jù)中大三種重跟蹤算法以及實(shí)測(cè)水位數(shù)據(jù)進(jìn)行對(duì)比,LRM的2級(jí)數(shù)據(jù)中包含的三種算法分別是改進(jìn)的OCOG、改進(jìn)的用戶定制法(customer furnished item,CFI)以及倫敦大學(xué)學(xué)院法(University College London,UCL)[17]。圖3中展示的為本文算法以及LRM的2級(jí)數(shù)據(jù)中的結(jié)果,橫坐標(biāo)的1代表本文算法,2代表改進(jìn)的CFI,3代表UCL,4代表改進(jìn)的OCOG。目視結(jié)果來(lái)看,改進(jìn)的CFI算法的結(jié)果存在明顯的異常值,本文算法的水位值分布較集中,僅低于改進(jìn)的OCOG算法。在剔除異常值后發(fā)現(xiàn),本文算法實(shí)際的有效觀測(cè)值更多,多于改進(jìn)的CFI以及UCL,與改進(jìn)的OCOG算法相比有效觀測(cè)值相差無(wú)幾,去除異常值后的最大水位值與當(dāng)天水位實(shí)測(cè)值的差值在0.12 m左右,基本上滿足測(cè)量要求。

(a)去除異常值前

在比較本文算法與L2中三種算法的結(jié)果,需要評(píng)估出已知算法的精度結(jié)果,因此用已有的實(shí)測(cè)水位數(shù)據(jù)進(jìn)行評(píng)估。由于實(shí)測(cè)水位數(shù)據(jù)只有單個(gè)測(cè)站點(diǎn),因此將與本文軌道數(shù)據(jù)重合度最高的巢湖忠廟站水位數(shù)據(jù)作為雷達(dá)高度計(jì)各個(gè)測(cè)量點(diǎn)的實(shí)測(cè)數(shù)據(jù)。已知的三種算法結(jié)果與實(shí)測(cè)水位評(píng)估的結(jié)果精度如表1所示,從表1可以看出改進(jìn)的OCOG算法的差值絕對(duì)值得最大值、最小值以及均方根誤差都是最小,因此本文將采用此數(shù)據(jù)作為主要驗(yàn)證數(shù)據(jù)與本文算法進(jìn)行對(duì)比。

表1 L2數(shù)據(jù)中三種算法精度對(duì)比單位：m

表2展現(xiàn)的是本文算法以及其他兩類算法與改進(jìn)的OCOG算法的精度對(duì)比,從表中可以看出本文算法在絕對(duì)值差值最小值有0.081 m、絕對(duì)值差值最大值有0.870 m,均方根誤差為0.770,相關(guān)系數(shù)為0.892,均優(yōu)于其他兩個(gè)算法。同時(shí)從水位圖中可以看出本文算法與改進(jìn)的OCOG算法在水位值上重合水位數(shù)量更多。所以,利用本文算法反演出的水位值具有比CFI和UCL具有更高的可信度且具有更多有效觀測(cè)值。

表2 本文算法與L2中算法精度對(duì)比

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

本文提出波形重跟蹤算法在貝塞爾曲線的基礎(chǔ)上,加入了主波峰提取這一步驟,更好地解決了擬合復(fù)雜的多波峰內(nèi)陸湖波形這一問(wèn)題,貝塞爾曲線擬合需要足夠的點(diǎn),因此在主波峰提取后向前緣和后緣多延伸了3到5個(gè)點(diǎn)以保證擬合的準(zhǔn)確度,為保證準(zhǔn)確找到重跟蹤點(diǎn),經(jīng)過(guò)多次實(shí)驗(yàn)將重跟蹤點(diǎn)確定為主波峰最大值的50%處。實(shí)驗(yàn)結(jié)果表明，本文算法可以很好地解決貝塞爾曲線擬合復(fù)雜內(nèi)陸湖水波形的問(wèn)題,且保存了更多的有效水位值,水位精度與現(xiàn)有結(jié)果相比基本滿足實(shí)際測(cè)量要求。本文算法的不足之處在于無(wú)法對(duì)一些復(fù)雜受損的波形進(jìn)行修復(fù),從而反演出更多的有效觀測(cè)值,另外本文算法在高原地區(qū)的實(shí)用性還有待探索,還需要更多的研究進(jìn)行優(yōu)化。