CN101566692A - 利用卫星遥感数据中的云影信息检测云高的方法 - Google Patents

Abstract

本发明涉及一种利用卫星遥感数据中的云影信息检测云高的方法，属于大气探测领域。包括插值计算、云和云阴影检测、MODIS数据的云和云影匹配、建立云高模型求得云层高度和云高计算和分析五部分。本方法对遥感图像数据中云和云影的信息进行检测；结合数字图像技术进行云和云影匹配，得到云和云影的之间的距离；最后，建立云高模型求得云层高度。本方法可以同一时刻对一片云多个点进行高度计算，以此求得的高度值将极大的提高云高的检测精度。通过本方法可以对符合测量条件的任何地区进行云高测量，而且不用在当地设置测量仪器。

Description

利用卫星遥感数据中的云影信息检测云高的方法

技术领域

本发明涉及大气探测领域，特别涉及一种利用卫星遥感数据中的云影信息检测云高的方法。

背景技术

云在天气系统的发展、降水的形成和大气辐射传输等物理过程中都扮演着极其重要的角色。如何准确、及时地获取云的信息，对于天气气候学研究、气象预报、人工影响天气以及国民经济和军事等诸多领域都十分重要。通常对云的气象观测主要有：云状、云量和云高三项，其中云高是指检测站地平面至云底的垂直距离，通常以米为单位。

目前，对于云量和云状的观测主要还是靠目测，但云高测量的各项技术业已成熟，出现了许多用于测量云高的仪器。因此，云高测量方法大致可分为以下三类：估测云高、仪器检测云高、利用卫星数据和空间几何反演获得云高。

(一)估测云高：包括目测和经验公式估算。

1.目测：

主要是根据云的形状来估测云高，所以必须正确地判定云的形状，然后再根据云的结构、大小、亮度、颜色、移动速度等情况，结合当地经常出现的云高范围进行估计。

2.经验公式估算：

除了目测云高的方法外，还可以利用经验公式估算云高，这主要针对某些特定云的估测，需要专业技术人员才能估算。

(二)仪器测云高。

1.云幕灯测云高：

云幕灯是一种云高指示器，是最早用来测云高的仪器，主要由光源和测角器两部分组成，通过测得对应的仰角、光源和测高仪之间的距离是来计算云高。

2.旋转光束云高仪测云高：

该系统由一个发射器、一个接收器和一台记录设备组成。通过测试，结合相应公式来计算云高。

3.气球测云高：

主要利用一只充满氢气或氦气的橡胶气球，从地面匀速至云底所经历的时间来测量云高。

4.微脉冲激光雷达(Micropulse Lidar，MPL)测云高：

MPL探测云高，主要是对激光回波数据进行采集，通过数据分析反演云层高度。由于云的种类繁多，形状千姿百态，每块云中云滴的尺度、浓度和含水量不同，具有强散射性，所以会引起后向散射的激光回波也不相同，所以由该仪器测量确定的云底始终是一个有争议的问题。

5.半导体激光云高仪：

半导体激光云高仪是激光雷达的一种，主要用来测量云的高度。半导体激光云高仪是近年来发展起来的，采用固态的近红外激光二极管和数字信号处理技术，独特的光学镜头设计，小型性便携，能使实时在线监测等优点是其它仪器不可比拟的，所以成为目前研究的一个热点领域。

通过对以上云高探测方法的分析比较得出：云幕灯和旋转光束云高仪采用的都是普通的光源，探测距离不高，只能用于低云探测，在降雨时测量不准确，但是，成本低、方法简单、容易操作。探空气球测云成本高，不容易操作，很少使用，但是能够同时探测到其它的一些气象参数。随着激光和电子技术的发展，测云激光雷达(激光测云仪)技术得到很大提高，很多国家都在花费大量的人力、物力进行研究，激光雷达是已成为测量云层高度的有效工具。目前，国内外对于云高的测量，主要采用半导体激光云高仪来测量。国外对云高仪的研制比较早，技术相对比较成熟，已经研发出很多测云的仪器；而国内在这方面起步较晚，主要还是靠进口的仪器，与国外还有一定技术差距。

(三)利用卫星数据信息和空间几何反演获得云层高度。

1.双站数字摄像测量云高方法

2005年，中国科学技术大学地球和空间科学学院与中国气象科学研究院发表双站数字摄像测量云高方法，其基本原理是利用双站测距法。不足之处在于该测试方法受观测场地和光学机械加工精度的限制，测试中需对测量系统(CCD等)进行精密的设计和严格的校准，否则会产生较大测试误差。

2.基于MODIS卫星数据的二氧化碳(CO2)薄片算法

MODIS的全称为中分辨率成像光谱仪，是搭载在Terra和Aqua卫星上的一个重要的传感器，MODIS是第一个拥有CO2薄片波段的高空间分辨率探测器。CO2薄片算法是使用被动遥感资料反演云顶性质的重要方法之一。CO2薄片技术推断云顶气压和有效云量或辐射，是利用15μm CO2吸收带的几个吸收不同的MODIS红外波段，这几个波段对于大气不同高度的层次都有不同的敏感度。CO2薄片技术可以反演出云顶气压和有效云量主要是中高层云，但对于云顶气压大于700hPa(接近地面)的云层，因为云信息的减少，所以不能再用这种方法。因此，CO2薄片技术对低层云的反演存在一些问题，对于云高检测有一定局限性。

3.双星立体观测云顶高度几何反演方法

2008年，北京大学物理学院和中国民用航空总局天津空中交通管理站气象台联合发表双星立体观测云顶高度几何反演方法。该方法存在的问题是：云顶高度反演方法实验中采用的“立体图像对”是卫星立体观测模拟资料。该方法中的立体测量依赖两颗卫星对给定区域的同步观测，然后对卫星资料进行云顶高度的反演。但是，目前该方法中提到的两颗卫星同步观测是很难达到的，即使按预定的相同时间开始摄像，也很难保证卫星所覆盖的区域相同。因此，两颗卫星测量过程中的不同步，就会导致云顶高度的计算出现误差。

综上所述，目前常用的云高测量方法都存在一些局限性：激光云高仪只能在测试区域中有限的地区进行测量，而且一台仪器同一时刻只能测天空中云的一个点的云高，存在测量值的局限性问题；双站数字摄像测量云高方法需要采用图像采集仪器进行本地测量，测量仪器调整过程比较复杂，仪器的调整会极大的影响检测结果的精度；基于MODIS遥感图像CO2薄片测云高的方法，对低云高度的检测有一定局限性；双星立体观测云顶高度几何反演方法对于卫星立体观测的模拟资料可以得到较好的检测结果，对于真实数据的检测精度还未证实；其他的云高测试方法也有一些局限性。因此，研究精度高、可靠性好、经济实用、自动化高的探测方法具有重要的科学意义和应用价值。

发明内容

本发明的目的在于提供一种利用卫星遥感数据中的云影信息检测云高的方法，以解决目前常用的云高测量方法都存在局限性、精度低、可靠性差、成本高及自动化成度低等问题。

MODIS(中分辨率成像光谱仪)最大空间分辨率可达250米，扫描宽度2330公里，是当前世界上新一代“图谱合一”的光学遥感仪器。MODIS仪器与其它卫星相比，有以下特点和优势：1)空间分辨率大幅提高。空间分辨率提高了一个量级，由NOAA气象卫星的千米级提高到了MODIS的百米级。2)时间分辨率有优势。一天可过境4次，对各种突发性、快速变化的自然灾害有更强的实时监测能力。3)光谱分辨率大大提高。具有36个波段，这种多通道观测大大增强了对地球复杂系统的观测能力和对地表类型的识别能力。4)MODIS数据全球免费接收。该数据可以从国内、外MODIS卫星数据地面接收站免费接收。

