反演离水辐亮度流程

离水辐亮度分成三部分进行求解，一是减去大气瑞利散射，然后减去气溶胶的影响，最后除以漫透射率得到离水辐亮度。流程图如下：

目前，我所做的工作主要有如下几块：

1.1遥感影像以及各类参数数据的搜集

前期的数据准备需要大量的光学遥感影像，所以选用了高时间分辨率的MODIS影像作为实验数据。我通过NASA网站，获取到合适时间尺度内（152天到175天）的MODIS光学影像数据。

反演过程中还需要用到气溶胶光学厚度AOD以及臭氧光学厚度，对于气溶胶光学厚度，我们可以用Modis的04级气溶胶产品代替。而臭氧光学厚度需要我们获取实验区域内的臭氧单位。

1.2数据预处理

数据预处理最主要的工作就是对影像几何校正，然后把各类数据重采样为同一分辨率，方便接下来的波段运算工作

1.3波段运算，得出反演结果

波段运算是将离水辐亮度的反演公式从数学表达式转换成idl可读表达式的过程，这也是反演过程中最主要的工作量之一。
波段运算我将其拆成三小部分完成。分别是反演大气瑞利散射辐亮度，气溶胶辐亮度，臭氧光学厚度，以及通过大气漫透过率反演离水辐亮度。
考虑到太阳耀斑相对于瑞丽散射以及气溶胶对离水辐亮度的影响较小，暂时没有进行太阳耀斑的提取工作。

1.3.1 反演大气瑞利散射辐亮度。

水色遥感中，瑞利散射占传感器接收的总辐射的较大份额，特别时在蓝光波段占80%-90%。太阳辐射经过大气层后，有一部分经大气分子后向散射到达卫星传感器，还有一部分到达海面，经海面漫反射后到达传感器，卫星接收到的大气瑞利散射为：

是太阳以三种不同的通道散射到卫星的概率因子， 是单次散射反照率（非水汽和氧气吸收带的波段近似为1），瑞利光学厚度 ，Fs是经太阳-地球距离校正后的辐照度。

IDL波段运算表达式：
Fs*(0.008569*(λ^(-4))(1+0.0113(λ^(-2))+0.00013*λ^(-4))(P/1013.25))((0.75*(1+(-cos (((B2)))cos((B3))-sin((B2))*sin((B3))*cos((B4)-(B5)))^2))+((0.5((sin((B2)-asin(sin ((B2)/1.34)))/sin ((B2)+asin (sin ((B2)/1.34))))^2+(tan ((B2)-asin (sin ((B2)/1.34)))/tan ((B2)+asin (sin ((B2)/1.34))))^2))+(0.5*((sin ((B3)-asin (sin ((B3)/1.34)))/sin ((B3)+asin (sin ((B3)/1.34))))^2+(tan ((B3)-asin (sin ((B3)/1.34)))/tan ((B3)+asin (sin ((B3)/1.34))))^2)))(0.75(1+(cos ((B2)) cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5)))^2)))(1/(cos ((B2))+0.15*(93.885-(B2))^(-1.253) *(P/1013.25))) /(4*3.1415926)

Figure 1瑞利散射辐亮度

1.3.2. 反演气溶胶辐射亮度值。

Gordon等人认为气溶胶的水平分布通常比较均匀，假设气溶胶类型、气溶胶光学厚度在小范围内(（50～100 km）基本不变 。气溶胶散射辐射亮度值为:

其中， Wa是气溶胶单次散射反照率，一般取值为0.9，Pa是经气溶胶散射到卫星的概率因子0.368863， 是气溶胶光学厚度为0.896， Fs是经太阳-地球距离校正后的辐照度。求取结果如下所示：
Fs*0.9*0.896*(((0.983*(1-0.82*0.82))/(1+0.82*0.82-0.82*0.82*(cos ((B2)) cos ((B3))-sin((B2))*sin((B3))*cos((B4)-(B5))))^1.5+((1-0.983)(1-(-0.55)^2))/(1+(-0.55)^2-2*(0.55)(cos ((B2)) *cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5))))^1.5)+(0.5((sin ((B2)-asin (sin ((B2)/1.34)))/sin ((B2)+asin (sin ((B2)/1.34))))^2+(tan ((B2)-asin (sin ((B2)/1.34)))/tan ((B2)+asin (sin ((B2)/1.34))))^2)+0.5*((sin ((B3)-asin (sin ((B3)/1.34)))/sin ((B3)+asin (sin ((B3)/1.34))))^2+(tan ((B3)-asin (sin ((B3)/1.34)))/tan ((B3)+asin (sin ((B3)/1.34))))^2))((0.983(1-0.82*0.82))/(1+0.82*0.82-0.82*0.82*(-cos ((B2)) cos ((B3))-sin ((B2))*sin((B3))*cos((B4)-(B5))))^1.5+((1-0.983)(1-(-0.55)^2))/(1+(-0.55)^2-2*(0.55)(-cos ((B2)) *cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5))))^1.5)) (1/(cos ((B2))+0.15*(93.885-(B2))^(-1.253) *(P/1013.25))) /(4*3.1415926)

