基于改進黑暗像元法的環境一號CCD數據海岸帶區域大氣校正研究＊

李 微，郭錫杰，劉 遠，2，李媛媛，牟 蒙，劉長發，2

（1.大連海洋大學 海洋科技與環境學院，遼寧 大連116023；2.遼寧省高校近岸帶海洋環境科學與技術重點實驗室，遼寧 大連116023）

環境一號（HJ－1）衛星是我國自主研發的專門用于環境與災害監測的新型衛星，具有高時間分辨率、高空間分辨率和覆蓋范圍廣等優勢，其不僅為環境與減災業務運行系統提供重要保障，還成為很多部門日常業務的重要數據源。由于受大氣的影響，衛星數據不能真實地反映地物的情況，制約了利用衛星數據進行地面參數反演和地面監測的應用。因此，對衛星數據進行高精度的大氣校正是十分必要的，大氣校正的精度決定了后續遙感衛星數據地面參數反演的精度。

目前，針對衛星遙感數據的大氣校正算法分為2類：1）基于輻射傳輸原理建立起來的大氣校正輻射傳輸模型，該算法精度較高，包括6S模型、MODTRAN模型、FLAASH模型等，但該類算法模型計算量大，需要很多不易獲得的地面同步參數；2）不需要大氣和地面的實測數據的黑暗像元法，該方法主要通過衛星數據本身的信息就可完成，校正精度能滿足一般的遙感預處理和研究，應用于只有衛星遙感數據，而缺乏必要大氣和地面參數數據的情況。

近幾年，一些學者針對HJ－1衛星CCD數據的大氣校正算法進行了研究。金鑫等［1］利用HJ－1A衛星CCD數據提出一種基于輻射傳輸模型的大氣校正方法，對太湖水域大氣校正結果較好。孫長奎等［2］利用查找表方法來實現HJ－1ACCD數據的大氣校正，得到了研究區域的高空間分辨率的地表反射率數據。方莉等［3］運用像元純凈指數來提取CCD數據上的純像元，并由 HJ－1A星和B星的多時相CCD觀測數據結合地表雙向反射率模型，確定了純像元的地表反射特性，進而較大程度地提高了HJ－1ACCD數據的精度。同時，部分學者也深入研究了適用于HJ－1衛星的黑暗像元大氣校正算法。針對HJ－1ACCD數據，王中挺等［4］將密集植被作為暗像元，提出一種基于密集植被的暗目標法，利用輻射傳輸模型構建查找表，獲取可見光范圍的地物反射率進而求算氣溶膠光學厚度。劉其悅［5］以HJ－1ACCD數據為例，深入探討了基于輻射傳輸模型大氣校正和基于黑暗像元法的大氣校正在不同下墊面條件下應用的可行性，發現不同氣候和土地類型對模型的選擇有較大影響。綜上所述，不同專家學者針對特定區域利用HJ－1CCD數據分別進行了不同大氣校正算法的研究，但針對海岸帶區域的研究尚不多見。

因此，本文以海岸帶為研究區域，針對HJ－1CCD數據，提出了一種改進的黑暗像元大氣校正算法。筆者以雙臺子河口區域為例，利用黑暗像元法大氣校正的原理，首先結合歸一化植被指數（Normalized Difference Vegetation Index，NDVI）和改進的歸一化差異水體指數（Improved Normal Differential Water Index，INDWI）確定初始黑暗像元區域，再利用區域增長法對該區域進行優化并確定最終黑暗像元，然后計算大氣的影響因素程輻射，進而去除該影響獲得大氣校正后的地表反射率R。最后，將該算法與真實地物光譜以及ENVI軟件下的Dark－subtract模塊和FLAASH大氣校正模塊的校正結果進行對比分析。

1 研究區域概況

我國海域面積大，海岸帶長，近海岸帶有著豐富的自然資源以及生物多樣性，因此海岸帶的監測和保護顯得尤為重要。近海岸區域存在大量水體，而且有的區域存在濃密的植被，如紅樹林、蘆葦或翅堿蓬等，另外，部分區域還會存在小島陰影，而水體、濃密的植被以及陰影在衛星遙感數據上都是理想的黑暗像元。