本发明以MODIS遥感图像处理技术为基础，提出平原区域云高测量的新方法。其应用条件是分辨率为250m的MODIS数据卫星遥感数据中的平原区域，使用的MODIS数据是经过常规纠正处理后的数据。例如选取250m的MODIS数据卫星遥感数据MOD02QKM.A2002179.1820.005.2007148021815.hdf。

本发明的技术方案包括以下步骤：

一、插值计算：由于需要云和云阴影点的具体地理信息数据，包括经度、纬度、太阳天顶角、太阳方位角、卫星天顶角、卫星方位角数据，因此需要MODIS数据星历表MOD03的数据，即MOD03.A2002179.1820.005.2007147225005.hdf。针对250m分辨率的MOD03数据中的某些数据量重叠的情况，需要从1000m分辨率的数据中重新计算获得250m分辨率对应的MOD03数据，才能获得云和云阴影上点的更加精确的信息。因此这里采用通过采用线性插值公式解决该问题。线性内插是假设在二个已知数据中的变化为线性关系，因此可由已知二点a和c，去计算内插点b的值，公式如下：

$f\left(b\right)=f\left(a\right)+\frac{b-a}{c-a}(f\left(c\right)-f\left(a\right))---\left(1\right)$

其中a＜b＜c在上式的b点即是代表要内插的点，f(b)则是要计算的内插函数值。

具体实现步骤如下：

(1)读取1000m MOD03中的地理信息数据，如经度。

(2)1000m MOD03数据格式是2030行，1354列；插值后的250m MOD03数据格式是8120行，5416列，其行和列大小分别为1000m数据的4倍。采用将1000mMOD03数据中每10行×1354列的数据通过公式(1)插值成为40行×5416列的数据。由于1000m MOD03数据有2030行，因此需要进行203次的插值运算。

(3)利用Matlab软件中的INTERP2函数可以实现线性内插，通过计算得到250mMOD03中的地理信息数据。

二、云和云阴影检测：从原图像中手动选取具有云和云阴影的图像，选取图像在原图中的位置信息是起始点在3960行、4480列，即(3960，4480)，范围是40×40大小的图像范围，选取的波段是可见光波段2。该步骤主要完成从选取图像中检测云和云阴影图像的具体位置，选取数据是MOD02QKM。

目前，绝大多数的遥感图像云特征建模都是从图像的光谱特征出发，考虑云与其它地物在可见光、近红外和热红外波段上的差异。云在可见光和红外波段与植被、土壤以及水域等下垫面介质的反射率和辐射亮温值存在差异，云具有较高的反射率而具有低的亮温值等特点进行云与地物的区分。在MODIS数据的可见光波段0.645μm处，晴空数据地物一般具有较低的反射率，反射值在0～0.3之间，云则有高的反射率，一般高于0.3，是进行云检测处理的首选波段。1.375μm处于强的水汽吸收波段，地面的反射因水汽吸收难以到达传感器，导致反射率低，而高云上方的水汽对云的反射值削弱较少，因此该波段对高层卷云敏感，具有高的反射率，低云则具有低的反射率，通常低于0.02。此波段可以用来进行卷云的检测。8.6μm和11μm波段对云特性敏感。短波中有强的水汽吸收作用，导致晴空中8.6μm处的辐射值低于11μm处的辐射值；冰晶在8.6μm的吸收远小于11μm处的吸收，所以，冰晶组成的卷云在8.6μm处的辐射值高于11μm处的辐射值；而水滴组成的光学厚度较大的云则在8.6μm处的辐射值低于11μm处的辐射值。

采用多光谱云检测主要是利用以上各个波段对云的敏感特性，以及MODIS数据的各通道的波谱特性实现云与地物的区分。在多光谱云检测结果的基础上，对于可见光波段的图像，结合常用的图像检测方法如区域标定法，区域生长法，分水岭法，阈值分割法等将云与非云区域进行更有效的区分。对于云阴影的检测采取同云检测一样，首先利用云影在各个波段的敏感特性，主要利用1.2、0.94、0.87、0.66μm这四个波段反射率来检测云阴影；另外，利用3.9μm和4.0μm波段的亮度温度差大于零也可以完成云阴影的检测。

这里采用双阈值分割法在可见光波段2中实现云和云影图像的检测。图像阈值分割是一种广泛使用的图像分割技术，它利用了图像中要提取的目标与其背景在灰度特性上的差异，把图像视为具有不同灰度级的两类区域(目标和背景)的组合，选取一个合适的阈值，以确定图像中的每一个像素点应该属于目标还是背景区域，从而产生相应的二值图像。采用的双阈值分割的数学表达式为(2)式所示：

$f(i,j)=\left\{\begin{array}{cc}70& T1<f(i,j)<T2\\ 255& f(i,j)\&GreaterEqual;T2\\ 0& 0<f(i,j)\&le;T1\end{array}\right.---\left(2\right)$

通过对云和云影图像的光谱分析，选取云图像的阈值是T2，阴影图像的T1，f(i，j)是图像中像素的值。然后采用区域标定法，在双阈值检测结果的基础上，将云图像的位置和云阴影图像提取出来。

区域标定是给连接在一起的像素(连接成分)附上相同的标记，不同的连接成分附上不同的标记。区域标定在二值图像处理中占有非常重要的地位。通过该处理将各个连接成分区分开，然后就可以调查各个连接成分的形状特征。

区域标定也有许多方法，本发明采用了一种简单的方法。步骤如下：

(1)扫描图像，遇到没加标记的目标像素(白像素)P时，添加一个新的标记(LABEL)。

(2)将与P连接在一起(即相同连接成分)的像素添加相同的标记。

(3)进一步给所有与加标记像素连接在一起的像素添加相同的标记。直到连接在一起的像素全部被添加标记。这样，一个连接成分就被添加了相同的标记。

(4)回到第1步，重新查找新的没加标记的像素，重复上述各个步骤。图像全部被扫描后，处理结束。

采用区域标定后，会检测出很多候选目标。通过设定面积参数，使得云区域和云阴影区域能够从候选目标中提取出来，同时去掉图像中噪声对检测结果的影响。

三、MODIS数据的云和云影匹配：

根据云和云影检测的结果，进行云和云影的图像匹配。可以采用数字图像处理学中的图像匹配方法进行云和云影匹配，通过匹配得到云和云影之间的距离。

在云影图像中，通过云影检测，云和阴影区域已确定。每个云都有对应的影子图像，而且通常云阴影的图像要小于云的图像。基于模板匹配的原理，求得云与影间的距离的匹配方法为：以云阴影图像为模板，搜索云图像中与阴影图像匹配的位置。图像匹配中，通过采用自相关算法来计算模板与搜索区域的相关程度，将具有最大相关度的位置定为匹配最佳位置。如S为云阴影图像，C为云图像，M×N是云阴影图像的大小，则相关系数Coef(coefficient自定义的参数)计算公式如(3)所示，(i，j)为云和云阴影图像的匹配相对位置矢量，i为水平相对位置矢量，j为垂直相对位置矢量。

$\mathrm{Coef}(i,j)=\frac{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))}^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}}---\left(3\right)$

其中各变量为(4)式和(5)式：

$\stackrel{\&OverBar;}{S}(m,n)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}S(m,n)---\left(4\right)$

$\stackrel{\&OverBar;}{C}(m+i,n+j)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}C(m+i,n+j)---\left(5\right)$

通过公式(6)得到Coef最大时对应的最佳匹配相对位置矢量(i，j)：

$\underset{(i,j)}{\max }\mathrm{Coef}(i,j)---\left(6\right)$

