戚甲偉,任照宇,趙金秀,朱金山*,3
(1.地理信息工程國家重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710054;2.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院,山東 青島 266590;3.自然資源部 海洋測繪技術(shù)重點(diǎn)實(shí)驗(yàn)室,山東 青島 266590)
水深信息對近海和海岸帶規(guī)劃與管理、航海和海洋環(huán)境科學(xué)研究具有重要意義。船舶航道的監(jiān)測以及水下沙洲、巖石、淺灘、暗礁等海洋特征的監(jiān)測都依賴于精確的水深測量[1-2]。一些近岸活動(如娛樂、漁業(yè)和水產(chǎn)養(yǎng)殖等)和近海工程(如電纜和管道鋪設(shè)、疏浚、石油鉆探等),也需要了解水深的分布狀況。研究海洋底質(zhì)分類時,可利用水深信息進(jìn)行水體校正,去除水體對底質(zhì)的影響從而提高底質(zhì)分類的準(zhǔn)確性。
根據(jù)光在水體中的輻射傳輸理論,光在水體中呈指數(shù)衰減,LYZENGA[3]給出了淺水區(qū)水深與海表輻亮度的關(guān)系。海表面的輻亮度除了與水深有關(guān)外,還與水底的底質(zhì)類型以及水體本身的光學(xué)屬性有關(guān)。如果底質(zhì)類型單一,水深值與單波段的輻亮度存在較好的相關(guān)關(guān)系,可根據(jù)單波段來反演水深;當(dāng)?shù)踪|(zhì)類型不單一時,需要利用多個波段來計算水深[3-4]。1985年LYZENGA[4]提出了一個雙波段線性回歸模型(Lyzenga 1985模型),利用多光譜數(shù)據(jù)與激光雷達(dá)數(shù)據(jù)進(jìn)行了水深反演,其研究證明利用多光譜遙感數(shù)據(jù)進(jìn)行水深反演是可行的。論文還討論了水底底質(zhì)的變化和水體光學(xué)屬性的變化對水深反演精度的影響。
針對多種水底底質(zhì),STUMPF et al[5]提出了對數(shù)變換模型(Stumpf 2003模型),此模型建立了2個波段對數(shù)變換后的比值與水深的關(guān)系,需要根據(jù)實(shí)測水深數(shù)據(jù)來確定模型中的3個未知參數(shù)。該研究利用高分辨率的IKONOS多光譜數(shù)據(jù)和激光雷達(dá)數(shù)據(jù)進(jìn)行了水深反演與驗(yàn)證。其研究表明,受底質(zhì)影響,如果遙感影像中某些地方淺水區(qū)的海表反射率比深水區(qū)的海表反射率還要小,則不能采用線性模型反演。
本文選擇Lyzenga 1985和Stumpf 2003兩種模型,在南海附近水域進(jìn)行水深反演比較,分析各波段與水深的相關(guān)性,使用Levenberg-Marquardt(L-M)[6-8]最小二乘回歸方法求解模型參數(shù),統(tǒng)計回歸點(diǎn)與真實(shí)水深點(diǎn)之間的偏移量,比較不同底質(zhì)情況下兩種模型的反演精度。
南海位于中國大陸的南方,是太平洋西部海域,為中國三大邊緣海之一。研究區(qū)域位于南海西沙群島的北島海域(圖1),水域面積約為4 km2,水深小于20 m,變化緩慢,層次性較好。研究區(qū)淺海水域海水清澈,適合利用光學(xué)方法進(jìn)行淺海水深反演研究。
海島被灌木和低地包圍,四周由寬約80 m左右的白色沙堤環(huán)繞,沙堤由珊瑚、貝殼、沙屑組成。其潮坪區(qū)在高潮時被海水淹沒,低潮時露出海面。
WorldView-2是DigitalGlobe公司的第3顆運(yùn)行衛(wèi)星,運(yùn)行在770 km高的太陽同步軌道上,可捕獲0.5 m分辨率的全色圖像和1.8 m分辨率的多光譜圖像,覆蓋的光譜范圍從400 nm到1 040 nm。本研究使用4個多光譜波段(紅、綠、藍(lán)、近紅外)的數(shù)據(jù),數(shù)據(jù)采集時間為2010年2月。整幅圖像云層覆蓋率約1%,研究區(qū)域內(nèi)無云層影響。
水深數(shù)據(jù)采用CGCS2000國家大地坐標(biāo)系、Gauss-Kruger 3°帶投影,基準(zhǔn)面為當(dāng)?shù)貚u礁潮高基準(zhǔn)面,由“海南省三沙市部分島礁周邊礁盤海域水下地形測量項(xiàng)目”獲得。測量時使用HY1601測深儀和廣州南方靈銳S86型GPS接收機(jī),將測深儀與接收機(jī)進(jìn)行連接,測深儀每秒(水平點(diǎn)位間距約3 m)采集1個點(diǎn),兩者實(shí)時同步獲得點(diǎn)位的坐標(biāo)與深度,同步記錄該點(diǎn)位潮汐高度信息。
圖1 研究區(qū)域地理位置示意圖Fig.1 The location of research area
研究流程如圖2所示,對WorldView-2衛(wèi)星影像進(jìn)行必要的幾何校正、輻射校正、太陽耀斑校正以改正遙感影像幾何誤差,削弱太陽輻射和水面局部強(qiáng)反射的影響,使圖像信息真實(shí)可靠[9]。
采用6S (Second Simulation of the Satellite Signal in the Solar Spectrum)數(shù)值模型進(jìn)行大氣校正[10],將大氣頂輻亮度轉(zhuǎn)換為海表輻亮度。6S模型由Fortran語言編寫,根據(jù)WorldView-2的波段信息給模型添加WorldView-2的光譜響應(yīng)。模型輸入的參數(shù)包括:(1)幾何參數(shù),如太陽天頂角、太陽方位角、衛(wèi)星天頂角、衛(wèi)星方位角等;(2)大氣模型,本文選用美國標(biāo)準(zhǔn)大氣;(3)氣溶膠類型,本文選用海洋性氣溶膠;(4)地球表面反射模型,本文選用海洋的反射率模型。由于研究區(qū)北島長度僅約1.6 km,空間跨度不大,校正時對整個研究區(qū)采用了相同的大氣校正參數(shù)。
圖2 研究流程圖Fig.2 Flow chart of research
對于寬視場角影像來說,無法避免波浪導(dǎo)致的水面局部強(qiáng)反射(Sun Glint),特別是對空間分辨率小于10 m的遙感影像[11-12]。GOODMAN et al[13]研究表明,遙感影像中30%以上的水深反演誤差都是由這種強(qiáng)反射造成。目前針對高分辨率衛(wèi)星遙感數(shù)據(jù)的太陽耀斑校正方法主要包括HOCHBERG et al[14]提出的方法以及在此基礎(chǔ)上改進(jìn)的HEDLEY et al[11]的方法和LYZENGA et al[15]的方法。本文采用HEDLEY et al[11]的方法對WorldView-2衛(wèi)星數(shù)據(jù)進(jìn)行太陽耀斑校正。具體算法公式如下:
Lwb(λ)=L(λ)-b(λ)[L(NIR)-Lmin(NIR)]
(1)
其中:Lwb(λ)為采樣波段耀斑校正后的像元值;L(λ)為此波段校正前的像元值;Lmin(NIR)為近紅外波段的像元最低值;L(NIR)為近紅外波段的像元值;b(λ)是以近紅外波段為橫軸,采樣波段為縱軸回歸出的直線的斜率。
圖3為太陽耀斑校正前后的對比圖(圖像進(jìn)行了增強(qiáng)處理,使噪聲便于觀察),校正后圖像中的干擾圖斑已經(jīng)被改正。
圖3 太陽耀斑校正前后圖像(圖像增強(qiáng)處理后)Fig.3 Before and after sun glint correction (after image enhancement)
本研究采集的水深數(shù)據(jù)是以島礁潮高基準(zhǔn)面為基準(zhǔn)的穩(wěn)態(tài)水深,由于水深反演模型需要的是影像過境時對應(yīng)的瞬時海水深度,故需進(jìn)行潮汐改正得到影像的實(shí)際水深值[16]。對照2010年2月10日(衛(wèi)星數(shù)據(jù)獲取時間)的潮汐表(永興潮高基準(zhǔn)面),獲取影像成像時刻的潮高為79 cm。由于潮汐表與實(shí)測水深數(shù)據(jù)的潮高基準(zhǔn)面一致,因此,遙感影像成像時刻的實(shí)際水深值即為實(shí)測水深值與成像時刻的潮高之和。
為實(shí)現(xiàn)遙感影像與水深信息的配準(zhǔn),首先對實(shí)測數(shù)據(jù)進(jìn)行投影轉(zhuǎn)換,將遙感數(shù)據(jù)和實(shí)測數(shù)據(jù)統(tǒng)一到UTM、WGS-84投影坐標(biāo)系中。再選擇3個清晰、易獲取的特征點(diǎn)作為控制點(diǎn),進(jìn)行精確配準(zhǔn)。由于衛(wèi)星影像的1個像元反映的是1.84 m×1.84 m范圍的水深平均值,以其中1個實(shí)測點(diǎn)作為整個像元的水深值誤差會很大,因此將其重采樣為1.84 m空間分辨率,對實(shí)測序列進(jìn)行統(tǒng)計分析,與遙感數(shù)據(jù)匹配[17]。
Lyzenga 1985模型[4]可表示為:
(2)
其中:a0和ai是常系數(shù),N是參加反演的波段數(shù)量,R(λi)是波段的海表反射率,R∞(λi)是波段i對應(yīng)的深水區(qū)反射率均值。可采用直方圖統(tǒng)計的方法確定深水區(qū)像元,其原理為在淺水區(qū)由于存在水底反射,其像元值大于深水區(qū)像元,隨深度加大,像元值會隨之減小。因此直方圖中像元值最小的部分,且遠(yuǎn)離島嶼的像元可確定為深水區(qū)像元。Lyzenga 1985模型中的自然對數(shù)變換使得水深與光譜波段輻射之間產(chǎn)生了一種線性關(guān)系,它是目前應(yīng)用較廣泛的多光譜遙感影像水深估算模型。
Lyzenga 1985模型中需要減去深水區(qū)輻亮度值,這有可能會造成進(jìn)行對數(shù)變換的值為負(fù)數(shù),進(jìn)而使結(jié)果出現(xiàn)無理數(shù),使得這些區(qū)域的反演為無效值。為避免這種情況,STUMPF et al[5]提出了一種基于對數(shù)變換波段比的非線性水深反演模型:
(3)
其中:m1,m0和n是常系數(shù);R(λ2)和R(λ1)是2個波段的海表反射率。
Lyzenga 1985和Stumpf 2003兩種水深反演模型都基于這樣的假設(shè):(1)水體的光學(xué)屬性如衰減系數(shù),在研究區(qū)內(nèi)無空間變化,也就是說要滿足在相當(dāng)大的區(qū)域內(nèi)水質(zhì)均勻。若水質(zhì)確實(shí)表現(xiàn)出顯著的空間變化,那么模型的深度估計將會有很大的誤差。(2)對于給定圖像中的不同的底質(zhì)類型,2個光譜波段的海底反射率之比是相同的。這一假設(shè)存在一定局限性,會對一些區(qū)域的反演精度造成影響。綜上所述,我們需要從多光譜圖像中提取水體清澈度高、水體水質(zhì)和底質(zhì)類型均勻的采樣區(qū)域,以滿足反演模型的假設(shè)前提。
L-M算法[7-8]通過最小化目標(biāo)函數(shù)得到問題的最優(yōu)解,可用于最小二乘問題的參數(shù)求解,也可用于非線性問題的求解。對于本研究中采用的2個模型,公式(2)中有3個未知數(shù),可轉(zhuǎn)化為求解線性方程的最小二乘問題;公式(3)中存在3個未知數(shù),可轉(zhuǎn)化為非線性方程的參數(shù)求解問題。因此,這2個模型均可采用L-M算法獲得模型參數(shù)。
L-M算法中,給定參數(shù)的初始值,通過迭代更新參數(shù),直到新的參數(shù)能夠滿足收斂條件(目標(biāo)函數(shù)最小)。其迭代公式如下:
Pi+1=Pi+[H+λdiag(H)]-1f(Pi)
(4)
其中:Pi表示當(dāng)前參數(shù)向量,Pi+1表示下一步迭代中采用的新的參數(shù)向量,H是目標(biāo)函數(shù)的黑森矩陣,λ是衰減因子,f(Pi)表示目標(biāo)函數(shù)值的梯度。
針對2個模型定義目標(biāo)函數(shù)為:
(5)
其中:N表示建模時采用的點(diǎn)的個數(shù);zi表示第i個實(shí)測水深值;zmi表示第i個由模型得到的估計水深值。當(dāng)目標(biāo)函數(shù)最小時,表示模型估計值與實(shí)測值最接近,此時的模型參數(shù)為最優(yōu)模型參數(shù)。
2.HBO組患者治療前后細(xì)胞免疫功能的變化:經(jīng)1個療程治療后,HBO組患者的細(xì)胞免疫指標(biāo)較治療前均有顯著改善,差異均有統(tǒng)計學(xué)意義(P<0.05)。見表2。
公式(5)中,根據(jù)水深反演模型的不同,zmi可由公式(2)或公式(3)得到。
L-M流程圖如圖4所示,首先給出模型的初始值,由初始值和反演模型[公式(2)或公式(3)]得到模型估計的水深值;根據(jù)模型值與實(shí)測值計算目標(biāo)函數(shù)值[公式(5)];然后根據(jù)目標(biāo)函數(shù)值判定是否達(dá)到收斂條件,如果未收斂則根據(jù)公式(4)更新模型參數(shù),并將新的參數(shù)輸入到反演模型重新估計水深值。經(jīng)多次迭代直到達(dá)到收斂條件,得到最優(yōu)的模型參數(shù)。
圖4 L-M算法流程圖Fig.4 Flow chart for the L-M algorithm
研究區(qū)域水深較淺,水體透明度較好,由于WorldView-2有較高的空間分辨率,根據(jù)真彩色影像(圖5),研究區(qū)中存在兩種明顯不同的底質(zhì)。圖中A點(diǎn)標(biāo)注的位置為底質(zhì)較暗的區(qū)域,經(jīng)現(xiàn)場拍攝的視頻確認(rèn)為珊瑚底質(zhì)(底質(zhì)區(qū)域A);B點(diǎn)所示位置底質(zhì)較亮,經(jīng)現(xiàn)場拍攝的視頻確認(rèn)為砂質(zhì)底質(zhì)(底質(zhì)區(qū)域B)。
圖5 研究區(qū)域劃分Fig.5 Study area division
由于研究區(qū)存在兩種明顯不同的底質(zhì),根據(jù)底質(zhì)不同建立了2組數(shù)據(jù)集。數(shù)據(jù)集每個樣本包含實(shí)測水深值和WorldView-2四個波段的像元值。其中區(qū)域A共有410個樣本,300個樣本用于建模,110個樣本用于模型的驗(yàn)證;區(qū)域B共125個樣本, 100個樣本用于建模,25個樣本用于模型的驗(yàn)證。
不同波長的光穿透水時,衰減程度不同,波長越長衰減越快。在清澈的水中,藍(lán)波段(440~540 nm)的最大穿透深度約為30 m,綠波段(500~600 nm)的約為15 m,紅波段(600~700 nm)的約為5 m,近紅外波段(700~800 nm)的約為0.5 m[2]。透光深度受水體渾濁度的影響,懸浮沉積物、葉綠素和溶解的有機(jī)化合物增加了渾濁度,從而降低了光的穿透深度。在建立反演模型之前,將大氣校正與太陽耀斑校正后所得到的每個波段的海水反射率與水深值進(jìn)行相關(guān)性分析(表1),確保建立模型的波段是合理可靠的。從表中可以看出,對于單波段來說,綠波段與水深值相關(guān)性系數(shù)最高,而近紅外波段由于在水中迅速衰減,其與水深數(shù)據(jù)相關(guān)性最差。選擇相關(guān)性較高的藍(lán)波段和綠波段作為模型的反演波段。
表1 各波段反射率與水深值相關(guān)系數(shù)Tab.1 Correlation coefficient between reflectivity of each band and water depth value
對于區(qū)域A,為了保證評估的可靠性,我們從已經(jīng)建立的數(shù)據(jù)集中選擇了300個樣本用于建模,110個樣本用于模型檢驗(yàn),檢驗(yàn)點(diǎn)專門用于檢驗(yàn)?zāi)P偷恼`差,不參加建模過程。根據(jù)公式(1)~(4),對兩種反演模型回歸所進(jìn)行的迭代過程如表2和表3所示。由表可知,迭代次數(shù)與選取的回歸初始參數(shù)有關(guān),隨著迭代次數(shù)的增加,殘差平方和逐漸減小并趨于穩(wěn)定,模型中的參數(shù)逐漸趨向最優(yōu)值。
表2 Stumpf 2003模型的回歸迭代結(jié)果(區(qū)域A)Tab.2 Regression iteration results of Stumpf 2003 model(area A)
表3 Lyzenga 1985模型的回歸迭代結(jié)果(區(qū)域A)Tab.3 Regression iteration results of Lyzenga 1985 model(area A)
為了比較兩種模型的反演精度和優(yōu)劣,將Stumpf 2003模型與Lyzenga 1985模型的結(jié)果與真實(shí)水深進(jìn)行對比和分析,統(tǒng)計檢驗(yàn)點(diǎn)水深反演值的決定系數(shù)(R2)和均方根誤差(RMSE)。結(jié)果如圖6和表4所示。
圖6 反演水深與實(shí)測水深對比(建模數(shù)據(jù)集,區(qū)域A)Fig.6 Comparison of the retrieved depth and the in-situ depth(modeling dataset, area A)
表4 兩種反演模型的誤差比較(建模數(shù)據(jù)集,區(qū)域A)Tab.4 Error comparison of the two inversion models (modeling dataset, area A)
從表4可以看出Stumpf 2003模型的決定系數(shù)為0.776; Lyzenga 1985模型的為0.860,并且均方根誤差也小于前者。所以在模型擬合方面,Lyzenga 1985反演模型優(yōu)于Stumpf 2003模型。
圖7 兩種模型的回歸標(biāo)準(zhǔn)化殘差圖 (建模數(shù)據(jù)集,區(qū)域A)Fig.7 Regression standardized residuals of the two models(modeling dataset, area A)
對于區(qū)域B,從數(shù)據(jù)集中選擇了100個樣本用于建模,25個樣本用于模型檢驗(yàn)。
迭代過程如表5和表6所示,與區(qū)域A結(jié)果類似,隨著迭代次數(shù)的增加,殘差平方和逐漸減小,模型中的參數(shù)逐漸趨向最優(yōu)值。
表5 Stumpf 2003模型的回歸迭代結(jié)果(區(qū)域B)Tab.5 Regression iteration results of Stumpf 2003 model(area B)
表6 Lyzenga 1985模型的回歸迭代結(jié)果(區(qū)域B)Tab.6 Regression iteration results of Lyzenga 1985 model(area B)
將Lyzenga 1985模型與Stumpf 2003模型的結(jié)果與真實(shí)水深進(jìn)行對比和分析,統(tǒng)計檢驗(yàn)點(diǎn)水深反演值的決定系數(shù)(R2)和均方根誤差(RMSE),結(jié)果如圖8和表7所示。
圖8 反演水深與實(shí)測水深對比(建模數(shù)據(jù)集,區(qū)域B)Fig.8 Comparison of the retrieved depth and the in-situ depth(modeling dataset, area B)
表7 兩種反演模型的誤差比較(建模數(shù)據(jù)集,區(qū)域B)Tab.7 Error comparison of the two inversion models (modeling dataset, area B)
從表7可以看出Stumpf 2003模型的決定系數(shù)為0.706,Lyzenga 1985模型的為0.878,二者均方根誤差相差不大。綜上,Lyzenga 1985模型擬合度仍然優(yōu)于Stumpf 2003模型。
圖9所示為由區(qū)域B數(shù)據(jù)集得到的兩種模型的回歸標(biāo)準(zhǔn)化殘差圖。同樣,兩種模型與高斯曲線吻合,殘差都基本符合正態(tài)分布。Stumpf 2003模型的標(biāo)準(zhǔn)偏差(1.427)顯著高于Lyzenga 1985模型(0.963),后者的殘差離散程度較小,分布更為集中??梢钥吹?,根據(jù)Stumpf 2003模型預(yù)測的水深與真實(shí)值的誤差集中于-2~1.5 m之間,而根據(jù)Lyzenga 1985模型預(yù)測的水深與真實(shí)值的誤差集中于-1~1.5 m之間。
圖9 兩種模型的回歸標(biāo)準(zhǔn)化殘差圖 (建模數(shù)據(jù)集,區(qū)域B)Fig.9 Regression standardized residuals of the two models(modeling dataset, area B)
本節(jié)利用未參與建模的數(shù)據(jù)集對4.1與4.2節(jié)中建立的模型進(jìn)行驗(yàn)證。
4.3.1 區(qū)域A的驗(yàn)證
利用未參與建模的110個樣本的數(shù)據(jù)集對區(qū)域A模型的驗(yàn)證結(jié)果如圖10和表8所示。
圖10 反演水深與實(shí)測水深對比(驗(yàn)證數(shù)據(jù)集,區(qū)域A)Fig.10 Comparison of the retrieved depth and the in-situ depth(validation dataset, area A)
表8 兩種反演模型的誤差比較(驗(yàn)證數(shù)據(jù)集,區(qū)域A)Tab.8 Error comparison of the two inversion models (validation dataset, area A)
由驗(yàn)證結(jié)果可知,Lyzenga 1985模型水深反演的決定系數(shù)R2和均方根誤差RMSE均優(yōu)于Stumpf 2003模型。
4.3.2 區(qū)域B的驗(yàn)證
利用未參與建模的25個樣本的數(shù)據(jù)集對區(qū)域B模型的驗(yàn)證結(jié)果如圖11和表9所示。
圖11 反演水深與實(shí)測水深對比(驗(yàn)證數(shù)據(jù)集,區(qū)域B)Fig.11 Comparison of the retrieved depth and the in-situ depth(validation dataset, area B)
表9 兩種反演模型的誤差比較(驗(yàn)證數(shù)據(jù)集,區(qū)域B)Tab.9 Error comparison of the two inversion models (validation dataset, area B)
由驗(yàn)證結(jié)果可知,Lyzenga 1985模型水深反演的決定系數(shù)R2較為顯著地優(yōu)于Stumpf 2003模型,均方根誤差RMSE小于Stumpf 2003模型。
結(jié)合兩種底質(zhì)區(qū)域的驗(yàn)證結(jié)果可知,在水體清澈的不同底質(zhì)區(qū)域,Lyzenga 1985模型普適性比Stumpf 2003模型更強(qiáng),能夠呈現(xiàn)出較為穩(wěn)定的反演效果。
本文利用高分辨率的WorldView-2影像反演了南海北島附近淺海水深。采用了Lyzenga 1985和Stumpf 2003兩種水深反演模型,并利用L-M算法求解模型參數(shù)。針對兩種不同的底質(zhì)進(jìn)行了水深反演與對比分析,根據(jù)A區(qū)(珊瑚底質(zhì))和B區(qū)(砂質(zhì)底質(zhì))的驗(yàn)證數(shù)據(jù),得到結(jié)論如下:
(1)由決定系數(shù)R2可知:對于A區(qū)和B區(qū),Lyzenga 1985模型的決定系數(shù)(0.902,0.897)均優(yōu)于Stumpf 2003模型(0.882,0.779)。
(2)由均方根誤差RMSE可知:對于A區(qū)和B區(qū),Lyzenga 1985模型的均方根誤差(1.651,0.529)均優(yōu)于Stumpf 2003模型(6.421,0.723)。說明利用Lyzenga 1985模型反演的水深值與實(shí)測值之間偏差較小。
綜上,在水體清澈的不同底質(zhì)區(qū)域,Lyzenga 1985模型普適性比Stumpf 2003模型更強(qiáng),能夠呈現(xiàn)出較為穩(wěn)定的反演效果。