本文以遼寧雙臺子河口為研究區域，該區域位于遼寧省盤錦市境內雙臺子河口國家級自然保護區河口濕地，雙臺子河入海口，是中國高緯度地區面積最大的濱海濕地。該地地貌為沖海積平原，海岸帶地勢低洼且平坦，同時分布著大面積的蘆葦沼澤區、翅堿蓬灘涂和淺海海域，海岸帶較長。因此該區域存在大片水體和濃密的植被，這為黑暗像元的選取提供了條件。

2 大氣校正過程

遙感成像過程中，由于大氣的影響遙感影像存在一定的輻射量失真現象，降低了遙感數據的精度，因此獲得衛星遙感數據后必須進行大氣校正。大氣校正的目的是消除大氣和光照等因素對地物反射率的影響，以獲得地物的真實反射率。

2.1 輻射定標

大氣校正前，為了建立遙感傳感器的數字量化輸出值DN與其所對應視場中輻射亮度值之間的定量關系，首先需要對數據進行輻射定標，將DN值轉化為表觀輻射度L。本文采用中國資源衛星應用中心2013－02－28公布的方法和定標系數（表1）對HJ－1衛星進行輻射定標［6］，其公式如下：

式中，A為定標系數增益，單位為DN／（W·m－2·sr－1·μm－1）；L0為定標系數偏移量，單位為 W·m－2·sr－1·μm－1。

表1 HJ－1衛星CCD數據的定標系數Table 1 Radiometric calibration coefficient of HJ－1CCD data

2.2 大氣校正

本文提出的算法有如下假定：1）待校正區域存在黑暗像元；2）地表為朗伯面反射和大氣性質均一；3）大氣多次散射輻照作用和鄰近像元漫反射作用可以忽略。因此反射率小的黑暗像元在大氣的影響下，反射率會相對增加［7］。

地表反射率R可通過式（2）計算：

式中，L為真實地物輻亮度，可由公式（1）得到；LP為程輻射，E0為大氣層外太陽光譜輻照度，參見表1；ED為天空光漫射到地表面的光譜輻照度，可利用Moran［8］的方法，即用πLP估算；θ為太陽天頂角；φ為傳感器天頂角，均可由遙感圖像頭文件中獲取；Tθ和Tφ為大氣透過率，從可見光到短波紅外過程中，當θ，φ小于70°時，圖像受到大氣散射和弱吸收的影響，大氣透過率Tθ和Tφ可由以下兩式估算［9］：

式中，τ為光學厚度，針對近海岸帶大氣環境，依據Kaufman［10］提出的簡化計算方法，利用式（5）估計得出τ：

式中，λ為波長，單位為μm。

由上可知，去除大氣的影響獲得真實的地表反射率，關鍵就是獲得程輻射，該值可以通過黑暗像元值來確定。因此準確、快速地提取黑暗像元值是黑暗像元法大氣校正的核心內容。本文提出的改進黑暗像元大氣校正算法具體步驟如下。

2.2.1 初始黑暗像元區域的確定

筆者利用歸一化植被指數（NDVI）、改進的歸一化差異水體指數（INDWI）對初始黑暗像元區域進行提取。

1）濃密植被備選區域確定

利用歸一化植被指數來確定植被備選區域。NDVI是植被生長狀態及植被覆蓋的最佳指示因子，NDVI值越大，植被覆蓋或長勢越好。NDVI表達式如下：

式中，Lband4和Lband3分別為HJ－1CCD數據的近紅外波段和紅光波段的表觀輻射度，其計算公式參見式（1）。

研究區域沒有濃密的林地，但是存在成片、長勢茂密的葦田和生長密集的翅堿蓬群落。通過大量試驗最終確定，當NDVI＞0.5時可以提取出濃密的蘆葦和翅堿蓬，并將其作為黑暗像元的備選區域。

2）純凈水體區域確定

通常認為，純凈的水體可以作為黑暗像元。但近岸水體中含有大量泥沙以及葉綠素，屬于二類水體，因此如何提取到純凈水體是研究重點。相關研究證明，采用歸一化的比值指數運算可以最大程度地抑制植被、城市和山地等造成的噪聲，從而達到突出水體信息的目的［11］。本文通過研究提出了改進的歸一化差異水體指數（INDWI）組合，計算公式如下：