将所有检测到的云图像的像素信息代入公式(5)的C，所有检测到的云阴影图像的像素信息代入S，根据公式(3)可以得到多个相关系数Coef。利用公式(6)选取最大的相关系数Coef，并据此得到最佳匹配相对位置矢量(i，j)。

四、建立云高模型求得云层高度：

根据云与云影之间的距离和卫星检测到的数据，通过云高模型的计算得到云层高度。由于空间中存在云，所以太阳的光线会被云完全或者部分遮挡，因此在地面上形成云的阴影图像。云阴影图像的位置与太阳相对于空间中云的位置有关，太阳-真云-云阴影共线；而当卫星观测地球时，空间中的云在MODIS探测器成像的位置与卫星的拍摄角度有关，即卫星-真云-投影云共线；因此MODIS数据中，云阴影图像和投影云图像都与空间中的云有着密切关系。可以根据云阴影图像、投影云图像和空间中的真云三者之间的空间位置关系建立模型，求得空间中的真云相对于地面的高度。

由于测试选择平原区域的云高测量，云和对应云阴影在同一平原区域，所以可以忽略地球是一个椭圆球体，并近似定义投影云和云阴影的地理信息在同一个平面上。

已知空间云上一点Point，太阳对其辐射形成的云阴影图像上对应的点是点Shadow，卫星对空间云中点Point投影后对应的点是点cloud。通过云和云影检测以及云和云影匹配处理，可以得到点shadow和点cloud之间的距离。从MOD03数据中读出点shadow和点cloud包含的地理信息，这里选取cloud点的位置信息是3971行、4502列，即(3971，4502)，则由最佳匹配相对位置(2，9)可知其对应的点shadow的坐标为(3969，4493)。根据这两点的位置从MOD03中读取点shadow的纬度lat1＝34.3604°、经度lon1＝-106.2759°，点cloud的纬度lat2＝34.3468°、经度lon2＝-106.2365°，将图像上云和云影像素之间的距离(单位：像素)转换为实际地理上的距离Distance(单位：公里)，所用计算公式为(7)-(9)：

$\mathrm{Dis}\tan \mathrm{ce}=2\arcsin \sqrt{{\sin }^2(a/2)+\cos \left(\mathrm{lat}1\right)\&times;\cos \left(\mathrm{lat}2\right)\&times;{\sin }^2(b/2)}\&times;6378.137---\left(7\right)$

其中：a＝lat1-lat2                (8)

b＝lon1-lon2                      (9)

经过公式计算，Distance＝3.9247(km)。由于测试选择平原地区的云和云阴影图像，因此可以忽略云和云阴影的高度差异，认为点shadow和点cloud是在同一平面上。从MOD03数据中读出点shadow的太阳方位角SunAZ＝134.905°＝2.3545rad、太阳天顶角SunZ＝14.954°＝0.261rad、点cloud的卫星方位角SensorAZ＝-77.2409°＝4.93508rad、卫星天顶角SensorZ＝41.2025°＝0.71912rad，将其结合云和云影距离等信息建立云高的空间几何模型，求得空间云中的一点Point相对于地面的高度Height。

令点Point为点1，其坐标为(x1，y1，z1)；点Shadow为点2，其坐标为(x2，y21，z2)；点Cloud为点3，其坐标为(x3，y3，z3)。

已知SunZ，SensorZ，因此A1，A2可知，即有公式(10)：

SunZ＝A1，SensorZ＝A2        (10)

已知SunAZ和Sensor AZ，因此B1，B2可求，即有公式(11)：

B1＝π-SunAZ，B2＝SensorAZ-π            (11)

令点1到点3的单位向量为

点1’的坐标为(x1’，y1’，z1’)，该向量在x、y、z轴的投影分别为

和则有公式(12)-(15)：

令点1到点2的向量为

该向量在x、y、z轴的投影分别为

和

且有公式(16)-(19)：

由向量内积公式有公式(20)：

$\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}=\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|\cos \mathrm{\&alpha;}---\left(20\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|=1---\left(21\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(22\right)$

α为

和

的夹角，0≤α≤π

则有： $\mathrm{\&alpha;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}}{\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|}\right)---\left(23\right)$

同理令点2到点3的单位向量为

点2’的坐标为(x2’，y2’，z2’)，该向量在x、y、z轴的投影分别为

和

则有公式(24)-(27)：

令点2到点1的向量为

该向量在x、y、z轴的投影分别为

和

且有公式(28)-(31)：

由向量内积公式有：

$\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}=\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|\cos \mathrm{\&beta;}---\left(32\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|=1---\left(33\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(34\right)$

β为

和的夹角，0≤β≤π

则有： $\mathrm{\&beta;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}}{\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|}\right)---\left(35\right)$

因此可得公式(36)：

γ＝π-α-β                    (36)

令D13为点1到点3的距离，令D23为点2到点3的距离，令D12为点1到点2的距离，由三角形正弦定理有

$\frac{D12}{\sin \mathrm{\&gamma;}}=\frac{D13}{\sin \mathrm{\&beta;}}=\frac{D23}{\sin \mathrm{\&alpha;}}---\left(37\right)$

且有D12＝‖r12‖                    (38)

则有：

$D13=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(39\right)$

$D23=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(40\right)$

点3的坐标3(x3，y3，z3)可用点1表示为：

x3＝x1+D13sinA1cosB1                (41)

y3＝y1+D13sinA1sinB1                (42)

z3＝z1+D13cosA1                     (43)

或者用点2表示为：

x3＝x2+D23sinA2cosB2                (44)

y3＝y2+D23sinA2sinB2                (45)

z3＝z2+D23cosA2                     (46)

理论上的云高Height的值计算为公式(47)：

Height＝z1+D13cosA1＝z2+D23cosA2    (47)

利用上述公式，将已知数据代入，可以求得云高。

五、云高计算和分析：寻找符合云高检测方法测试条件的时间和地区，确定测试地点；通过本发明提出的云高检测方法，将检测云的MODIS空间数据信息与云和云影之间的距离代入云高模型，对该检测云的多个点进行云高计算，并将计算结果求均值，作为检测云的云层高度，采用一至四的步骤，对周围9点进行云高计算，最后求得的云的云层高度的均值。

本发明的有益效果在于：可以利用MODIS遥感图像数据来实现满足检测条件的云高的检测，测试条件主要是MODIS卫星拍摄的平原区域中有云和云影的信息，从而完成对拍摄区域中单层云高的检测。该方法可以在不借助其它仪器的条件下，利用MODIS遥感图像数据的信息获得云层高度。与传统激光探测云高的方法相比，本发明中提出的方法可以实现从多点进行云高检测，进而提高云高检测结果的可靠性。该方法研究的最重要意义是可以完成对全球区域内满足测试条件的云进行高度测量，解决了云高测量区域局限性的问题。

附图说明：

图1是本发明实施例1选取出的云和云阴影图像；

图2是本发明实施例1检测得到的云和云阴影图像；

图3是本发明实施例1匹配的云和云阴影图像；

图4是本发明云高模型建立的理论基础图；

图5是本发明云高模型建立的几何模型；

图6是本发明云高的具体计算几何模型；

图7是本发明实施例2选取出的云和云阴影图像；

图8是本发明实施例2检测得到的云和云阴影图像；

图9是本发明实施例2匹配的云和云阴影图像。

具体实施方式：

实施例1：

本发明利用卫星遥感数据中的云影信息检测云高的方法，该方法的应用条件是分辨率为250m的MODIS卫星遥感数据中的平原区域的数据。

本发明以MODIS遥感图像处理技术为基础，提出平原区域云高测量的新方法。其应用条件是分辨率为250m的MODIS数据卫星遥感数据中的平原区域，使用的MODIS数据是经过常规纠正处理后的数据。例如选取250m的MODIS数据卫星遥感数据MOD02QKM.A2002179.1820.005.2007148021815.hdf。

