資源一號02D衛星高光譜數據水體透明度反演研究

周毅 劉瑤 田淑芳

(1 中國地質大學(北京) 地球科學與資源學院，北京 100083)(2 自然資源部國土衛星遙感應用中心，北京 100048)

透明度作為反映水體清澈/渾濁程度的水質參數，被廣泛應用于湖泊學和海洋學研究中，對評估水體內水生植物多樣性、生產力和水體營養程度等方面具有重要意義[1]。過去幾十年來，受人口增長、經濟發展中不合理的開發利用及污染等因素影響，我國許多內陸水體(湖泊、河流和水庫)呈富營養化狀態，水體質量急劇下降，帶來水體生態系統的結構破壞與功能退化，對人民生活及國家經濟造成巨大損失和不良影響[2]。近年來，國家為改善水生態環境，開展了重點流域保護修復、重點工程水質保障等一系列工作。有必要定期開展內陸水體透明度等水質參數的監測工作，服務于生態修復措施的水環境改善效果評價，以及為后續水資源的保護和管理提供決策數據支持。透明度通常定義為塞氏盤(Secchi Disk)在水中垂直下放至剛剛在視野中消失的深度。由于塞氏盤獲取成本低且操作簡便，使用塞氏盤進行水體透明度測量已有上百年的歷史。然而，通過布設地面監測點或者船載測量的方式難以實現對大區域和多湖泊的透明度重復觀測。相比而言，能夠周期性穩定獲取和實現空間覆蓋的遙感衛星，成為了開展長期大規模水體透明度反演的重要數據源。

國內外已有眾多學者基于遙感數據開展水體透明度反演研究，常見方法有經驗/半經驗模型和半分析模型。文獻[3]基于陸地衛星-5(Landsat-5)專題成像儀(TM)數據的藍波段、紅藍波段反射率比值，使用相關系數矩陣、多元逐步回歸等方法完成了1973～1998年間美國大型湖泊的透明度評價研究；文獻[4]使用Landsat-5和陸地衛星-7(Landsat-7)增強型專題成像儀(ETM+)數據的藍紅波段，通過建立經驗模型對1990～2010年間美國緬因州湖泊開展透明度反演研究；文獻[5-6]建立了一種基于半分析算法(Quasi Analytical Algorithm，QAA)的水體透明度反演模型，可應用在包括中分辨率成像光譜儀(MODIS)、中等分辨率成像頻譜儀(MERIS)、陸地衛星-8(Landsat-8)等衛星的多種傳感器上實現水體透明度的反演。國內方面，文獻[7]使用Landsat-5中紅藍波段建立自然對數與透明度自然對數的線性模型，對2004年6個時期鄱陽湖的水體開展透明度反演，取得了較好精度；文獻[8]應用Landsat-8 陸地成像儀(OLI)數據的紅、綠波段建立對數模型，反演了2013年的東平湖水體透明度，取得了19.77%的總體相對誤差；文獻[9]基于透明度和水體漫衰減系數關系建立半分析模型，并通過MODIS、近海高光譜成像儀(HICO)數據對近岸水體透明度進行反演。

目前，水體透明度反演研究多以多光譜數據為數據源，采用星載高光譜數據開展研究相對較少。相較于多光譜數據，高光譜數據具有更高的光譜分辨率，在用于水體透明度反演時，能夠提供更豐富的波段選擇，且具備更精準的光譜特征，因此在水質參數反演精度的提升方面具有優勢和潛力。

我國于2018年5月9日成功發射高分5號(GF-5)衛星，其搭載新一代高光譜相機(Advanced Hyper Spectral Imager，AHSI)的光譜覆蓋范圍為0.40～2.50 μm，空間分辨率為30 m，可見光近紅外(VNIR)與短波紅外(SWIR)分別有150和180個波段，光譜分辨率分別為5 nm和10 nm。此外，資源一號02D衛星(又稱為5米光學業務衛星)也于2019年9月12日成功升空。該星搭載的AHSI，具有與GF-5衛星相同的光譜覆蓋范圍和空間分辨率，可見光近紅外波段有76個通道(光譜分辨率為10 nm)，短波紅外波段包括90個通道(光譜分辨率為20 nm)。我國后續還將發射搭載高光譜相機的衛星，通過多星組網，未來有望能夠基于星載高光譜數據開展內陸水體水質監測。

本文以北京市官廳水庫為研究區，以資源一號02D衛星高光譜影像為數據源，應用半經驗和半分析模型方法反演研究區的水體透明度。通過對結果精度進行評價，分析資源一號02D衛星高光譜數據在反演水體透明度方面的表現與應用潛力。

1 材料與方法

1.1 研究區概況

官廳水庫位于北京西北部永定河上游(見圖1)，水庫面積約130 km2。1970年以來，受入庫流量下降、污染負荷增加等影響，官廳水庫水質不斷惡化，于1997年退出北京飲用水水源地之列；后經治理，水庫在2007年再度成為北京備用飲用水水源地。近年來，隨永定河綜合治理和生態修復推進，官廳水庫水質逐年好轉[10]。