式中，Lband1，Lband3和Lband4分別為HJ－1CCD數據的藍波段、紅波段和近紅外波段的表觀輻射度，其計算公式參見式（1）。

通過大量試驗最終確定，當INDWI1＞0.1且INDWI2＞0.2時即可認為提取到較純凈的水體，作為黑暗像元的備選區域。

2.2.2 黑暗像元區域優化與提取

區域增長法的基本思想是將具有相似性質的像元整合成為一個區域。本文利用區域增長法對黑暗像元區域進行優化提取，具體算法參見李微等［12］。具體實現過程如下：1）確定初始黑暗像元生長點：即將黑暗像元備選區域中第一個使LP大于零的非零像元值所對應的像元作為生長點，可以是一個或者多個生長點；2）確定生長規則：將每個生長點8鄰域中值和標準差以及整幅圖像平均值作為特征值，可降低隨機噪聲的影響，避免中值過大的情況；3）進行多閾值區域增長：采用8鄰域增長，當符合生長規則時，將其作為新的生長點并更新中值和標準差，直至無法生長為止。如果區域生長范圍過小，則認定原生長點是噪聲點，同時為了節省運算時間，當區域增長范圍過大時終止運算，并將該范圍作為黑暗像元區。

統計每一個波段所有增長區域的像素平均值和增長區域個數，將所有增長區域像素平均值相加，除以增長區域個數，即獲得黑暗像元值DW′。

2.2.3 程輻射LP的計算

由于程輻射的影響，在傳感器上黑暗像元反射率很小但并不為零，由此可以計算出程輻射。一些研究學者例如 Kaufman等［13］假設黑暗像元的反射率為（0.02±0.01），還有一些學者［8，14］假設黑暗像元的反射率為0.01。本文結合研究區域現狀，假設黑暗像元具有1%的地表反射率，因此程輻射LP計算方法如下：

式中，L′為黑暗像元區域的輻亮度，計算公式如式（1）所示。將程輻射LP帶入到式（2）中，即可得到真實的地表反射率R。

3 大氣校正結果與分析

本文利用成像日期為2013－09－01的HJ－1ACCD2級數據為例，驗證算法的有效性。研究區域原始數據真彩色合成影像顏色較暗，圖1為2%線性拉伸后研究區域真彩色合成影像，由圖可見圖像對比度差。

圖1 2013－09－01研究區域 HJ－1ACCD 3，2，1波段真彩色合成影像Fig.1 True color composite image of HJ－1ACCD with bands 3，2and 1in study area on September 1，2013

現以該數據的第三波段為例說明黑暗像元區域優化提取的過程（圖2）。其中：圖2a為利用NDVI和INDWI確定的黑暗像元備選區域，圖中白色表示黑暗像元備選區域；圖2b為第三波段數據初始黑暗像元生長點的位置，用紅色標示；圖2c為區域增長優化后黑暗像元最終分布圖，用紅色標示。

圖2 黑暗像元區域提取過程圖Fig.2 The schematic drawing of extracted dark pixel areas

圖3為3種大氣校正算法校正后的真彩色合成圖方法，其中圖3a和3b分別為ENVI軟件下的黑暗像元法（Dark－Subtract）和FLAASH方法的校正結果，圖3c為改進的黑暗像元大氣校正算法（IMDS）的校正結果。三者與圖1對比可見圖像對比度均明顯提高，其中，IMDS方法目視效果更好。

圖3 大氣校正以后3，2，1波段合成真彩色影像圖Fig.3 The true color image with bands 3，2，1after atmospheric correction

為了驗證和定量評價改進的黑暗像元大氣校正算法，本文對不同大氣校正方法校正前后的研究區域內4種典型地物（蘆葦、翅堿蓬、水體和灘涂）光譜曲線以及地面實測光譜進行分析對比，實測光譜利用美國ASD地物光譜儀測得的光譜，按照CCD波段帶寬取平均得到，如圖4所示。從圖中可以看出，相比大氣校正前，IMDS、Dark－Subtract法和FLAASH法在一定程度上均可以消除大氣的影響，相同地物的光譜曲線趨勢大致相同，并且相比FLAASH方法，IMDS方法與地面實測值更接近，效果較好。