本发明的技术方案包括以下步骤：

一、插值计算：由于需要云和云阴影点的具体地理信息数据，包括经度、纬度、太阳天顶角、太阳方位角、卫星天顶角、卫星方位角数据，因此需要MODIS数据星历表MOD03的数据，即MOD03.A2002179.1820.005.2007147225005.hdf。针对250m分辨率的MOD03数据中的某些数据量重叠的情况，需要从1000m分辨率的数据中重新计算获得250m分辨率对应的MOD03数据，才能获得云和云阴影上点的更加精确的信息。因此这里采用通过采用线性插值公式解决该问题。线性内插是假设在二个已知数据中的变化为线性关系，因此可由已知二点a和c，去计算内插点b的值，公式如下：

$f\left(b\right)=f\left(a\right)+\frac{b-a}{c-a}(f\left(c\right)-f\left(a\right))---\left(1\right)$

其中a＜b＜c在上式的b点即是代表要内插的点，f(b)则是要计算的内插函数值。

具体实现步骤如下：

(1)读取1000m MOD03中的地理信息数据，如经度。

(2)1000m MOD03数据格式是2030行，1354列；插值后的250m MOD03数据格式是8120行，5416列，其行和列大小分别为1000m数据的4倍。采用将1000mMOD03数据中每10行×1354列的数据通过公式(1)插值成为40行×5416列的数据。由于1000m MOD03数据有2030行，因此需要进行203次的插值运算。

(3)利用Matlab软件中的INTERP2函数可以实现线性内插，通过计算得到250mMOD03中的地理信息数据。

二、云和云阴影检测：从原图像中手动选取具有云和云阴影的图像，选取图像在原图中的位置信息是起始点在3960行、4480列，即(3960，4480)，范围是40×40大小的图像范围，选取的波段是可见光波段2。具体选取图像如图1所示。该步骤主要完成从选取图像中检测云和云阴影图像的具体位置，选取数据是MOD02QKM。

目前，绝大多数的遥感图像云特征建模都是从图像的光谱特征出发，考虑云与其它地物在可见光、近红外和热红外波段上的差异。云在可见光和红外波段与植被、土壤以及水域等下垫面介质的反射率和辐射亮温值存在差异，云具有较高的反射率而具有低的亮温值等特点进行云与地物的区分。在MODIS数据的可见光波段0.645μm处，晴空数据地物一般具有较低的反射率，反射值在0～0.3之间，云则有高的反射率，一般高于0.3，是进行云检测处理的首选波段。1.375μm处于强的水汽吸收波段，地面的反射因水汽吸收难以到达传感器，导致反射率低，而高云上方的水汽对云的反射值削弱较少，因此该波段对高层卷云敏感，具有高的反射率，低云则具有低的反射率，通常低于0.02。此波段可以用来进行卷云的检测。8.6μm和11μm波段对云特性敏感。短波中有强的水汽吸收作用，导致晴空中8.6μm处的辐射值低于11μm处的辐射值；冰晶在8.6μm的吸收远小于11μm处的吸收，所以，冰晶组成的卷云在8.6μm处的辐射值高于11μm处的辐射值；而水滴组成的光学厚度较大的云则在8.6μm处的辐射值低于11μm处的辐射值。

采用多光谱云检测主要是利用以上各个波段对云的敏感特性，以及MODIS数据的各通道的波谱特性实现云与地物的区分。在多光谱云检测结果的基础上，对于可见光波段的图像，结合常用的图像检测方法如区域标定法，区域生长法，分水岭法，阈值分割法等将云与非云区域进行更有效的区分。对于云阴影的检测采取同云检测一样，首先利用云影在各个波段的敏感特性，主要利用1.2、0.94、0.87、0.66μm这四个波段反射率来检测云阴影；另外，利用3.9μm和4.0μm波段的亮度温度差大于零也可以完成云阴影的检测。

这里采用双阈值分割法在可见光波段2中实现云和云影图像的检测。图像阈值分割是一种广泛使用的图像分割技术，它利用了图像中要提取的目标与其背景在灰度特性上的差异，把图像视为具有不同灰度级的两类区域(目标和背景)的组合，选取一个合适的阈值，以确定图像中的每一个像素点应该属于目标还是背景区域，从而产生相应的二值图像。采用的双阈值分割的数学表达式为(2)式所示：

$f(i,j)=\left\{\begin{array}{cc}70& T1<f(i,j)<T2\\ 255& f(i,j)\&GreaterEqual;T2\\ 0& 0<f(i,j)\&le;T1\end{array}\right.---\left(2\right)$

通过对云和云影图像的光谱分析，选取云图像的阈值是T2＝0.315μm，阴影图像的T1＝0.134μm，f(i，j)是图像中像素的值。然后采用区域标定法，在双阈值检测结果的基础上，将云图像的位置和云阴影图像提取出来。

区域标定是给连接在一起的像素(连接成分)附上相同的标记，不同的连接成分附上不同的标记。区域标定在二值图像处理中占有非常重要的地位。通过该处理将各个连接成分区分开，然后就可以调查各个连接成分的形状特征。

区域标定也有许多方法，本实施例采用了一种简单的方法。步骤如下：

(1)扫描图像，遇到没加标记的目标像素(白像素)P时，添加一个新的标记(LABEL)。

(2)将与P连接在一起(即相同连接成分)的像素添加相同的标记。

(3)进一步给所有与加标记像素连接在一起的像素添加相同的标记。直到连接在一起的像素全部被添加标记。这样，一个连接成分就被添加了相同的标记。

(4)回到第1步，重新查找新的没加标记的像素，重复上述各个步骤。图像全部被扫描后，处理结束。

采用区域标定后，会检测出很多候选目标。通过设定面积参数，使得云区域和云阴影区域能够从候选目标中提取出来，同时去掉图像中噪声对检测结果的影响。经过云和云阴影检测后的图像如图2所示。其中灰色的像素代表背景，黑色的像素代表云阴影图像，白色的像素代表云图像。从图像可见，该方法可以较好地从原图像中检测出云和云影图像。

三、MODIS数据的云和云影匹配：

根据云和云影检测的结果，进行云和云影的图像匹配。可以采用数字图像处理学中的图像匹配方法进行云和云影匹配，通过匹配得到云和云影之间的距离。

在云影图像中，通过云影检测，云和阴影区域已确定。每个云都有对应的影子图像，而且通常云阴影的图像要小于云的图像。基于模板匹配的原理，求得云与影间的距离的匹配方法为：以云阴影图像为模板，搜索云图像中与阴影图像匹配的位置。图像匹配中，通过采用自相关算法来计算模板与搜索区域的相关程度，将具有最大相关度的位置定为匹配最佳位置。如S为云阴影图像，C为云图像，M×N是云阴影图像的大小，则相关系数Coef(coefficient自定义的参数)计算公式如(3)所示，(i，j)为云和云阴影图像的匹配相对位置矢量，i为水平相对位置矢量，j为垂直相对位置矢量。

$\mathrm{Coef}(i,j)=\frac{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))}^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}}---\left(3\right)$

其中各变量为(4)式和(5)式：

$\stackrel{\&OverBar;}{S}(m,n)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}S(m,n)---\left(4\right)$

$\stackrel{\&OverBar;}{C}(m+i,n+j)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}C(m+i,n+j)---\left(5\right)$

通过公式(6)得到Coef最大时对应的最佳匹配相对位置矢量(i，j)：

$\underset{(i,j)}{\max }\mathrm{Coef}(i,j)---\left(6\right)$

将所有检测到的云图像的像素信息代入C，所有检测到的云阴影图像的像素信息代入S，根据公式(3)可以得到多个相关系数Coef。利用公式(6)选取最大的相关系数Coef，并据此得到最佳匹配相对位置矢量(i，j)。通过计算，最大的相关系数是0.95，对应的(i，j)为(2，9)，即云和云影图像的对应匹配点水平方向相差9个像素，垂直方向相差2个像素。匹配效果如图3所示。

