張春玲,李由之,張 剛,邱振戈

(1. 河南省地圖院,河南 鄭州 450008; 2. 南京大學地理與海洋科學學院,江蘇 南京 210023;3. 北京四維遠見信息技術有限公司,北京 100070; 4. 中國測繪科學研究院,北京 100830;5. 上海海洋大學海洋科學學院,上海 201306)

淺海水深數據是淺海環(huán)境治理、淺海資源開發(fā)利用、淺海航行、登島作戰(zhàn)等經濟軍事應用所需要的重要基礎地理空間信息,準確、高效、經濟地獲取淺海水深數據是海洋測繪目標[1]。與傳統(tǒng)聲學手段的水深測量相比,光學衛(wèi)星遙感技術具有低成本、可重復觀測、不受通航條件限制的優(yōu)勢,已成為淺海水深測量的重要手段。雙介質立體攝影測量和激光雷達測深是現階段主要的兩種光學衛(wèi)星遙感測深手段。其中,雙介質立體攝影測量屬于被動光學遙感測量技術,能夠快速、高效、大范圍地獲取遙感影像覆蓋區(qū)域的水深信息,但水深反演精度相對較低[2];而激光雷達測深技術屬于典型的主動光學遙感手段,雖然水深測量精度較高,但覆蓋范圍有限且重訪周期長。為了充分發(fā)揮兩種技術手段的優(yōu)勢,需要探討一種適用性好、反演精度高的主被動融合測深方法。

在衛(wèi)星遙感數據源方面,被動光學遙感數據來源豐富,可以從IKONOS、Pleiades、GeoEye-1、WorldView-3、GF-2、ZY-3等系列衛(wèi)星獲取中高分辨率的立體像對數據[3-6]。主動光學衛(wèi)星起步較晚,美國NASA在2003年發(fā)射了全球首顆對地觀測的激光測高衛(wèi)星ICESat[7]。后續(xù)衛(wèi)星ICESat-2也于2018年發(fā)射并繼續(xù)執(zhí)行激光測高任務[8]。而我國2019年發(fā)射的GF-7衛(wèi)星也搭載了全波形星載激光測高載荷[9]。在雙介質立體攝影測量方面,文獻[10]提出了一種基于高分辨率衛(wèi)星多光譜立體像對的淺水水深測量方法,淺海島礁水深反演的相對誤差小于20%;文獻[11]利用WorldView-2衛(wèi)星影像,解決了雙介質立體攝影測量折射改正算法適用性差的問題,在甘泉島的淺海測深精度達2.08 m。在星載激光測深方面,文獻[12]驗證了ICESat-2的激光載荷ATLAS具備淺海測深能力,并推導了水底高程折射改正方法,試驗中探測的最大水深達38 m,均方根誤差小于0.6 m;文獻[13]針對ICESat-2(ATL03數據)提出了一種自適應橢圓濾波測深方法,在我國南海甘泉島和珊瑚島開展測深試驗,結果表明,測深的均方根誤差分別為0.48和0.79 m。通過上述研究發(fā)現,雙介質立體攝影測量水深的反演結果難以達到激光雷達測深的精度水平,其中一個原因是對水面高程估算的不確定性。

本文以我國南海淺海島礁為例,利用ICESat-2激光測高數據和WorldView-3立體像對數據,構建主被動光學衛(wèi)星遙感數據融合的測深模型,開展高精度、無實測控制點的淺海水深測量試驗。

1 研究區(qū)域與數據

研究區(qū)域位于中國南海西沙群島東部,該區(qū)域四周沙堤環(huán)繞,水體透明度高,漫射衰減系數[14](Kd)約為0.11 m-1,太陽光透水深度大于20 m,屬于典型的一類水體,適合開展光學遙感水深反演試驗。

被動遙感影像數據采用的是WorldView-3立體像對數據(見表1)。該衛(wèi)星搭載推掃式線陣相機,影像成像模型穩(wěn)定,可采集0.31 m分辨率的全色影像和1.24 m分辨率的多光譜影像。目前廣泛應用的是4個波段(藍、綠、紅、近紅外)的多光譜立體像對,全色光譜范圍為400～850 nm,在分辨率、透水波長及影像角度上均滿足雙介質立體攝影測量的要求。研究區(qū)的立體像對影像為2020年9月30日拍攝,影像中有少量的云和耀斑,對水深反演結果影響較小。

表1 WorldView-3和ICESat-2衛(wèi)星參數

主動激光雷達數據采用ICESat-2衛(wèi)星搭載激光載荷ATLAS的二級產品ATL03,該數據產品為全球定位光子數據,記錄了接收光子的經緯度、高程、激光束的采樣時間等信息。ATLAS激光測高系統(tǒng)采用光子計數探測體制,具有微脈沖、高重頻、窄脈寬的特點,激光波長為532 nm,單脈沖能量為48～170 μJ,重復頻率為10 kHz,脈沖寬度為1.5 ns,具體的衛(wèi)星參數見表1,能夠按照沿軌方向連續(xù)獲取星下點的高程[8]。研究區(qū)的ATLAS軌跡如圖1所示,圓點為激光束由南向北經過研究區(qū)的一軌數據(GT1L),采集時間為2019年1月17日。由于ATLAS采用光子計數探測體制,在極大地提高探測靈敏度的同時,回波信號周圍也充滿大量的背景噪聲,如圖2所示。