由圖4a可知，IMDS和FLAASH法得到的蘆葦光譜曲線在綠光波段（0.55μm）有一小反射峰，紅光波段（0.68μm）有一光合作用吸收谷［15］，曲線趨勢相近，而Dark－subtract法在綠光波段沒有出現明顯的反射峰，效果較差；影像獲取時值9月，翅堿蓬群落呈現紫紅色，從圖4b可見，3種方法獲得光譜曲線在紅光波段反射率都較高，在綠光波段也沒有出現反射峰［16］，但IMDS與FLAASH法更接近；圖4c為大氣校正前后水體的光譜曲線，由圖中可見，IMDS和FLAASH法得到水體光譜曲線在綠光波段出現小的反射峰，之后反射率逐漸下降［17］，Dark－subtract法沒有出現明顯的反射峰；圖4d中IMDS和FLAASH法得到的土壤的光譜曲線均優于Dark－subtract法。總的來說，IMDS更接近基于輻射傳輸模型的FLAASH法，與真實地物實測光譜特征擬合較好，精度較高。

圖4 不同方法大氣校正前后典型地物光譜與實測光譜值對比Fig.4 Comparison between in situ－measured spectra and satellite－derived ones by different atmospheric correction methods

地物反射率的正常范圍為0～1，而大氣校正后會出現負值，我們利用“數據可利用率”表示非負值像元個數所占的比例。IMDS與FLAASH法對HJ－1CCD數據可利用率對比見表2。由表2可知，IMDS能有效的改善這種現象，提高了HJ－1CCD數據的可利用率。

表2 大氣校正后IMDS與FLAASH數據可利用率（%）對比Table 2 The usable percentage（%）comparison between IMDS and FLAASH after atmospheric correction

本文在計算大氣透過率時，IMDS算法只是考慮了大氣分子的衰減（式（5）），沒有考慮氣溶膠的衰減作用；在估算地表下行天空光輻照度時，也只假設了天空光為完全漫射光，沒有考慮天空光存在顯著的各向異性。這是該算法的不足之處，也是今后需要改進的重點。

4 結論

本文基于黑暗像元大氣校正模型，針對海岸帶的地表狀況，提出了一種改進的黑暗像元大氣校正模型。首先結合歸一化植被指數（NDVI）、改進的歸一化差異水體指數（INDWI）來確定黑暗像元，同時利用區域增長法對黑暗像元區域進行優化，計算出程輻射，實現了大氣校正過程。結果表明，IMDS得到地物反射率曲線相近，與真實地物光譜特征擬合較好，且進行大氣校正之后的圖像地物邊緣清晰度明顯增強，圖像亮度增加，圖像對比度較高，目視效果好。同時與FLAASH法相比，IMDS能改善FLAASH法校正過程中出現負值的現象，可以有效的保證了HJ－1CCD數據的可利用率。

致謝：文中HJ－1A數據來源于國家環境保護部衛星環境應用中心（http：∥www.secmep.cn），在此表示感謝！

（References）：

［1］JIN X，LI Y M，WANG Q，et al.Atmospheric correction of satellite HJ－1CCD image based on aerosol distribution in Lake Taihu［J］.Journal of Lake Sciences，2010，22（4）：504－512.金鑫，李云梅，王橋，等.基于太湖氣溶膠類型分區的環境一號衛星 CCD大氣校正［J］.湖泊科學，2010，22（4）：504－512.

［2］SUN C K，SUN L，MA S F，et al.Atmospheric correction method based on HJ－1CCD data［J］.Journal of Remote Sensing，2012，16（4）：831－836.孫長奎，孫林，麻盛芳，等.HJ－1CCD數據大氣校正方法研究［J］.遙感學報，2012，16（4）：831－836.

［3］FANG L，YU T，GU X F，et al.Aerosol retrieval and atmospheric correction of HJ－1CCD data over Beijing［J］.Journal of Remote Sensing，2013，17（1）：158－164.方莉，余濤，顧行發，等.北京地區 HJ－1衛星CCD數據的氣溶膠反演及在大氣校正中的應用［J］.遙感學報，2013，17（1）：158－164.