圖1 研究區示意圖Fig.1 Map of study area

1.2 數據及處理

2020年8月13日在官廳水庫開展了與資源一號02D衛星同步的地表試驗，獲取了地表實地的光譜數據和水質參數(葉綠素a、懸浮物濃度、透明度等)。高光譜遙感數據的成像時間為2020年8月13日上午11時21分，地面試驗則在當天10～14時內完成，以保證星地數據獲取的一致性。本次試驗在研究區范圍內布設了18個均勻分布的點位(見圖2)，其中實測水體透明度范圍為0.4～0.7 m，均值為0.58 m。

圖2 水體透明度實測位置圖Fig.2 Location map of measured water transparency

對獲取的星載高光譜數據，進行了輻射定標、大氣校正、水體范圍提取和遙感反射率(Rrs)計算等一系列預處理工作，具體的處理過程如圖3所示。

圖3 高光譜數據處理過程Fig.3 Hyperspectral data process flow

1)輻射定標

通過ENVI5.3軟件，將原始影像的VNIR數據(76個波段)和SWIR數據(90個波段)合并為一個文件，應用Apply Gain and Offset工具導入絕對定標參數后進行輻射定標。

2)大氣校正

FLAASH大氣校正模型是目前在多/高光譜影像中使用效果最好的大氣校正工具之一，對于地物波譜信息可以高保真地恢復，簡潔快速地獲取比較準確的物理模型參數[11]。因此本文在ENVI5.3中應用FLAASH模塊對影像數據進行大氣校正。綜合考慮數據特點、成像時間及研究區位置，本文在FLAASH模塊中選擇中緯度夏季(Mid-Latitude Summer)大氣模型，水體反演波段設置為940 nm，并選用鄉村(Rural)氣溶膠模型，以及2-波段(考夫曼-泰倫)(2-Band(K-T))氣溶膠反演方法。

3)水體范圍提取

本文選擇歸一化水體指數(Normalized Difference Water Index，NDWI)方法進行水體范圍提取[12]，其公式為

(1)

式中：ρGreen、ρNIR分別代表綠波段和近紅外波段的地表反射率。

式(1)最初以Landsat TM數據進行計算，而資源一號02D衛星高光譜數據信噪比相對較低，選擇單一波段計算容易受到噪聲的影響[13]。因此，本文參照Landsat TM數據的綠光(0.52～0.60 μm)和近紅外波段范圍(0.76～0.90 μm)，通過波段合并，分別計算資源一號02D衛星數據在Landsat TM這兩個波段范圍內的地表反射率平均值，用于歸一化水體指數的計算。在此基礎上，運用直方圖雙峰法(Histogram Bimodal Method，HBM)分離水體和其他類型地物，并將小于1000像元(約0.9 km2)的微小水體消除，只保留主體水域，從而實現影像中水體范圍的提取。

4)遙感反射率計算

遙感反射率計算公式為

(2)

式中：λ為波長，ρ(λ)為經過FLAASH大氣校正后的地表反射率，ρsky(λ)為天空光反射率，Lu(λ)為水面上行輻亮度，Ld(λ)為水面下行輻亮度，Lsky(λ)為天空光輻亮度。

即使是反射率相對較高的渾濁水體，在短波紅外的反射率也趨于0，因此通常認為圖像中短波紅外的反射率是來自天空光的影響，將短波紅外的反射率用于天空光的去除[14]。但是，由于AHSI的輻射定標系數是基于地表亮目標計算得到，而水體反射率相對較低，造成水體像元在短波紅外波段噪聲相對較大，不適宜用于天空光的糾正。因此，本文采用忽略天空光的影響，以ρ(λ)/π作為遙感反射率的近似值。盡管采用這樣近似計算的結果可能略大于實際的遙感反射率，但由于通常認為天空光的影響與波長無關，因此，近似計算的遙感反射率與真實反射率具有一致的光譜形狀?；诖?，本文采用波段比值的經驗模型，和半分析模型進行透明度反演。這兩類方法對于透明度的反演，主要取決于光譜的波形特征，受遙感反射率的量級影響較小[15]。

1.3 水體透明度反演方法

本文分別運用半經驗模型和半分析模型方法進行官廳水庫水體透明度反演。

1.3.1 半經驗模型

根據文獻[16]研究成果，半經驗模型方法主要通過對地面同步點位的影像光譜和地表測量的透明度進行相關性分析，選擇合適的波段或波段組合進行建模。本文通過分析資源一號02D衛星圖像和地表采集的水體光譜特征，選擇中心波長為653 nm的波段，采用對數回歸模型(見圖4)得到水體透明度反演公式為

Zsd=-0.438×ln (Rrs(653))-1.335 8

(3)

式中：Rrs(653)為中心波長為653 nm波段處的遙感反射率。

圖4 對數回歸模型Fig.4 Logarithmic regression model

1.3.2 半分析模型

通過半分析模型反演水體透明度主要步驟為：①估算固有光學量(IOPs)，即總吸收系數a和后向散射系數bb；②基于a和bb反演漫衰減系數(Kd，m-1)；③用Kd和Rrs確定Zsd。