四、建立云高模型求得云层高度：

根据云与云影之间的距离和卫星检测到的数据，通过云高模型的计算得到云层高度。图4为云高模型建立的理论基础图。由于空间中存在云，所以太阳的光线会被云完全或者部分遮挡，因此在地面上形成云的阴影图像。云阴影图像的位置与太阳相对于空间中云的位置有关，太阳-真云-云阴影共线；而当卫星观测地球时，空间中的云在MODIS探测器成像的位置与卫星的拍摄角度有关，即卫星-真云-投影云共线；因此MODIS数据中，云阴影图像和投影云图像都与空间中的云有着密切关系。可以根据云阴影图像、投影云图像和空间中的真云三者之间的空间位置关系建立模型，求得空间中的真云相对于地面的高度。

图5是云高模型建立的几何模型：由于测试选择平原区域的云高测量，投影云和对应云阴影在同一平原区域，所以可以忽略地球是一个椭圆球体，并近似定义云和云阴影的地理信息在同一个平面上。

已知空间云上一点Point，太阳对其辐射形成的云阴影图像上对应的点是点Shadow，卫星对空间云中点Point投影后对应的点是点cloud。通过云和云影检测以及云和云影匹配处理，可以得到点shadow和点cloud之间的距离。从MOD03数据中读出点shadow和点cloud包含的地理信息，这里选取cloud点的位置信息是3971行、4502列，即(3971，4502)，则由最佳匹配相对位置(2，9)可知其对应的点shadow的坐标为(3969，4493)。根据这两点的位置从MOD03中读取点shadow的纬度lat1＝34.3604°、经度lon1＝-106.2759°，点cloud的纬度lat2＝34.3468°、经度lon2＝-106.2365°，将图像上云和云影像素之间的距离(单位：像素)转换为实际地理上的距离Distance(单位：公里)，所用计算公式为(7)-(9)：

$\mathrm{Dis}\tan \mathrm{ce}=2\arcsin \sqrt{{\sin }^2(a/2)+\cos \left(\mathrm{lat}1\right)\&times;\cos \left(\mathrm{lat}2\right)\&times;{\sin }^2(b/2)}\&times;6378.137---\left(7\right)$

其中：a＝lat1-lat2                    (8)

b＝lon1-lon2                          (9)

经过公式计算，Distance＝3.9247(km)。由于测试选择平原地区的云和云阴影图像，因此可以忽略云和云阴影的高度差异，认为点shadow和点cloud是在同一平面上。从MOD03数据中读出点shadow的太阳方位角SunAZ＝134.905°＝2.3545rad、太阳天顶角SunZ＝14.954°＝0.261rad、点cloud的卫星方位角SensorAZ＝-77.2409°＝4.93508rad、卫星天顶角SensorZ＝41.2025°＝0.71912rad，将其结合云和云影距离等信息建立云高的空间几何模型，求得空间云中的一点Point相对于地面的高度Height。

为了便于表达云高模型，将图5的模型表示为图6。令点Point为点1，其坐标为(x1，y1，z1)；点Shadow为点2，其坐标为(x2，y2，z2)；点Cloud为点3，其坐标为(x3，y3，z3)。

已知SunZ，SensorZ，因此A1，A2可知，即有公式(10)：

SunZ＝A1，SensorZ＝A2          (10)

已知SunAZ和Sensor AZ，因此B1，B2可求，即有公式(11)：

B1＝π-SunAZ，B2＝SensorAZ-π                (11)

令点1到点3的单位向量为点1’的坐标为(x1’，y1’，z1’)，该向量在x、y、z轴的投影分别为

和

则有公式(12)-(15)：

令点1到点2的向量为

该向量在x、y、z轴的投影分别为

和

且有公式(16)-(19)：

由向量内积公式有公式(20)：

$\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}=\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|\cos \mathrm{\&alpha;}---\left(20\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|=1---\left(21\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(22\right)$

α为

和

的夹角，0≤α≤π

则有： $\mathrm{\&alpha;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}}{\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|}\right)---\left(23\right)$

同理令点2到点3的单位向量为

点2’的坐标为(x2’，y2’，z2’)，该向量在x、y、z轴的投影分别为

和则有公式(24)-(27)：

令点2到点1的向量为

该向量在x、y、z轴的投影分别为

和

且有公式(28)-(31)：

由向量内积公式有：

$\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}=\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|\cos \mathrm{\&beta;}---\left(32\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|=1---\left(33\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(34\right)$

β为

和

的夹角，0≤β≤π

则有： $\mathrm{\&beta;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}}{\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|}\right)---\left(35\right)$

因此可得公式(36)：

γ＝π-α-β                    (36)

令D13为点1到点3的距离，令D23为点2到点3的距离，令D12为点1到点2的距离，由三角形正弦定理有

$\frac{D12}{\sin \mathrm{\&gamma;}}=\frac{D13}{\sin \mathrm{\&beta;}}=\frac{D23}{\sin \mathrm{\&alpha;}}---\left(37\right)$

且有D12＝‖r12‖

则有：

$D13=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(39\right)$

$D23=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(40\right)$

点3的坐标3(x3，y3，z3)可用点1表示为：

x3＝x1+D13sinA1cosB1                (41)

y3＝y1+D13sinA1sinB1                (42)

z3＝z1+D13cosA1                     (43)

或者用点2表示为：

x3＝x2+D23sinA2cosB2                (44)

y3＝y2+D23sinA2sinB2                (45)

z3＝z2+D23cosA2                     (46)

理论上的云高Height的值计算为公式(47)：

Height＝z1+D13cosA1＝z2+D23cosA2    (47)

利用上述公式，将已知数据代入，可以求得云高Height＝3.547(km)。

五、云高计算和分析：寻找符合云高检测方法测试条件的时间和地区，确定测试地点；通过本发明提出的云高检测方法，将检测云的MODIS空间数据信息与云和云影之间的距离代入云高模型，对该检测云的多个点进行云高计算，并将计算结果求均值，作为检测云的云层高度，采用一至四的步骤，对周围9点进行云高计算，最后求得的云的云层高度的均值是3.554km。

实施例2：

本发明利用卫星遥感数据中的云影信息检测云高的方法，该方法的应用条件是分辨率为250m的MODIS卫星遥感数据中的平原区域的数据。

本发明以MODIS遥感图像处理技术为基础，提出平原区域云高测量的新方法。其应用条件是分辨率为250m的MODIS数据卫星遥感数据中的平原区域，使用的MODIS数据是经过常规纠正处理后的数据。例如选取250m的MODIS数据卫星遥感数据MOD02QKM.A2002179.1820.005.2007148021815.hdf。

一、插值计算：由于需要云和云阴影点的具体地理信息数据，包括经度、纬度、太阳天顶角、太阳方位角、卫星天顶角、卫星方位角数据，因此需要MODIS数据星历表MOD03的数据，即MOD03.A2002179.1820.005.2007147225005.hdf。针对250m分辨率的MOD03数据中的某些数据量重叠的情况，需要从1000m分辨率的数据中重新计算获得250m分辨率对应的MOD03数据，才能获得云和云阴影上点的更加精确的信息。因此这里采用通过采用线性插值公式解决该问题。线性内插是假设在二个已知数据中的变化为线性关系，因此可由已知二点a和c，去计算内插点b的值，公式如下：

$f\left(b\right)=f\left(a\right)+\frac{b-a}{c-a}(f\left(c\right)-f\left(a\right))---\left(1\right)$

其中a＜b＜c在上式的b点即是代表要内插的点，f(b)则是要计算的内插函数值。