［4］WANG Z T，LI Q，TAO J H，et al.Monitoring of aerosol optical depth over land surface using CCD camera on HJ－1satellite［J］.China Environmental Science，2009，29（9）：902－907.王中挺，厲青，陶金花，等.環境一號衛星 CCD 相機應用于陸地氣溶膠的監測［J］.中國環境科學，2009，29（9）：902－907.

［5］LIU Q Y.The atmospheric correction for visible－neainfrared region of HJ－1constellation data［D］.Henan：Henan Polytechnic University，2010.劉其悅.基于環境星數據可見光－近紅外波段的大氣校正方法研究［D］.河南：河南理工大學，2010.

［6］Chinese Center for Resources Satellite Data and Application.HJ－1AB，ZY－3and 02Cabsolute radiometric calibration coefficient in 2013.［EB／OL］.［2013－09－17］.http：∥www.cresda.com／n16／n1115／n1522／n2103／index.html.中國資源衛星應用中心.2013年 HJ－1AB，ZY－3和02C絕對輻射定標系數［EB／OL］.［2013－09－17］.http：∥www.cresda.com／n16／n1115／n1522／n2103／index.html.

［7］DENG S B.ENVI RS image processing methods［M］.Beijing：Science Press，2010.鄧書斌.ENVI遙感圖像處理方法［M］.北京：科學出版社，2010.

［8］MORAN M S，JACKSON R D，SLATER P.Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output［J］.Remote Sensing of Environment，1992，41（2）：169－184.

［9］ZHENG W，ZENG Z Y.Dark－object methods for atmospheric correction of remote sensing image［J］.Remote Sensing for Land ＆ Resources，2005，63（3）：8－11.鄭偉，曾志遠.遙感圖像大氣校正的黑暗像元法［J］.國土資源遙感，2005，63（3）：8－11.

［10］KAUFMAN Y J，WALD A E，REMER L A，et al.The MODIS 2.1μm channel－correlation with visible reflectance for use in remote sensing of aerosol［J］.IEEE Transactions on Geoscience and Remote Sensing，1997，35（5）：1286－1296.

［11］XU H Q.A study on information extraction of water body with the modified normalized difference water index（MNDWI）［J］.Journal of Remote Sensing，2005，9（5）：589－595.徐涵秋.利用改進的歸一化差異水體指數（GNDWI）提取水體信息的研究［J］.遙感學報，2005，9（5）：589－595.

［12］LI W，TIAN Y，LIU Y，et al.The algorithm to extract automatically dark objects in coastal zone based on region growing［J］.Journal of Dalian Ocean University，2013，28（5）：502－505.李微，田彥，劉遠，等.基于區域增長的海岸帶黑暗像元自動提取算法研究［J］.大連海洋大學學報，2013，28（5）：502－505.

［13］KAUFMAN Y J，SENDRA C.Algorithm for automatic atmospheric corrections to visible and near－IR satellite imagery［J］.International Journal of Remote Sensing，1988，9（8）：1357－1381.

［14］CHAVEZ Jr P S.Image－based atmospheric correction revisited and improved［J］.Photogrammetric Engineering and Remote Sensing，1996，62：1025－1035.

［15］MEI A X，PENG W L，QIN Q M，et al.An introduction to remote sensing［M］.Beijing：High Education Press，2001.梅安新，彭望琭，秦其明，等.遙感導論［M］.北京：高等教育出版社，2001.

［16］WU T，ZHAO D Z，KANG J C，et al.Research on remote sensing inversion biomass method based on the Suaeda Salsa＇s Measured Spectrum［J］.Spectroscopy and Spectral Analysis，2010，30（5）：1336－1341.吳濤，趙冬至，康建成，等.基于現場光譜的翅堿蓬生物量遙感反演方法研究［J］.光譜學與光譜分析，2010，30（5）：1336－1341.

［17］ZHANG Y，ZHANG Y，WANG J J.Correlations between the spectral reflectance of water mass containing suspended sediments and the concentration and particle size of the suspended sediments［J］.Advances in Marine Science，2008，26（3）：340－346.張蕓，張鷹，王晶晶.懸沙水體光譜反射率與質量濃度、粒徑的相關關系［J］.海洋科學進展，2008，26（3）：340－346.