圖1 研究區(qū)域影像與ATLAS軌跡(20190117GT1L)

圖2 ATL03數據光子點云沿軌剖面

2 研究方法

2.1 雙介質立體攝影測量模型

雙介質立體攝影測量是根據立體攝影測量原理,同時考慮光線在大氣和水體兩種介質中折線傳播的特點,利用水底點與對應同名像點的幾何光學關系確定水底點三維坐標[1]。光學衛(wèi)星在雙介質立體成像的幾何關系如圖3所示。(點A為攝站S1、S2垂直面上的理論交會點,點P為S1P2和S2P2反向延長線的交點,即水底點。p1、p2分別為左右圖像上的像點,其中,像點p1在直線S1P1上,S1為左影像相機鏡頭的中心位置,P1為成像光線與水面的交點,折線S1p1P為左影像上像點p1的成像光線。同理,右影像上同名像點p2的成像光線為折線S2p2P)。

圖3 光學衛(wèi)星在雙介質立體成像的幾何結構

假設水面近似水平面,為了得到折射改正后的水深h,按照單介質攝影測量原理先求解出水面到理論交匯點A的距離hA,則有

hA=Z0-ZA

(1)

式中,ZA為理論交匯點A的高程;Z0為估算的水面高程,即水面點P1和P2的高程。一般Z0是通過識別水陸岸線內插出對應像元的高程,取平均值作為水面高程。

由于水面點的連線P1P2平行于飛行方向,水深h的折射改正由同名光線S1P2和S2P2的入射角r1和r2確定,而折射角i1和i2可根據入射角r1和r2的幾何關系計算得到。設P1和P2之間的距離為k,根據同名光線的幾何關系可得[15]

k=(tanr1+tanr2)·hA=(tani1+tani2)·h

(2)

式中,n為水的折射率,則有

(3)

同時考慮左右兩個三角形結構,由點A垂線兩側的三角關系可知

(4)

(5)

(6)

(7)

折射角i1和i2可進一步推導為

(8)

(9)

將式(8)和式(9)代入式(3)中,雙介質立體攝影測量的折射改正水深h為

(10)

2.2 激光雷達測深模型

激光雷達測深的基本原理是通過主動發(fā)射一束532 nm波長的激光脈沖,依次穿過大氣、水面、水體及水底,激光探測器接收并記錄該脈沖在幾何路徑上光子的飛行時間,通過檢測水面和水底的回波信號往返的時間差,結合光速反演水深。由于ATL03中包含了大量的噪聲光子,有效識別水面和水底的信號光子是ICESat-2水深反演的核心問題。假設噪聲光子是隨機分布的,信號光子相較于噪聲光子在空間分布上更加密集。根據光子點云的離散分布特征,采用自適應橢圓濾波測深方法[13]分別提取水面光子(Zs)和水底光子(Zb)。

由于ATL03數據中地理坐標沒有考慮氣-水界面處折射率的變化,因此水底光子的高程需要進行折射改正。激光折射改正幾何結構關系如圖4所示,根據斯涅爾定律,折射角θ推導公式為

圖4 激光折射改正幾何結構

(11)

入水光線R和S之間的關系可表示為

(12)

(13)

而水底光子校正前與校正后的光子之間的距離P利用余弦定理可得[12]

(14)

其中,α和β的表達式分別為

(15)

(16)

因此,在垂軌方向和垂直高程上的偏移量ΔY和ΔZ可推導為

(17)

式中,na、nw分別為空氣和水的折射率;φ為激光入射角;θ為折射角,則夾角ψ=φ-θ;hA為直線傳播條件下水面與水底光子的高程差;S為未折射改正的水下斜距;R為改正后的水下斜距;P為水底光子折射改正前后的距離。

由于ICESat-2幾乎垂直對地觀測,折射角θ的最大值為0.38°,因此,偏移量ΔY的大小可以忽略,僅計算垂直方向的折射改正量ΔZ。則激光雷達的折射改正水深公式為

h=hA-ΔZ

(18)

2.3 主被動融合測深模型

雙介質立體攝影測量模型是建立在平靜水面的基本假設,但實際情況中的風速、涌浪等環(huán)境要素會引起水面高程的變化,改變立體雙介質中同名光線的幾何結構,從而導致測深誤差。而激光雷達測深模型中可以獲得沿軌方向上的瞬時高程,提取的水面光子具備高程的動態(tài)范圍。因此,主被動融合測深模型是將激光雷達反演的水面高程作為雙介質立體攝影測量水面高程糾正的起算位置。則水面激光點到雙介質理論交匯點的距離hA為

hA=Zs-ZA

(19)