Figure 2气溶胶辐射亮度

1.3.3. 反演臭氧的光学厚度 和大气漫透射率t(λ,θ)

臭氧的光学厚度 求法如下：

其中： 为臭氧吸收系数（cm-1），为大气臭氧含量， 以cm为单位时，一般在0.25-0.35之间，或   ，（DU为多普逊臭氧单位）

波段 波段宽度/nm 中心波长/nm F0 臭氧吸收系数 空间分辨率
1 620-670 645 163.00 0.07 250m
2 841 858 100.58 0 250m
8 405~420 412 180.50 0 1000m
9 438~448 443 194.50 0.003 1000m
10 483~498 488 187.50 0.019 1000m
11 526~536 531 196.00 0.064 1000m
12 546~556 551 186.50 0.085 1000m
13 662~672 667 154.00 0.049 1000m
14 673~683 678 150.50 0.04 1000m
15 743~753 748 128.05 0.0092 1000m
16 862~877 865 99.05 0 1000m
DU=346 =0.346，吸收系数如上表，臭氧光学厚度 =0.02422。

大气漫透射率t(λ,θ)的求法如下：
大气漫透射率与大气中各粒子的单次和多次散射有关，计算在任意波长上的大气漫透射率的目的是为了获得大气校正方程中的离水辐亮度，属于大气校正工作中必不可少的一环。t(λ,θ) 的化简，在漫透射光的条件下，fRωR≈0.5，fAωA≈0，t(λ,θ)可近似地简化为：

IDL波段运算表达式：

(exp(-((0.008569*λ^(-4)(1+0.0113λ^(-2)+0.00013*λ^(-4))*(P/1013.25))/2+0.896+0.02422)/cos (B2)))

Figure 3大气漫透射率

1.3.4. 反演离水辐亮度

水色遥感大气校正方程：

其中：L(λ)是星探测到的辐亮度；LR(λ) 是大气分子瑞利散射的辐亮度；LA(λ) 是气溶胶散射的辐亮度，Lw(λ)是离水辐亮度，t是大气漫透射率。对上述公式进行转换，获取离水辐亮度Lw(λ)，求法如下：

IDL波段运算表达式：
(B1-(Fs*(0.008569*(λ^(-4))(1+0.0113(λ^(-2))+0.00013*λ^(-4))(P/1013.25))((0.75*(1+(-cos (((B2)))cos((B3))-sin((B2))*sin((B3))*cos((B4)-(B5)))^2))+((0.5((sin((B2)-asin(sin ((B2)/1.34)))/sin ((B2)+asin (sin ((B2)/1.34))))^2+(tan ((B2)-asin (sin ((B2)/1.34)))/tan ((B2)+asin (sin ((B2)/1.34))))^2))+(0.5*((sin ((B3)-asin (sin ((B3)/1.34)))/sin ((B3)+asin (sin ((B3)/1.34))))^2+(tan ((B3)-asin (sin ((B3)/1.34)))/tan ((B3)+asin (sin ((B3)/1.34))))^2)))(0.75(1+(cos ((B2)) cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5)))^2)))(1/(cos ((B2))+0.15*(93.885-(B2))^(-1.253) (P/1013.25))) /(4*3.1415926))-( Fs*0.9*0.896(((0.983*(1-0.82*0.82))/(1+0.82*0.82-0.82*0.82*(cos ((B2)) cos ((B3))-sin((B2))*sin((B3))*cos((B4)-(B5))))^1.5+((1-0.983)(1-(-0.55)^2))/(1+(-0.55)^2-2*(0.55)(cos ((B2)) *cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5))))^1.5)+(0.5((sin ((B2)-asin (sin ((B2)/1.34)))/sin ((B2)+asin (sin ((B2)/1.34))))^2+(tan ((B2)-asin (sin ((B2)/1.34)))/tan ((B2)+asin (sin ((B2)/1.34))))^2)+0.5*((sin ((B3)-asin (sin ((B3)/1.34)))/sin ((B3)+asin (sin ((B3)/1.34))))^2+(tan ((B3)-asin (sin ((B3)/1.34)))/tan ((B3)+asin (sin ((B3)/1.34))))^2))((0.983(1-0.82*0.82))/(1+0.82*0.82-0.82*0.82*(-cos ((B2)) cos ((B3))-sin ((B2))*sin((B3))*cos((B4)-(B5))))^1.5+((1-0.983)(1-(-0.55)^2))/(1+(-0.55)^2-2*(0.55)(-cos ((B2)) *cos ((B3))-sin ((B2))*sin ((B3))*cos((B4)-(B5))))^1.5)) (1/(cos ((B2))+0.15*(93.885-(B2))^(-1.253) (P/1013.25))) /(4*3.1415926)))/ (exp(-((0.008569λ^(-4)(1+0.0113λ^(-2)+0.00013*λ^(-4))*(P/1013.25))/2+0.896+0.02422)/cos (B2)))