具体实现步骤如下：

(1)读取1000m MOD03中的地理信息数据，如经度。

(2)1000m MOD03数据格式是2030行，1354列；插值后的250m MOD03数据格式是8120行，5416列，其行和列大小分别为1000m数据的4倍。采用将1000mMOD03数据中每10行×1354列的数据通过公式(1)插值成为40行×5416列的数据。由于1000m MOD03数据有2030行，因此需要进行203次的插值运算。

(3)利用Matlab软件中的INTERP2函数可以实现线性内插，通过计算得到250mMOD03中的地理信息数据。

二、云和云阴影检测：从原图像中手动选取具有云和云阴影的图像，选取图像在原图中的位置信息是起始点在1930行、2985列，即(1930，2985)，范围是30×30大小的图像范围，选取的波段是可见光波段2。具体选取图像如图7所示。该步骤主要完成从选取图像中检测云和云阴影图像的具体位置，选取数据是MOD02QKM。

目前，绝大多数的遥感图像云特征建模都是从图像的光谱特征出发，考虑云与其它地物在可见光、近红外和热红外波段上的差异。云在可见光和红外波段与植被、土壤以及水域等下垫面介质的反射率和辐射亮温值存在差异，云具有较高的反射率而具有低的亮温值等特点进行云与地物的区分。在MODIS数据的可见光波段0.645μm处，晴空数据地物一般具有较低的反射率，反射值在0～0.3之间，云则有高的反射率，一般高于0.3，是进行云检测处理的首选波段。1.375μm处于强的水汽吸收波段，地面的反射因水汽吸收难以到达传感器，导致反射率低，而高云上方的水汽对云的反射值削弱较少，因此该波段对高层卷云敏感，具有高的反射率，低云则具有低的反射率，通常低于0.02。此波段可以用来进行卷云的检测。8.6μm和11μm波段对云特性敏感。短波中有强的水汽吸收作用，导致晴空中8.6μm处的辐射值低于11μm处的辐射值；冰晶在8.6μm的吸收远小于11μm处的吸收，所以，冰晶组成的卷云在8.6μm处的辐射值高于11μm处的辐射值；而水滴组成的光学厚度较大的云则在8.6μm处的辐射值低于11μm处的辐射值。

采用多光谱云检测主要是利用以上各个波段对云的敏感特性，以及MODIS数据的各通道的波谱特性实现云与地物的区分。在多光谱云检测结果的基础上，对于可见光波段的图像，结合常用的图像检测方法如区域标定法，区域生长法，分水岭法，阈值分割法等将云与非云区域进行更有效的区分。对于云阴影的检测采取同云检测一样，首先利用云影在各个波段的敏感特性，主要利用1.2、0.94、0.87、0.66μm这四个波段反射率来检测云阴影；另外，利用3.9μm和4.0μm波段的亮度温度差大于零也可以完成云阴影的检测。

这里采用双阈值分割法在可见光波段2中实现云和云影图像的检测。图像阈值分割是一种广泛使用的图像分割技术，它利用了图像中要提取的目标与其背景在灰度特性上的差异，把图像视为具有不同灰度级的两类区域(目标和背景)的组合，选取一个合适的阈值，以确定图像中的每一个像素点应该属于目标还是背景区域，从而产生相应的二值图像。采用的双阈值分割的数学表达式为(2)式所示：

$f(i,j)=\left\{\begin{array}{cc}70& T1<f(i,j)<T2\\ 255& f(i,j)\&GreaterEqual;T2\\ 0& 0<f(i,j)\&le;T1\end{array}\right.---\left(2\right)$

通过对云和云影图像的光谱分析，选取云图像的阈值是T2＝0.382μm，阴影图像的T1＝0.136μm，f(i，j)是图像中像素的值。然后采用区域标定法，在双阈值检测结果的基础上，将云图像的位置和云阴影图像提取出来。

区域标定是给连接在一起的像素(连接成分)附上相同的标记，不同的连接成分附上不同的标记。区域标定在二值图像处理中占有非常重要的地位。通过该处理将各个连接成分区分开，然后就可以调查各个连接成分的形状特征。

区域标定也有许多方法，本发明采用了一种简单的方法。步骤如下：

(1)扫描图像，遇到没加标记的目标像素(白像素)P时，添加一个新的标记(LABEL)。

(2)将与P连接在一起(即相同连接成分)的像素添加相同的标记。

(3)进一步给所有与加标记像素连接在一起的像素添加相同的标记。直到连接在一起的像素全部被添加标记。这样，一个连接成分就被添加了相同的标记。

(4)回到第1步，重新查找新的没加标记的像素，重复上述各个步骤。图像全部被扫描后，处理结束。

采用区域标定后，会检测出很多候选目标。通过设定面积参数，使得云区域和云阴影区域能够从候选目标中提取出来，同时去掉图像中噪声对检测结果的影响。经过云和云阴影检测后的图像如图8所示。其中灰色的像素代表背景，黑色的像素代表云阴影图像，白色的像素代表云图像。从图像可见，该方法可以较好地从原图像中检测出云和云影图像。

三、MODIS数据的云和云影匹配：

根据云和云影检测的结果，进行云和云影的图像匹配。可以采用数字图像处理学中的图像匹配方法进行云和云影匹配，通过匹配得到云和云影之间的距离。

在云影图像中，通过云影检测，云和阴影区域已确定。每个云都有对应的影子图像，而且通常云阴影的图像要小于云的图像。基于模板匹配的原理，求得云与影间的距离的匹配方法为：以云阴影图像为模板，搜索云图像中与阴影图像匹配的位置。图像匹配中，通过采用自相关算法来计算模板与搜索区域的相关程度，将具有最大相关度的位置定为匹配最佳位置。如S为云阴影图像，C为云图像，M×N是云阴影图像的大小，则相关系数Coef(coefficient自定义的参数)计算公式如(3)所示，(i，j)为云和云阴影图像的匹配相对位置矢量，i为水平相对位置矢量，j为垂直相对位置矢量。

$\mathrm{Coef}(i,j)=\frac{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))}^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}}---\left(3\right)$

其中各变量为(4)式和(5)式：

$\stackrel{\&OverBar;}{S}(m,n)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}S(m,n)---\left(4\right)$

$\stackrel{\&OverBar;}{C}(m+i,n+j)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}C(m+i,n+j)---\left(5\right)$

通过公式(6)得到Coef最大时对应的最佳匹配相对位置矢量(i，j)：

$\underset{(i,j)}{\max }\mathrm{Coef}(i,j)---\left(6\right)$

将所有检测到的云图像的像素信息代入C，所有检测到的云阴影图像的像素信息代入S，根据公式(3)可以得到多个相关系数Coef。利用公式(6)选取最大的相关系数Coef，并据此得到最佳匹配相对位置矢量(i，j)。通过计算，最大的相关系数是0.96，对应的(i，j)为(4，8)，即云和云影图像的对应匹配点水平方向相差8个像素，垂直方向相差4个像素。匹配效果如图9所示。

四、建立云高模型求得云层高度：

根据云与云影之间的距离和卫星检测到的数据，通过云高模型的计算得到云层高度。图4为云高模型建立的理论基础图。由于空间中存在云，所以太阳的光线会被云完全或者部分遮挡，因此在地面上形成云的阴影图像。云阴影图像的位置与太阳相对于空间中云的位置有关，太阳-真云-云阴影共线；而当卫星观测地球时，空间中的云在MODIS探测器成像的位置与卫星的拍摄角度有关，即卫星-真云-投影云共线；因此MODIS数据中，云阴影图像和投影云图像都与空间中的云有着密切关系。可以根据云阴影图像、投影云图像和空间中的真云三者之间的空间位置关系建立模型，求得空间中的真云相对于地面的高度。