此外,雙介質攝影測量中海面曲率變化也會引起同名光線S1P2和S2P2的入射角r1和r2計算誤差,但這部分的幾何模型誤差與水面高程無關。因此,在海面波浪較小的情況下,可以利用區(qū)域網平差的簡單多項式擬合水深改正誤差。試驗采用雙線性二次多項式的誤差改正模型,在各點深度中加入適當的改正數δh,公式為

δh=C0+C1X+C2Y+C3X2+C4XY

(20)

式中,C0、C1、C2、C3、C4為待定系數,各待定系數對于一個局部小區(qū)域而言是常數。由式(20)可知,只需5個以上改正后的水深點作為控制點,即可以列出5個以上方程式,求解5個待定系數。在求出各待定系數后,代入式(20)可計算出各水下點的改正數,從而求得其他點的水深。

3 試驗結果與分析

為了實現高精度、無實測控制點的淺海水深測量,利用WorldView-3立體像對影像范圍內的ICESat-2激光點云提取水面高程,輸入主被動融合測深模型,獲得研究區(qū)的淺海水深反演結果?？紤]研究區(qū)內的實測水深難以收集,本文將ICESat-2折射改正后的水底高程作為精度驗證的參考依據。

由于ATLAS軌跡由南向北直接穿過島礁,水面光子也被島礁分割成南北兩部分。島礁部分的陸地光子仍需人工標記,再根據前文方法提取出水面光子的高程,結果如圖5所示。需要說明的是,兩個區(qū)域提取的水面平均高程為3.72 m。其中,最大高程為5.61 m,最小高程為1.82 m,且均在島礁北部,這說明北部的風浪更大,導致水面高程相較于南部變化較大。

圖5 水面光子提取結果

主被動融合測深結果如圖6所示,地理坐標統(tǒng)一為WGS-84坐標系,高程采用大地高。從定性角度上看,水深結果層次分明,區(qū)域內的最大水深為22.51 m,與影像中目視透水深度基本一致。定量評估采用與激光水面高程同測線折射改正后的水底高程作為精度驗證的依據。圖7分別為島礁南北兩個區(qū)域激光水下點的折射改正結果,藍色點表示水面光子,紅色點表示折射改正前的水底光子,綠色點為改正后的水底光子。由于光子計數的離散特性,采用三次樣條曲線獲取連續(xù)的激光水底高程,擬合后的最大水深為32.13 m,能夠滿足精度驗證的需求。

圖6 主被動融合測深結果

圖7 激光水下點折射改正結果

將島礁南北兩個區(qū)域的主被動融合測深結果與雙介質立體攝影測量結果同步比較,統(tǒng)計了平均絕對誤差(MAE),標準差(SD),均方根誤差(RMSE)和R2,結果見表2。主被動融合測深模型在島礁南北的RMSE分別為0.69和0.28 m,而立體雙介質模型的RMSE分別為0.99和2.47 m。此外,主被動融合測深模型和立體雙介質模型在南部的絕對精度分別為0.48±0.50和0.78±0.61 m,在北部的絕對精度分別為0.23±0.16和2.25±1.02 m。相比之下,主被動融合測深模型的總體測深精度更高,計算的R2分別為0.94和0.98,特別是在島礁北部的精度有顯著提升。這說明該模型適用于水面高程變化較大的區(qū)域,精度上要優(yōu)于傳統(tǒng)的雙介質立體攝影測量,能夠有效彌補立體雙介質模型水面高程不確定性帶來的測深誤差。

表2 主被動融合測深精度 m

4 結 論

本文基于ICESat-2激光測高數據和WorldView-3立體像對數據,對我國南海島礁開展光學衛(wèi)星遙感數據融合的淺海水深測量,得到以下結論:

(1)本文模型在研究區(qū)的RMSE分別為0.68和0.28 m,表現出較好的精度水平,能在一定程度上解決雙介質立體攝影測量由于水面高程不確定帶來的測深精度問題。

(2)從研究區(qū)的ICESat-2水深反演結果來看,ATLAS反演的水深能夠達30 m以上,在一定程度上表現出較強的水體透射能力。

(3)由于ICESat-2激光測高數據和WorldView-3立體像對數據不一定在同一時間采集,ICESat-2反演的水面高程不一定表示WorldView-3立體像對影像中的瞬時海面。但現有的在軌衛(wèi)星難以獲得同一海域的時間同步衛(wèi)星遙感立體影像和多波束激光雷達測深數據,建議國家研制主被動融合淺海測繪衛(wèi)星,同時搭載高靈敏度、高動態(tài)范圍海洋立體測繪相機和多波束海洋測深激光雷達。

綜上所述,利用ICESat-2和WorldView-3主被動光學衛(wèi)星遙感數據融合的淺海水深測量,結合了主被動遙感手段的特點,在技術路線上合理可行,能夠實現高精度淺海遙感三維測繪,為滿足淺海測繪、大區(qū)域珊瑚礁遙感調查監(jiān)測、海岸帶和淺海遙感地質調查的需求提供技術手段。