遥感图像大气校正的黑暗像元法

书书书第１期，总第６３期国　土　资　源　遥　感Ｎｏ．１，２００５　２００５年３月１５日ＲＥＭＯＴＥＳＥＮＳＩＮＧＦＯＲＬＡＮＤ＆ＲＥＳＯＵＲＣＥＳＭａｒ．，２００５　遥感图像大气校正的黑暗像元法郑　伟，曾志远（南京师范大学地理科学学院，南京　２１００９７）摘要：介绍黑暗像元法大气校正的思路，阐述了不同黑暗像元法的大气校正模型，选择其中２种模型对一景ＳＰＯＴ图像进行了大气校正，最后指出了黑暗像元法需要进一步改进的地方。关键词：遥感图像；黑暗像元法；大气校正中图分类号：ＴＰ７５１　文献标识码：Ａ　　　文章编号：１００１－０７０Ｘ（２００５）０１－０００８－０４０　引言卫星遥感图像大气校正方法很多，校正精度较高的方法是利用基于复杂的辐射传输原理建立起来的大气校正模型，如６Ｓ模型、ＬＯＷＴＲＡＮ模型及ＭＯＲＴＲＡＮ模型等，但是，这些方法计算量大，而且需要参数多，而这些参数的获得对于大多数卫星图像应用人员来说是难以实现的。因此，研究者提出了一些不需大气和地面的实测数据，尤其不需要卫星同步观测数据的大气校正方法，其中被广泛应用的就是黑暗像元法。用黑暗像元法进行大气校正主要是依靠图像本身的信息，一些不能直接在图像上获得的信息，可在相应的文献资料中找到。这种方法直接、简易，其校正精度可以满足一般遥感研究和应用，具有较强的实用性。该方法的基本原理就是在假定待校正的遥感图像上存在黑暗像元区域，地表为朗伯面反射，大气性质均一，大气多次散射辐照作用和邻近像元漫反射作用可以忽略的前提下，反射率或辐射亮度很小的黑暗像元由于大气的影响，而使得这些像元的亮度值相对增加，可以认为这部分增加的亮度是由于大气的程辐射影响产生的。利用黑暗像元值计算出程辐射，并代入适当的大气校正模型，获得相应的参数后，通过计算就得到了地物真实的反射率。１　大气校正简化理论模型假设天空辐照度各向同性，地表面是一个理想的朗伯体，并忽略大气的折射、湍流和偏振，由遥感方程，地物表面反射辐亮度为Ｌ０＝（Ｒ／π）（Ｅ０ｃｏｓθＴθ＋ＥＤ）（１）　　式中，Ｒ为地物表面反射率；Ｅ０为大气层外相应波长的太阳光谱辐照度；θ为太阳天顶角；Ｔθ为在太阳辐射入射方向上的大气透射率；ＥＤ为由天空光漫射到地表面的光谱辐照度。由卫星接收到的表观辐亮度是地面反射辐亮度、大气透射率、大气反照率和由于大气散射造成的向上大气光谱辐亮度Ｌｐ（即程辐射）的函数，即Ｌ＝Ｌ０Ｔ／（１－Ｓ·Ｒ）＋Ｌｐ（２）　　式中，Ｔ是地物到传感器反射方向的大气透射率；是卫星传感器天顶角；Ｓ是大气反照率。考虑到上式中Ｓ非常小，即Ｓ·Ｒ＝１［１］，因此，由（１）式和（２）式可得到大气校正模型为Ｒ＝π（Ｌ－Ｌｐ）／［Ｔ（ＴθＥ０ｃｏｓθ＋ＥＤ）］（３）　　在式（３）中，Ｌ可由星上或地面定标结果求得；θ由日期和时间计算得；可从遥感图像头文件中读出；Ｅ０可由探测器响应函数计算求得。因此，要求Ｒ，还有Ｌｐ、Ｔ、Ｔθ和ＥＤ４个未知数。而对Ｔ、Ｔθ、ＥＤ的不同简化假设［２，３］，以及对Ｌｐ、Ｔ、Ｔθ和ＥＤ的不同计算方法，得到了不同的黑暗像元法。近些年来，很多研究者对黑暗像元法的改进也是基于此。２　程辐射Ｌｐ的确定２．１　黑暗像元的确定遥感图像中的黑暗像元可以是云或山体形成的阴影，也可以是水体或某些类型的土壤（黑土），或是收稿日期：２００４－０８－３１；修订日期：２００４－１１－１５基金项目：国家自然科学基金项目（４０３７１０５３）。第１期郑伟，等：　遥感图像大气校正的黑暗像元法高密度的植被区（例如针叶林）等。Ａｈｅｒｎ等［４］和Ｇｏｒｄｏｎ［５］就采用了把深且清洁的水体像元作为黑暗像元的方法进行大气校正。这种方法被认为在海洋上具有可接受的精度。Ｃｈａｖｅｚ［６］根据遥感图像不同波段的直方图，特别是对于可见光波段，从零亮度值对应的像元到非零亮度值对应的像元数量有一个很快的上升，像元数量突然增加处所对应的亮度值就可以当作黑暗像元值；或者是图像直方图由非零亮度值开始，逐渐随像元个数增高，则非零亮度值开始处对应的亮度值就可以当做黑暗像元值。曾志远［７］指出，由于水体辐射性质与一般地物不同，用它的最低值来校正大多数地物，可能并不理想。而直方图最低值法中的最低值，也是水体和其它低辐射地物的混合物，可能也不理想。因此，他设想用没有水体像元的低辐射地物的统计值来校正一般地物。根据研究区的具体特点，提出了用山地阴影直方图主体法确定黑暗像元值。Ｈｏｌｂｅｎ等［８］提出可以利用中红外波段来确定可见光波段的黑暗像元。Ｋａｕｆｍａｎ等［９］发现，可见光波段黑暗像元的反射率与这些像元在中红外波段反射率存在很好的线性关系，根据这种关系，利用中红外波段发现了可见光波段的黑暗像元，这种方法比较客观，精度也很好。２．２　程辐射Ｌｐ的计算根据上述方法可以确定遥感图像中的黑暗像元值，设为ＤＮｍｉｎ。利用传感器定标系数Ｇａｉｎ和Ｂｉａｓ可以计算出程辐射Ｌｐ，即Ｌｐ＝Ｇａｉｎ·ＤＮｍｉｎ＋Ｂｉａｓ（４）　　但要注意的是，虽然这种黑暗像元的反射率或辐射亮度很低，但它并不等于零［６］，如果把上述方法确定的黑暗像元值完全当作是大气的影响，利用它计算程辐射，将产生很大误差。因此，利用上面方法得到的黑暗像元值计算出的辐射亮度，应该为该光谱段向上程辐射亮度（即程辐射）和黑暗像元反射辐射亮度之和，而我们要计算的程辐射应该为黑暗像元值计算出的辐射亮度减去黑暗像元反射辐射亮度的贡献。考虑到上述因素，一些研究者根据具体情况，对黑暗像元反射率进行了一定的假设。如Ｋａｕｆｍａｎ等［１０］假设黑暗像元反射率为０．０２±０．０１，还有一些研究者［６，１１，１２］，假设黑暗像元反射率为０．０１，则程辐射Ｌｐ可表示为Ｌｐ＝Ｇａｉｎ·ＤＮｍｉｎ＋Ｂｉａｓ－（０．０２±０．０１）（ＴθＥ０ｃｏｓθ＋ＥＤ）Ｔ／π（５）３　黑暗像元法中的几种不同大气校正模型确定程辐射Ｌｐ后，再对大气校正模型（３）中的Ｔ、Ｔθ和ＥＤ进行不同的简化假设和不同的计算，这样得到了下面不同的大气校正模型。３．１　第１种大气校正模型这是一种最简单的大气校正模型，该模型仅考虑到了θ、Ｌｐ的影响，没有考虑Ｔθ、Ｔ和ＥＤ的影响［１３］。大气校正模型（３）变为Ｒ＝π（Ｌ－Ｌｐ）／（Ｅ０ｃｏｓθ）（６）　　在可见光波段，大气散射影响大，该模型考虑了程辐射的影响，因此，模型对可见光波段大气校正的效果较好，且波长越短的波段校正效果越好。随着波长的增加，程辐射的影响逐渐减弱，大气吸收影响逐渐增强，而该模型没有考虑Ｔθ、Ｔ的影响，使得该模型的校正精度下降。另外，该模型也没有考虑到的影响，因此，该模型适合于传感器天顶角较小的遥感图像大气校正。３．２　第２种大气校正模型该模型考虑了θ、Ｌｐ及Ｔθ的影响。Ｃｈａｖｅｚ［１３］提出了两种估算Ｔθ的方法：①直接利用已有且较成熟的辐射传输模型中对应波段的大气透射率的默认值来替代Ｔθ；②指出在大多数图像中，大气光学厚度在０．０８～０．３之间，太阳高度角在３０°～５５°之间，此时，Ｔθ可用太阳天顶角的余弦值表示，即Ｔθ＝ｃｏｓθ。则大气校正模型（３）变为Ｒ＝π（Ｌ－Ｌｐ）／（ｃｏｓθＥ０ｃｏｓθ）（７）　　方法①中利用辐射传输模型得到的Ｔθ较方法②中直接用Ｔθ＝ｃｏｓθ更接近实际值，因此精度高。但该模型没有考虑到Ｔ的影响，因此，也只适合于传感器天顶角较小的遥感图像大气校正。与模型１比较，该模型考虑了吸收的影响，因此对可见光波段中的长波段及近红外波段的大气校正精度较高。３．３　第３种大气校正模型该模型考虑了θ、、Ｌｐ、Ｔθ、Ｔ及ＥＤ的影响。在可见光到短波红外波段中，当θ、小于７０°时，遥感图像主要受到大气散射和弱吸收影响，Ｔθ、Ｔ可利用大气光学厚度值τ，采用如下公式计算Ｔθ＝ｅｘｐ（－τｓｅｃθ）（８）Ｔ＝ｅｘｐ（－τｓｅｃ）（９）　　对于大气光学厚度值τ的确定也有不同的方法。有些研究者直接采用已有的美国标准大气光学·９·国　土　资　源　遥　感２００５年厚度值［１４］计算Ｔθ、Ｔ，但是，τ值实际是有时间和地域的差异，直接采用标准的大气光学厚度将存在一定的误差，因此，在实际应用中要对标准的τ值进行修正。还有些研究者利用不同的模型计算τ。例如Ｋａｕｆｍａｎ［１０］在仅考虑瑞利散射因素，没有考虑气溶胶因素的情况下，用式（１０）估算光学厚度τｒτｒ＝０．００８５６９λ－４（１＋０．０１１３λ２＋０．０００１３λ－４）（１０）　　Ｍｏｒａｎ［１２］在假设天空辐照度各向同性的情况下，用４πＬｐ估算了大气层外相应波长的太阳光谱辐照度的损失量，在考虑瑞利散射因素，气溶胶因素的情况下，用式（１１）、（１２）式推算大气光学厚度Ｔθ＝ｅｘｐ（－τｓｅｃθ）＝１－４πＬｐ／Ｅ０ｃｏｓθ（１１）通过上式解大气光学厚度τ，得如下公式τ＝－ｃｏｓθｌｎ（１－４πＬｐ／Ｅ０ｃｏｓθ）（１２）　　利用公式（１１）、（１２）迭代运算，得到稳定的τ后，再通过τ计算Ｔθ、Ｔ。Ｇｉｌａｂｅｒｔ［１］、王庆安等［１５］在推算大气光学厚度τ时，不但考虑了瑞利散射τｒ、气溶胶τａ的影响，还考虑了臭氧光学厚度τｏ的作用，即τ＝τｒ＋τａ＋τｏ（１３）　　τｒ的计算采用了与（１０）式相似的方法。臭氧吸收光学厚度τｏ可认为相对稳定，通过相关文献提供的数据内插得到［１６］。τａ的推算也是个关键。Ｇｉｌａｂｅｒｔ等［１］利用黑暗像元法首先推算出程辐射Ｌｐ，并估算出瑞利散射辐亮度Ｌｒ，再利用Ｌｐ、Ｌｒ计算出气溶胶辐射亮度Ｌａ，最后，根据Ｌａ与气溶胶光学厚度的关系推算出τａ。如果模型不考虑由天空光漫射到地表面的光谱辐照度ＥＤ，则算出的地物反射率偏大。ＥＤ也有不同的确定方法。例如，Ｍｏｒａｎ［１２］用πＬｐ估算ＥＤ，Ｖｅｒｍｏｔｅ等［１７］利用６Ｓ模型估算了ＥＤ，Ｇｉｌａｂｅｒｔ［１］利用算出的τｒ、τａ、τｏ得到ＥＤ，用上面的方法得到Ｔθ、Ｔ、ＥＤ，代入公式（３）得到大气校正模型Ｒ＝π（Ｌ－Ｌｐ）／｛ｅｘｐ（－τｓｅｃ）［ｅｘｐ（－τｓｅｃθ）Ｅ０ｃｏｓθ＋ＥＤ］｝（１４）３．４　第４种大气校正模型该模型也考虑了θ、、Ｌｐ、Ｔθ、Ｔ及ＥＤ的影响，但Ｔ、Ｔθ和ＥＤ都是通过较成熟的大气辐射传输模型得到的［１８］。具体过程为，首先，利用图像中的黑暗像元计算出Ｌｐ。这里要注意，因为可见光波段中的较短波段（例如蓝波段）受程辐射影响大，因此一般选用较短的波段计算Ｌｐ。Ｔ、Ｔθ和ＥＤ可根据计算的Ｌｐ，利用大气辐射传输模型（６Ｓ、ＬＯＷＴＲＡＮ、ＭＯＲＴＲＡＮ）内设的模拟大气状态模型，得到模拟的Ｌｐ，再把Ｌｐ与利用黑暗像元算出的Ｌｐ比较，以最接近Ｌｐ为准，确定最合适的大气模式，这样推导出了Ｔ、Ｔθ和ＥＤ，代入大气校正模型（３）即可。该模型充分利用了已有的大气辐射传输模型得到Ｔ、Ｔθ和ＥＤ，计算速度快、精度高。４　应用实例为直观说明大气校正后的效果，选择模型４（利用６Ｓ模型）对南京地区一景１１月份的ＳＰＯＴ图像（红波段）进行大气校正。从图１可以看出，校正后的图像较原图像更加清晰。图１　大气校正前后图像对比（左：原始图像；右：校正后图像）为了说明不同大气校正模型对大气校正结果的影响，选择模型２（Ｔθ的计算采用方法１）和模型４对图像进行大气校正，并在图像中选择３种典型地物（湖水、草地、沥青路面），与利用光谱仪准同步测量的反射率结果及根据图像地理位置和成像时间查阅资料［１９，２０］得到的３种地物的参考反射率进行比较，结果如表１所示。表１　大气校正前后典型地物反射率（红波段）比较典型地物校正前反射率模型２校正后反射率模型４校正后反射率参考反射率湖水０．０７５０．０２５０．０５１０．０４５草地０．１０６０．０８２０．０９３０．０９０沥青路面０．１２９０．１４３０．１７６０．１８０由表１可见，模型２和模型４都对图像起到了大气校正的作用，但模型４校正效果较好，校正后的值和参考值的误差不超过１５％。这主要因为模型４不但考虑了Ｌｐ、Ｔθ的影响，还考虑了Ｔ、ＥＤ的影响，而且通过６Ｓ模型得到的Ｔ、Ｔθ和ＥＤ精度很高，可以满足目前大部分定量化遥感的应用要求。·０１·第１