图5是云高模型建立的几何模型：由于测试选择平原区域的云高测量，投影云和对应云阴影在同一平原区域，所以可以忽略地球是一个椭圆球体，并近似定义云和云阴影的地理信息在同一个平面上。

已知空间云上一点Point，太阳对其辐射形成的云阴影图像上对应的点是点Shadow，卫星对空间云中点Point投影后对应的点是点cloud。通过云和云影检测以及云和云影匹配处理，可以得到点shadow和点cloud之间的距离。从MOD03数据中读出点shadow和点cloud包含的地理信息，这里选取cloud点的位置信息是1942行、2998列，即(1942，2998)，则由最佳匹配相对位置(4，8)可知其对应的点shadow的坐标为(1938，2990)。根据这两点的位置从MOD03中读取点shadow的纬度lat1＝39.7052°、经度lon1＝-109.8255°，点cloud的纬度lat2＝39.6927°、经度lon2＝-109.8047°，将图像上云和云影像素之间的距离(单位：像素)转换为实际地理上的距离Distance(单位：公里)，所用计算公式为(7)-(9)：

$\mathrm{Dis}\tan \mathrm{ce}=2\arcsin \sqrt{{\sin }^2(a/2)+\cos \left(\mathrm{lat}1\right)\&times;\cos \left(\mathrm{lat}2\right)\&times;{\sin }^2(b/2)}\&times;6378.137---\left(7\right)$

其中：a＝lat1-lat2                    (8)

b＝lon1-lon2                          (9)

经过公式计算，Distance＝2.2606(km)。由于测试选择平原地区的云和云阴影图像，因此可以忽略云和云阴影的高度差异，认为点shadow和点cloud是在同一平面上。从MOD03数据中读出点shadow的太阳方位角SunAZ＝137.1709°＝2.3941rad、太阳天顶角SunZ＝20.9475°＝0.3656rad、点cloud的卫星方位角SensorAZ＝-80.4722°＝4.8787rad、卫星天顶角SensorZ＝6.4925°＝0.1133rad，将其结合云和云影距离等信息建立云高的空间几何模型，求得空间云中的一点Point相对于地面的高度Height。

为了便于表达云高模型，将图5的模型表示为图6。令点Point为点1，其坐标为(x1，y1，z1)；点Shadow为点2，其坐标为(x2，y2，z2)；点Cloud为点3，其坐标为(x3，y3，z3)。

已知SunZ，SensorZ，因此A1，，A2可知，即有公式(10)：

SunZ＝A1，SensorZ＝A2                            (10)

已知SunAZ和Sensor AZ，因此B1，B2可求，即有公式(11)：

B1＝π-SunAZ，B2＝SensorAZ-π                    (11)

令点1到点3的单位向量为点1’的坐标为(x1’，y1’，z1’)，该向量在x、y、z轴的投影分别为

和

则有公式(12)-(15)：

令点1到点2的向量为

该向量在x、y、z轴的投影分别为

和

且有公式(16)-(19)：

由向量内积公式有公式(20)：

$\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}=\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|\cos \mathrm{\&alpha;}---\left(20\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\left|\right|=1---\left(21\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(22\right)$

α为和

的夹角，0≤α≤π

则有： $\mathrm{\&alpha;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{11}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}}{\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|}\right)---\left(23\right)$

同理令点2到点3的单位向量为点2’的坐标为(x2’，y2’，z2’)，该向量在x、y、z轴的投影分别为

和则有公式(24)-(27)：

令点2到点1的向量为

该向量在x、y、z轴的投影分别为

和且有公式(28)-(31)：

由向量内积公式有：

$\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}=\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|\cos \mathrm{\&beta;}---\left(32\right)$

其中

$\left|\right|\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\left|\right|=1---\left(33\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(34\right)$

β为

和

的夹角，0≤β≤π

则有： $\mathrm{\&beta;}=\arccos \left(\frac{\stackrel{\&RightArrow;}{{r_{22}}^{\&prime;}}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}}{\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|}\right)---\left(35\right)$

因此可得公式(36)：

γ＝π-α-β                (36)

令D13为点1到点3的距离，令D23为点2到点3的距离，令D12为点1到点2的距离，由三角形正弦定理有

$\frac{D12}{\sin \mathrm{\&gamma;}}=\frac{D13}{\sin \mathrm{\&beta;}}=\frac{D23}{\sin \mathrm{\&alpha;}}---\left(37\right)$

且有D12＝‖r12‖            (38)

则有：

$D13=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(39\right)$

$D23=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(40\right)$

点3的坐标3(x3，y3，z3)可用点1表示为：

x3＝x1+D13sinA1cosB1                (41)

y3＝y1+D13sinA1sinB1                (42)

z3＝z1+D13cosA1                     (43)

或者用点2表示为：

x3＝x2+D23sinA2cosB2                (44)

y3＝y2+D23sinA2sinB2                (45)

z3＝z2+D23cosA2                     (46)

理论上的云高Height的值计算为公式(47)：

Height＝z1+D13cosA1＝z2+D23cosA2    (47)

利用上述公式，将已知数据代入，可以求得云高Height＝4.7181(km)。

五、云高计算和分析：寻找符合云高检测方法测试条件的时间和地区，确定测试地点；通过本发明提出的云高检测方法，将检测云的MODIS空间数据信息与云和云影之间的距离代入云高模型，对该检测云的多个点进行云高计算，并将计算结果求均值，作为检测云的云层高度，采用一至四的步骤，对周围9点进行云高计算，最后求得的云的云层高度的均值是4.6736km。

Claims (2)

1、一种利用卫星遥感数据中的云影信息检测云高的方法，该方法的应用条件是分辨率为250m的MODIS卫星遥感数据中的平原区域的数据，该方法包括以下四个部分：

一、插值计算：线性内插是由已知二点a和c，去计算内插点b的值，公式如下：

$f\left(b\right)=f\left(a\right)+\frac{b-a}{c-a}(f\left(c\right)-f\left(a\right))---\left(1\right)$

其中a＜b＜c，上式的b点是代表要内插的点，f(b)则是要计算的内插函数值；

具体实现步骤如下：

1)读取1000m MOD03中的地理信息数据；

2)1000m MOD03数据格式是2030行，1354列；插值后的250m MOD03数据格式是8120行，5416列，其行和列大小分别为1000m数据的4倍；采用将1000m MOD03数据中每10行×1354列的数据通过公式(1)插值成为40行x5416列的数据；由于1000m MOD03数据有2030行，因此需要进行203次的插值运算；

3)利用Matlab软件中的INTERP2函数实现线性内插，通过计算得到250m MOD03中的地理信息数据；

二、云和云阴影检测：

从气象卫星提供的数据中选取具有云和云阴影的图像；

采用双阈值分割法在可见光波段2中实现云和云影图像的检测，采用的双阈值分割的数学表达式为(2)式所示：