本文分別運用第5版及第6版半分析算法(QAAv5、QAAv6)[17]計算總吸收系數a和后向散射系數bb；并采用更新后的半分析模型[18]對漫衰減系數Kd進行反演

m1×(1-m2×e-m3×a(λ))bb(λ)

(4)

式中：θs為太陽天頂角，bbw(λ)為純水的后向散射系數，m0-3與γ分別取0.005、4.26、0.52、10.8和0.265的固定值，a(λ)為總吸收系數，bb(λ)為后向散射系數。

對于水體透明度的反演，本文通過公式[19]進行計算

(5)

式中：Kd(tr)表示水體在可見光譜范圍(410～665 nm)上的漫衰減系數，Rrs(tr)表示該波長下相應的遙感反射率。

1.4 評價方法

本文通過建立線性回歸模型來分析反演結果與實測值之間的相關性，并計算絕對誤差(AE)、相對誤差(RE)以及均方根誤差(RMSE)用于對結果的精度評價為

IAE=|dm-dp|

(6)

(7)

(8)

式中：dm和dp分別表示透明度的實測值和反演值；n為樣本數。

2 結果與分析

2.1 不同反演結果與實測值對比

本文分別運用1.3節中的方法采用資源一號02D衛星數據對研究區水體透明度進行計算，并提取了與實測點位匹配的相應反演值。應用線性回歸模型對反演結果與實測值進行了相關性分析并得到相應的相關系數值(見圖5)。

圖5 線性回歸分析圖Fig.5 Linear regression analysis results

此外，運用1.4節中各誤差分析方法計算得到各反演結果的誤差分析情況(見表1)。

表1 反演結果誤差分析表Table 1 Error analysis of inversion results

根據線性回歸分析結果可以看出，運用半經驗模型方法反演得到的水體透明度結果與實測結果相關性最高，R2為0.473 3；而基于QAAv5與QAAv6的半分析模型方法反演結果相關性較低，其中基于QAAv5計算的反演結果與實測值相關性最低，R2為0.419 3。圖5(b)和(c)中，由于QAAv5與QAAv6算法主要是針對一類大洋水體提出，是以550～560 nm作為標準波長；將其應用于渾濁水體，會造成對總吸收系數和后向散射系數的低估，進而低估漫衰減系數，導致反演的水體透明度偏高[20-21]，因此擬合線低于1∶1線。通過表1可知，半經驗模型取得的反演結果，相較2種半分析模型方法結果，整體精度較高，絕對誤差平均值為0.05 m，平均相對誤差為9.5%，而RMSE僅為0.07 m；兩種半分析模型得到的反演結果精度相對較低，其中，基于QAAv5計算的反演結果精度較差，最大相對誤差達到了174.0%，其平均相對誤差為94.7%，RMSE為0.55 m。而基于QAAv6反演結果平均相對誤差及RMSE分別為31.4%、0.19 m，相較QAAv5結果精度稍高，主要原因是QAAv6算法對渾濁水體的情況進行了適當的優化，對渾濁的內陸水體的適用性稍有提升。

上述3種方法反演的透明度結果與實測結果的R2都比較低，這主要是由于本研究區的同步點位數量有限；另一方面，反演結果的精度分析可以看到，表明該方法在應用于官廳水庫透明度反演中具有較好的效果。

2.2 官廳水庫水體透明度分析

通過上述分析，運用半經驗模型方法反演得到的水體透明度具有較高精度，因此本文基于該方法得到官廳水庫水體透明度(Zsd)狀態圖(見圖6)。

圖6 官廳水庫Zsd狀態圖Fig.6 Zsd state of Guanting reservoir

從圖6可以看出，官廳水庫水體透明度范圍在0～1 m之間，相比東北部水體，西南部水體透明度更高，說明該區域水體質量更好。該結果與實地測量得到的水質狀態基本一致。

3 結論

本文以北京市官廳水庫為研究區，采用資源一號02D衛星高光譜數據，應用半經驗模型、半分析模型方法進行研究區的水體透明度反演。通過反演結果與實地測量結果的相關性分析和精度評價，得到主要結論如下：

(1)本文分別使用半經驗模型、基于QAAv5與QAAv6的半分析模型反演得到水體透明度，通過結果與實測數據相關性及精度評價，半經驗模型方法得到的反演值與實測值相關性最高，相關系數為0.473 3，且結果精度也最高，平均絕對誤差、平均相對誤差及RMSE分別為0.05 m、9.5%、0.07 m；其次為基于QAAv6半分析模型方法的反演結果，平均絕對誤差、平均相對誤差及RMSE分別為0.17 m、31.4%、0.19 m；而基于QAAv5的半分析模型反演結果精度較差。

(2)基于資源一號02D衛星的高光譜數據，可實現水體范圍的提取并用于內陸渾濁水體的透明度反演。相比較常規的水色衛星，資源一號02D衛星具有更高的空間分辨率，可應用于中小型湖庫的水質監測。下一步的研究是驗證該數據應用于更多不同光學特性水體中的透明度反演能力。