专利名称：非均匀大气程辐射遥感数字图像的计算机生成方法
技术领域：
本发明涉及一种遥感数字图像的获取方法，具体地说，是涉及一种非均匀大气程辐射遥感数字图像的计算机生成方法。
背景技术：
大气程辐射是太阳辐射在大气传输过程中经大气中各组分以及气溶胶微粒散射后直接到达遥感器的辐射。由于大气成分在空间分布的不均匀和大气的污染，随着空气组分和气溶胶含量、形状、性质、粒径、浓度等因素的空间分布差异，卫星遥感数字图像上各点的大气程辐射值是不同的。这些不同的大气程辐射值是太阳辐射在大气传输过程中与不同地点的不同空气组分和气溶胶相互作用的结果，表征着该点地面上空的大气状况和质量。
由于大气程辐射是与地面状况无关，且表征大气质量状况的天空光亮度遥感分量，所以对于作为对地探测技术的遥感而言，大气程辐射值是叠加在地面遥感值上与地面信息无关的大气干扰。因此，大气程辐射遥感数字图像可作为是研究区域大气环境、反演和监测城市大气环境质量的理想的遥感信息源，同时也是研究太阳辐射的大气传输机制、机理和进行遥感数字图像大气修正的重要遥感资料，具有广泛的应用价值。例如，中国发明专利申请公开了“一种大气程辐射遥感数字图像的计算机生成方法”(专利申请号为200510027276.X)，该方法采用遥感数字图像中的大气程辐射值构建传统遥感手段无法获得的大气遥感数字图像。但是，该方法在生成大气程辐射遥感数字图像计算过程中，没有具体考虑非标准朗伯体和标准朗伯体之间的物理机制区别，即对所有图斑内部像元采用相同公式完成内插，产生了在非标准朗伯体上大气程辐射值的计算误差，降低了大气程辐射值的计算精度。再者，对于传统的大气程辐射计算模型，如回归分析法、直方图法，都在不同程度上将各个像素点上的大气状况假定为相同情况，且某些大气参数无法精确的获得，造成大气程辐射值计算精度低，效果差。
发明内容
鉴于以上所述现有技术存在的问题和不足，本发明的目的是提供一种非均匀大气程辐射遥感数字图像的计算机生成方法，该方法能提高大气程辐射值计算精度，实现准确地生成大气程辐射遥感数字图像。
为达到上述目的，本发明采用的技术解决方案如下 上述一种非均匀大气程辐射遥感数字图像的计算机生成方法，其步骤包括 A、在IDL/ENVI系统支持下，读取选定地区的MODIS02和MODIS09的遥感数字图像，参阅图4； B、根据MODIS09地面反射率数据和MODIS02表观反射率数据，分别计算水体掩模图像和云掩模图像，进而可得陆地图像，即陆地图像是从MODIS09图像(选定地区的整幅图像)中筛除水体掩模图像和云掩模图像后的图像； C、计算水体图像、陆地图像和云层图像边界的大气程辐射值； D、分别逐点计算水体图像、陆地图像和云层图像内部大气程辐射值，完成所有像元大气程辐射值的计算； E、完成大气程辐射图像的可视化，从而完成大气程辐射遥感数字图像的计算机生成。
本发明的非均匀大气程辐射遥感数字图像的计算机生成方法优点在于该方法根据MODIS09地面反射率数据和MODIS02表观反射率数据，分别计算水体掩模图像和云掩模图像，采用逐点像素计算的方法，实现对非均匀大气状况的大气程辐射遥感数字图像的生成；根据MODIS09地面反射率数据和MODIS02辐射亮度数据，采用反演计算方法，可以在不需要大气参数的条件下，精确反演大气程辐射遥感数字图像，适用性广。
图1为本发明的一个实施例(即中国长江三角洲地区)的MODIS09遥感数字图像； 图2为与图1相匹配的计算机生成的MODIS大气程辐射遥感数字图像； 图3为本发明的非均匀大气程辐射遥感数字图像的计算机生成方法的流程图； 图4为图3中步骤B中步骤B1计算水体掩模图像的流程图； 图5为图3中步骤D中步骤D1计算陆地图像内部大气程辐射值的流程图； 图6为图3中步骤D中步骤D2和D3计算水体和云层内部大气程辐射值的流程图。
具体实施例方式 下面根据图1至6给出本发明一种非均匀大气程辐射遥感数字图像的计算机生成方法的实施例，要指出的是，所给出的实施例是为了说明本发明方法的技术特点和功能特点，能更易于理解本发明，而不是用来限制本发明的范围。
本发明的非均匀大气程辐射遥感数字图像的计算机生成方法的实施例是在IDL/ENVI系统支持下进行的，先请参阅图3、图4、图5和图6，具体步骤如下 A、在IDL/ENVI系统支持下，读取选定地区的MODIS02和MODIS09的遥感数字图像，如图1所示； B、根据MODIS09地面反射率数据和MODIS02表观反射率数据，分别计算水体掩模图像和云掩模图像，进而可得陆地图像，即陆地图像是从MODIS09图像(选定地区的整幅图像)中筛除水体掩模图像和云掩模图像后的图像， 遥感数字图像上任一像元点的大气程辐射值计算是计算机生成大气程辐射遥感数字图像的关键，将遥感数字图像上像元点分为同一地物的内部像元点和边界像元点。以MODIS图像为例，将其分为水体图像、陆地图像和云层图像，即分别建立水体掩模和云掩模图像，具体步骤如下 B1、根据水体各波段反射率的性质，利用NDVI(归一化植被指数)计算水体掩模图像，其具体如下 首先，计算NDVI遥感数字图像，其计算表达式为 同时生成其直方图，再利用数学求导方法计算出NDVI直方图上第一个陡升点坐标，记为阈值K1； 其次，根据水体的物理特性，即水体在近红外波段上具有低反射率以及在蓝光波段上具有高反射率的特性，利用MODIS09的蓝光波段(BAND3)和近红外波段(BAND2)数据，分别计算蓝光波段(BAND3)和近红外波段(BAND2)直方图，再利用数学求导方法分别计算出蓝光波段(BAND3)和近红外波段(BAND2)直方图的陡降和陡升拐点坐标值，蓝光波段(BAND3)计算陡降拐点坐标值，记为阈值K2，近红外波段(BAND2)计算陡升拐点坐标值，记为阈值K3，即确定水体掩模其他两个提取阈值； 最后，利用布尔计算， RIVER_MASK＝(NDVI＜K1)AND(BAND3＞K2)AND(BAND2＜K3)，得到水体掩模图像，参阅图4。
B2、根据MODIS02表观反射率数据，计算云掩模图像。
如果表观反射率值BAND1_DN＞0.14，则此像素为云层，此像元值赋值为1；若BAND1_DN＜0.14，则此像素不是云层，赋值为0，得到云掩模图像。
C、基于水体掩模图像和云掩模图像，计算水体图像、陆地图像和云层图像边界的大气程辐射值，其具体步骤如下 设立水体图像、陆地图像和云层图像边界上的大气程辐射值的逐点成像计算表达式 假设在小范围内边界上的大气状况相同，根据遥感数字图像的点成像通用公式可得 DNi＝k×Hi×ri×τi/π+DNAi DNj＝k×Hj×rj×τj/π+DNAj 其中，DN为遥感像元的辐射亮度值，H为卫片像元的照度，r为地面反射率，τ为大气透过率，DNA为大气程辐射值，若在小范围内边界上，大气视为均匀大气，即DNAi＝DNAj，τi＝τj，由上式联立方程可得大气程辐射逐点成像公式 DNA＝DNi-ri·(DNi-DNj)/(ri-rj) D、分别逐点计算陆地图像、水体图像和云层图像内部大气程辐射值，完成所有像元大气程辐射值的计算，具体步骤如下 D1、设置一个用于控制计算精度和运算速度的阈值ε2，计算陆地内部大气程辐射值，如图5所示，其具体如下 设定一个与计算结果的精度以及运算速度相关的阈值ε2，α≤ε2≤β，若阈值ε2小于α，则计算陷入死循环，运算速度变慢；若ε2值大于β，则运算速度提高，但计算结果不是最优，整体效果精度低，参阅图5，其具体算法步骤如下 D11、利用陆地边界上的大气程辐射值，计算A0值，并进行初始化赋值Ka＝1，Aka＝A0，DNAka为图像边界上第k个像元的大气程辐射值，L为陆地边界上总像元数； D12、计算DN′ka＝DNka-Aka，得到MODIS02辐射亮度值与陆地边界大气程辐射值之间的差值； D13、建立回归方程DNAka＝aln DN′ka+b，具体步骤如下 采用线性对数拟合方法，建立大气程辐射值与MODIS02辐射亮度值两者之间的回归方程 DNAk＝a*ln(DNk-DNAk)+b＝a*ln(DN′k)+b 令Y＝DNAk，X＝ln(DNk-DNAk)＝ln(DN′k) 建立回归方程Y＝a*X+b 其中b＝Y′-aX 其中L为陆地边界上总像元数，k为第k个像元的下标，DNk为图像边界上第k个像元的辐射亮度值，DNAk为图像边界上第k个像元的大气程辐射值； D14、设定阈值ε2，若|DNAka-Aka|＞ε2为真，则Ka＝Ka+1，Aka＝Aka±1，返回步骤D12，重新建立回归方程；若|DNAka-Aka|＞ε2为假，则转步骤D15； D15、输出DNAka，即为所求的大气程辐射值，DNAij＝DNAka。
D2、计算水体内部的大气程辐射值，水体为非标准朗伯体，主要为海洋，如图6所示，其具体步骤如下 D21、在计算水体大气程辐射值前，首先需推导遥感数字图像的点成像通用公式 DN＝k*H*r*τ+DNA＝k*E0*r*τ上*τ下/π+DNA＝k*E0*r*τ2/π+DNA 其中，τ上和τ下分别为上行大气透过率和下行大气透过率，若在小范围内，视为τ上＝τ下。
水体在近红外波段具有低反射率，即DN1＝k*H*r*τ几乎等于0，所以，MODIS02近红外波段的辐射亮度值(DN)近似等于近红外波段上的大气程辐射值(DNA)。利用近红外波段和需计算大气程辐射值(DNA)波段之间的线性相关性，按下面步骤联立方程 由方程组(1)可得 由方程组(2)可得 推导出 若假设DNA3＝DN3，进而得到水体上空大气程辐射值计算方程 其中，DN1，DN2为需计算大气程辐射值波段上两个像素的辐射亮度值，DN3，DN4为近红外波段图像上两个像素的辐射亮度值。
D22、为保证MODIS09图像中的(r1-r2)和(r3-r4)差值足够大，设定一个限定反射率差值范围的阈值ε1，γ≤ε1≤η，若阈值ε1小于γ，则计算结果可出现负值，计算速度变慢；若ε1值大于η，则计算结果可出现非法的极大值，同时结果受薄云影响增大，大气程辐射图像生成效果不理想。
若|r1-r2|＞ε1and|r3-r4|＞ε1为真，则转步骤D23，若|r1-r2|＞ε1and|r3-r4|＞ε1为假，则按照均值内插方法计算水体内部大气程辐射值，转步骤D32。
D23、使用公式(3)完成水体图像内部的大气程辐射值计算，从而完成水体的大气辐射值计算。
D3、利用插值算法计算云层或水体图像内部大气程辐射值，其具体步骤如下 D31、计算云层图像内部大气程辐射值。采用均值内插的方法计算云层图像内部大气程辐射值，由于云层为非标准朗伯体，且云层在遥感数字图像上为高亮度的物体，若完全利用步骤D1中陆地内部大气程辐射值计算方法，则计算结果为云层上方大气程辐射值，而利用邻近像元均值内插方法则可计算云层下方的大气程辐射值； D32、建立3*3矩阵(‘a’)，利用均值内插方法完成云层或水体图像内部的大气程辐射值计算； D33、计算云层或水体图像内部大气程辐射值 基于上述步骤B中生成的云掩模图像和水体掩模图像，利用3*3矩阵(‘a’)分别遍历整幅大气程辐射图像和云掩模图像，或分别遍历整幅大气程辐射图像和水体掩模图像，按下面公式计算云层或水体图像内部大气程辐射值 其中，N为3*3矩阵(‘a’)遍历云掩模图像或水体掩模图像时m a s k＝1的个数。
E、当陆地图像、水体图像和云层图像逐点大气程辐射值计算完成后，将全部大气程辐射值输入ENVI/IDL遥感处理软件，利用ENVI/IDL软件本身程序完成大气程辐射图像的可视化，从而完成大气程辐射遥感数字图像的计算机生成，实施例即为可视化过程，具体可视化结果见图2。
权利要求
1.一种非均匀大气程辐射遥感数字图像的计算机生成方法，其步骤包括
A、在IDL/ENVI系统支持下，读取选定地区的MODIS02和MODIS09的遥感数字图像；
B、根据MODIS09地面反射率数据和MODIS02表观反射率数据，分别计算水体掩模图像和云掩模图像，进而可得陆地图像；
C、计算水体图像、陆地图像和云层图像边界的大气程辐射值；
D、分别逐点计算水体图像、陆地图像和云层图像内部大气程辐射值，完成所有像元大气程辐射值的计算；
E、完成大气程辐射图像的可视化，从而完成大气程辐射遥感数字图像的计算机生成。
2.根据权利要求1所述的非均匀大气程辐射遥感数字图像的计算机生成方法，其特征在于，
上述步骤B所述的分别计算水体掩模图像和云掩模图像，其具体步骤如下
B1、利用NDVI计算水体掩模图像，其具体如下
首先计算NDVI遥感数字图像，其计算表达式为
同时生成其直方图，再利用数学求导方法计算出NDVI直方图上第一个陡升点坐标，记为阈值K1；
其次，分别计算蓝光波段(BAND3)和近红外波段(BAND2)直方图，再利用数学求导方法分别计算出蓝光波段(BAND3)和近红外波段(BAND2)直方图的陡降和陡升拐点坐标值，蓝光波段(BAND3)计算陡降拐点坐标值，记为阈值K2，近红外波段(BAND2)计算陡升拐点坐标值，记为阈值K3，即确定水体掩模其他两个提取阈值；
最后，利用布尔计算，
RIVER_MASK＝(NDVI＜K1)AND(BAND3＞K2)AND(BAND2＜K3)，得到水体掩模图像。
B2、根据MODIS02表观反射率数据，计算云掩模图像。
如果表观反射率值BAND1_DN＞0.14，则此像素为云层，此像元值赋值为1；若BAND1_DN＜0.14，则此像素不是云层，赋值为0，得到云掩模图像。
3.根据权利要求2所述的非均匀大气程辐射遥感数字图像的计算机生成方法，其特征在于，上述步骤C所述的计算水体图像、陆地图像和云层图像边界的大气程辐射值，其具体步骤如下
设立水体图像、陆地图像和云层图像边界上的大气程辐射值的逐点成像计算表达式
假设在小范围内边界上的大气状况相同，根据遥感数字图像的点成像通用公式可得
DNi＝k×Hi×ri×τi/π+DNAi
DNj＝k×Hj×rj×τj/π+DNAj
其中，DN为遥感像元的辐射亮度值，H为卫片像元的照度，r为地面反射率，τ为大气透过率，DNA为大气程辐射值，若在小范围内边界上，大气视为均匀大气，即DNAi＝DNAj，τi＝τj，由上式联立方程可得大气程辐射逐点成像公式
DNA＝DNi-ri·(DNi-DNj)/(ri-rj)。
4.根据权利要求3所述的非均匀大气程辐射遥感数字图像的计算机生成方法，其特征在于，上述步骤D所述的分别逐点计算水体图像、陆地图像和云层图像内部大气程辐射值，完成所有像元大气程辐射值的计算，具体步骤如下
D1、设置一个用于控制计算精度和运算速度的阈值ε2，计算陆地内部大气程辐射值，其具体如下设定一个与计算结果的精度以及运算速度相关的阈值ε2，α≤ε2≤β，若阈值ε2小于α，则计算陷入死循环，运算速度变慢；若ε2值大于β，则运算速度提高，但计算结果不是最优，整体效果精度低，其具体算法步骤如下
D11、利用陆地边界上的大气程辐射值，计算A0值，
并进行初始化
赋值Ka＝1，Aka＝A0，DNAka为图像边界上第k个像元的大气程辐射值，L为陆地边界上总像元数。
D12、计算DN′ka＝DNka-Aka，得到MODIS02辐射亮度值与陆地边界大气程辐射值之间的差值。
D13、建立回归方程DNAka＝a ln DN′ka+b，具体步骤如下
采用线性对数拟合方法，建立大气程辐射值与MODIS02辐射亮度值两者之间的回归方程
DNAk＝a*ln(DNk-DNAk)+b＝a*ln(DN′k)+b
令Y＝DNAk，X＝ln(DNk-DNAk)＝ln(DN′k)
建立回归方程Y＝a*X+b
其中
其中
L为陆地边界上总像元数，k为第k个像元的下标，DNk为图像边界上第k个像元的辐射亮度值，DNAk为图像边界上第k个像元的大气程辐射值。
D14、设定阈值ε2，若|DNAka-Aka|＞ε2为真，则Ka＝Ka+1，Aka＝Aka±1，返回步骤D12，重新建立回归方程；若|DNAka-Aka|＞ε2为假，则转步骤D15；
D15、输出DNAka，即为所求的大气程辐射值，DNAij＝DNAka。
D2、计算水体内部的大气程辐射值，水体为非标准朗伯体，其具体步骤如下
D21、在计算水体大气程辐射值前，首先需推导遥感数字图像的点成像通用公式DN＝k*H*r*τ+DNA＝k*E0*r*τ上*τ下/π+DNA＝k*E0*r*τ2/π+DNA其中，τ上和τ下分别为上行大气透过率和下行大气透过率，若在小范围内，视为τ上＝τ下。利用近红外波段和需计算大气程辐射值(DNA)波段之间的线性相关性，按下面步骤联立方程
由方程组(1)可得
由方程组(2)可得
推导出
若假设DNA3＝DN3，进而得到水体上空大气程辐射值计算方程
其中，DN1，DN2为需计算大气程辐射值波段上两个像素的辐射亮度值，DN3，DN4为近红外波段图像上两个像素的辐射亮度值。
D22、设定一个限定反射率差值范围的阈值ε1，γ≤ε1≤η，若阈值ε1小于γ，则计算结果可出现负值，计算速度变慢；若ε1值大于η，则计算结果可出现非法的极大值，同时结果受薄云影响增大，大气程辐射图像生成效果不理想。若|r1-r2|＞ε1and|r3-r4|＞ε1为真，则转步骤D23，若|r1-r2|＞ε1and|r3-r4|＞ε1为假，则按照均值内插方法计算水体内部大气程辐射值，转步骤D32。
D23、使用公式(3)完成水体图像内部的大气程辐射值计算，从而完成水体的大气辐射值计算。
D3、利用插值算法计算云层或水体图像内部大气程辐射值，其具体步骤如下
D31、计算云层图像内部大气程辐射值，采用均值内插的方法计算云层图像内部大气程辐射值，利用邻近像元均值内插方法计算云层下方的大气程辐射值；
D32、建立3*3矩阵(‘a’)，利用均值内插方法完成云层或水体图像内部的大气程辐射值计算；
D33、计算云层或水体图像内部大气程辐射值基于上述步骤B中生成的云掩模图像和水体掩模图像，利用3*3矩阵(‘a’)分别遍历整幅大气程辐射图像和云掩模图像，或分别遍历整幅大气程辐射图像和水体掩模图像，按下面公式计算云层或水体图像内部大气程辐射值
其中，N为3*3矩阵(‘a’)遍历云掩模图像或水体掩模图像时mask＝1的个数。
全文摘要
本发明公开了一种基于颜色分割的多目摄像立体匹配获取深度图像的方法，它包括步骤(1)对所有输入图像进行图像规正；(2)对参考图像进行颜色分割，提取图像中的颜色一致性区域；(3)对多幅输入图像分别进行局部窗口匹配得到多幅视差图；(4)应用双向匹配策略来消除在匹配过程中产生的误匹配点；(5)将多幅视差图合成为一幅视差图，填充误匹配点的视差信息；(6)对视差图进行后处理优化，得到一个稠密视差图；(7)根据视差与深度之间的关系，将视差图转换为深度图。该方法从多个视点图像中获取深度信息，利用多个视点图像提供的图像信息，不仅能解决图像中周期性重复纹理特征、遮挡等带来的误匹配，还能提高匹配精度，得到一个准确的深度图像。
文档编号G01W1/00GK101713832SQ20091019954
公开日2010年5月26日 申请日期2009年11月26日 优先权日2009年11月26日
发明者时向勇, 李先华, 朱瑞芳, 朱姗 申请人:上海大学