在实际的反演过程中遇到了很多难题，总的来说，有数据获取上的难题，反演过程上的难题，还有实际操作过程中的难题

2.1 数据获取中存在的问题

2.1.1 臭氧数据

ppt里提到的臭氧的一些资料现在已经无法使用，不能获得准确的臭氧单位。

所以第一个难题就是获取臭氧单位。后来经过多次查找，在NASA的GSFC中找到臭氧的获取网址。
臭氧数据下载：：https://ozoneaq.gsfc.nasa.gov/data/ozone/
NASA以前OMI数据已经不再提供，被新的臭氧传感器OMPS所代替。但是历史数据仍然能被检索到。

Figure 4臭氧数据下载

我一开始想得到臭氧的分布文件，但是下载后不知道如何打开，后来直接在GSFC的Tools里面找到臭氧分布的地图，可以得到全球范围内任意地区的臭氧单位DU。
臭氧map查询（建议）：可以直接得到多普逊臭氧单位
臭氧map网址：https://ozoneaq.gsfc.nasa.gov/tools/ozonemap/

2.1.2 气溶胶数据产品。

暂时用MODIS气溶胶04级的气溶胶光学厚度产品。但是对于MODIS04级产品不熟悉，没有了解到04级各波段所代表的产品是什么。
考虑到水汽对气溶胶，对大气可能存在的影响，接下来的任务之一便是如何消除水汽的影响（垂直订正）。

2.2 波段运算中存在的问题

最开始几次反演离水辐亮度得到的图都是无效值NaN，和老师探讨多次，终于找到一些解决的办法。

(1)MODIS中的角度都是以特定格式存储的，在波段运算时要乘比例系数scale=0.01，然后转换成弧度制。一开始对波段运算了解的不多，计算用的是角度制，得到的概率因子都是无效值NaN。

Figure 5 MODIS太阳高度角的存储格式

(2)波段运算时，系数尽量以小数形式表示。
举例：比如3/4尽量以0.75表示。用分子分母形式表示容易出错。
(3)反演离水辐亮度波段运算公式较长，不能一次性将所有的公式全部输入。尽量分步走，将公式拆分成几个部分分块处理。这样即可以得到各部分的结果图，又便于检查波段运算过程中的错误。

3离水辐射亮度反演结果分析
将得到的辐亮度文件打开，以rainbow色表显示。
可以在十字光标附近看到明显的白色轮廓区域，这片区域疑似油膜（１６５天）的油膜位置，它的地理范围是在山东长岛以西的渤海海域。
对应的新闻报道中的油膜覆盖区域：
5
Figure 6油膜影响区域
实际反演得到的离水辐亮度影像：

Figure 7实际反演的离水辐亮度图

Figure 8 离水辐亮度图中疑似油膜放大图
如第一张图所示，油膜实际影响的范围主要处于长岛西北附近海域，渤海的中部。而通过我们反演得到的离水辐亮度图来看，在长岛以西附近海域有一明显区域，其离水辐射亮度高于附近位置，大致轮廓如Figure 8所示。通过放大图，可以清楚的看到油膜的轮廓形状。

然后我做了一项工作就是针对疑似油膜区域以及周边海水的离水辐亮度信息做了一个对比。

油膜区域的离水辐亮度值

Figure 10疑似油膜b点 离水辐亮度172.77
海水的离水辐亮度值（渤海内，渤海口，太平洋三处）

Figure 11渤海内海水 离水辐亮度131.75

Figure 12渤海口海水 离水辐亮度95.31

Figure 13太平洋海水 离水辐亮度55.21

对该区域是否为油膜的判断依据：
１通过图像目视判读，
疑似油膜与周围云在形状上有明显区别，而且石油平台的地理位置在疑似油膜的附近。所以判定该区域是油膜。

２通过离水辐亮度值。
疑似油膜区域的离水辐亮度值普遍在１４８以上，最高值达到了２６７，
而海水的离水辐亮度值在近海区域以及远海区是不同的，近海区域在１２５－１３３之间，而远海区是在８０－１００左右。深海区一类水体（太平洋上）的离水辐亮度３０－６０左右。

工作展望：

1.进一步完善细节.气溶胶部分的反演过程较为粗糙，只是简单的利用MODIS气溶胶产品，没有考虑水汽对大气的影响，下一步要展开更深层次的气溶胶散射机制研究。

2.考虑利用交叉验证的方法来验证该反演方法的准确程度。
利用landsat该时间段内渤海区域的光学影像，反演出离水辐射亮度。与MODIS反演得到的离水辐射亮度进行交叉验证。