$f(i,j)=\left\{\begin{array}{cc}70& T1<f(i,j)<T2\\ 255& f(i,j)\&GreaterEqual;T2\\ 0& 0<f(i,j)\&le;T1\end{array}\right.---\left(2\right)$

通过对云和云影图像的光谱分析，选取云图像的阈值是T2，阴影图像的T1，f(i，j)是图像中像素的值；然后采用区域标定法，在双阈值检测结果的基础上，将云图像的位置和云阴影图像提取出来；

区域标定法的步骤如下：

(a)扫描图像，遇到白像素P时，添加一个新的标记LABEL；

(b)将与白像素P连接在一起，即相同连接成分的像素添加相同的标记；

(c)进一步给所有与加标记像素连接在一起的像素添加相同的标记，直到连接在一起的像素全部被添加标记；这样，一个连接成分就被添加了相同的标记；

(d)回到第(a)步，重新查找新的白像素，重复上述各个步骤；图像全部被扫描后，处理结束；

采用区域标定后，对检测出的候选目标通过设定面积参数，使得云区域和云阴影区域能够从候选目标中提取出来，同时去掉图像中噪声对检测结果的影响；

三、MODIS数据的云和云影匹配：

根据云和云影检测的结果，采用数字图像处理学中的图像匹配方法进行云和云影匹配，通过匹配得到云和云影之间的距离；

具体步骤如下：

设S为云阴影图像，C为云图像，MxN是云阴影图像的大小，则相关系数Coef计算公式如(3)所示，(i，j)为云和云阴影图像的匹配相对位置矢量，i为水平相对位置矢量，j为垂直相对位置矢量；

$\mathrm{Coef}(i,j)=\frac{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}{\sqrt{\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{(S(m,n)-\stackrel{\&OverBar;}{S}(m,n))}^2\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}(C(m+i,n+j)-\stackrel{\&OverBar;}{C}(m+i,n+j))}}---\left(3\right)$

其中各变量为(4)式和(5)式：

$\stackrel{\&OverBar;}{S}(m,n)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}S(m,n)---\left(4\right)$

$\stackrel{\&OverBar;}{C}(m+i,n+j)=\frac{1}{M\&times;N}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}C(m+i,n+j)---\left(5\right)$

通过公式(6)得到Coef最大时对应的最佳匹配相对位置矢量(i，j)：

$\underset{(i,j)}{\max }\mathrm{Coef}(i,j)---\left(6\right)$

根据上述公式进行云和云阴影的匹配，将所有检测到的云图像的像素信息代入公式(5)的C，所有检测到的云阴影图像的像素信息代入S，根据公式(3)可以得到多个相关系数Coef；利用公式(6)选取最大的相关系数Coef，并据此得到最佳匹配相对位置矢量(i，j)；

四、建立云高模型求得云层高度：

根据云与云影之间的距离和卫星检测到的数据，通过云高模型的计算得到云层高度；根据云阴影图像、投影云图像和空间中的真云三者之间的空间位置关系建立模型，求得空间中的真云相对于地面的高度；

已知空间云上一点Point，太阳对其辐射形成的云阴影图像上对应的点是点Shadow，卫星对空间云中点Point投影后对应的点是点cloud；通过云和云影检测以及云和云影匹配处理，可以得到点shadow和点cloud之间的距离；从MOD03数据中读出点shadow和点cloud包含的地理信息：点shadow的纬度lat1、经度lon1，点cloud的纬度lat2、经度lon2，将图像上云和云影像素之间的单位为像素的距离，转换为实际地理上的单位为公里的距离Distance，所用计算公式为(7)-(9)：

$\mathrm{Dis}\tan \mathrm{ce}=2\arcsin \sqrt{{\sin }^2(a/2)+\cos \left(\mathrm{lat}1\right)\&times;\cos \left(\mathrm{lat}2\right)\&times;{\sin }^2(b/2)}\&times;6378.137---\left(7\right)$

其中：a＝lat1-lat2         (8)

b＝lon1-lon2               (9)

从MOD03数据中读出点shadow的太阳方位角SunAZ、太阳天顶角SunZ、点cloud的卫星方位角SensorAZ、卫星天顶角SensorZ，将其结合云和云影距离等信息建立云高的空间几何模型，求得空间云中的一点Point相对于地面的高度Height；

令点Point为点1，其坐标为(x1，y1，z1)；点Shadow为点2，其坐标为(x2，y2，z2)；点Cloud为点3，其坐标为(x3，y3，z3)；

已知SunZ，SensorZ，因此A1，A2可知，即有公式(10)：

SunZ＝A1，SensorZ＝A2      (10)

已知SunAZ和Sensor AZ，因此B1，B2可求，即有公式(11)：

B1＝π-SunAZ，B2＝SensorAZ-π   (11)

令点1到点3的单位向量为

点1’的坐标为(x1’，y1’，z1’)，该向量在x、y、z轴的投影分别为

和

则有公式(12)-(15)：

令点1到点2的向量为该向量在x、y、z轴的投影分别为

和且有公式(16)-(19)：

由向量内积公式有公式(20)：

${\stackrel{\&RightArrow;}{r_{11}}}^{\&prime;}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}=\left|\right|{\stackrel{\&RightArrow;}{r_{11}}}^{\&prime;}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|\cos \mathrm{\&alpha;}---\left(20\right)$

其中

$\left|\right|{\stackrel{\&RightArrow;}{r_{11}}}^{\&prime;}\left|\right|=1---\left(21\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(22\right)$

α为

和

的夹角，0≤α≤π，则有：

$\mathrm{\&alpha;}=\arccos \left(\frac{{\stackrel{\&RightArrow;}{r_{11}}}^{\&prime;}\&CenterDot;\stackrel{\&RightArrow;}{r_{12}}}{\left|\right|\stackrel{\&RightArrow;}{r_{12}}\left|\right|}\right)---\left(23\right)$

同理令点2到点3的单位向量为

点2’的坐标为(x2’，y2’，z2’)，该向量在x、y、z轴的投影分别为

和

则有公式(24)-(27)：

令点2到点1的向量为

该向量在x、y、z轴的投影分别为

和且有公式(28)-(31)：

由向量内积公式有：

${\stackrel{\&RightArrow;}{r_{22}}}^{\&prime;}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}=\left|\right|{\stackrel{\&RightArrow;}{r_{22}}}^{\&prime;}\left|\right|\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|\cos \mathrm{\&beta;}---\left(32\right)$

其中

$\left|\right|{\stackrel{\&RightArrow;}{r_{22}}}^{\&prime;}\left|\right|=1---\left(33\right)$

$\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|=\sqrt{{(x1-x2)}^2+{(y1-y2)}^2+{(z1-z2)}^2}---\left(34\right)$

β为和

的夹角，0≤β≤π，则有：

$\mathrm{\&beta;}=\arccos \left(\frac{{\stackrel{\&RightArrow;}{r_{22}}}^{\&prime;}\&CenterDot;\stackrel{\&RightArrow;}{r_{21}}}{\left|\right|\stackrel{\&RightArrow;}{r_{21}}\left|\right|}\right)---\left(35\right)$

因此可得公式(36)：

γ＝π-α-β    (36)

令D13为点1到点3的距离，令D23为点2到点3的距离，令D12为点1到点2的距离，由三角形正弦定理有

$\frac{D12}{\sin \mathrm{\&gamma;}}=\frac{D13}{\sin \mathrm{\&beta;}}=\frac{D23}{\sin \mathrm{\&alpha;}}---\left(37\right)$

且有D12＝‖r12‖      (38)

则有：

$D13=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(39\right)$

$D23=\frac{D12\sin \mathrm{\&beta;}}{\sin \mathrm{\&gamma;}}---\left(40\right)$

已知点3的坐标3(x3，y3，z3)可用点1表示为：

x3＝x1+D13sinA1cosB1          (41)

y3＝y1+D13sinA1sinB1          (42)

z3＝z1+D13cosA1               (43)

或者用点2表示为：

x3＝x2+D23sinA2cosB2          (44)

y3＝y2+D23sinA2sinB2          (45)

z3＝z2+D23cosA2               (46)

云高Height的值计算为公式(47)：

Height＝z1+D13cosA1＝z2+D23cosA2    (47)

利用上述公式，将已知数据代入，求得云高。

2、根据权利要求1所述的利用卫星遥感数据中的云影信息检测云高的方法，其特征在于：该方法包括云高计算和分析，将检测云的MODIS空间数据信息与云和云影之间的距离代入云高模型，对该检测云的多个点进行云高计算，并将计算结果求均值。

Images