本发明涉及遥感领域，特别是涉及一种光学遥感图像估算一年生海冰厚度方法。
背景技术：
：高纬度近海海域是一年生海冰的多发区(如渤海、波罗的海)，尽管一年生海冰厚度大多为几厘米至几十厘，但是每年冬季的海冰给海上交通安全、海洋养殖、海上石油平台造成了严重危害。另一方面，中国研究人员正在开展的海冰淡化研究可以通过脱盐处理使海冰(微咸水)转化为淡水，为沿岸地区的社会经济发展提供一种新的淡水资源。因此，不论是从灾害方面考虑，还是从资源方面衡量考虑将海冰变害为宝，及时掌握海冰面积和厚度等海冰参数信息对于局地海-气热量交换机理、海冰工程设计、海冰灾害防灾减灾和海冰资源化利用都具有极其重要的意义。遥感技术可以大范围地、快速地获取渤海海冰电磁波反射、辐射或散射图像，通过这些图像可以进一步估算或推测诸如海冰面积、厚度等海冰参数。从上个世纪80年代开始，遥感技术开始被应用于海冰面积提取和海冰厚度估算，包括多光谱遥感图像、热红外遥感图像、微波遥感图像和高光谱遥感图像先后被广泛应用于一年生海冰面积和厚度估算中。国内外学者建立了许多用于估算海冰厚度的模型。然而，这些以光学、热红外和微波遥感图像为数据源，用于估算渤海海冰厚度的模型定量化略显不足，偏定性和经验。很多模型中的参数需要凭经验确定，不具有广泛适用性；或者模型中的参数随时间和空间而变化，不具有稳定性；抑或确定模型中的参数不具可操作性，导致模型不够实用。因此，为了提高一年生海冰估算精度，开发出具有可操作性、实用化、高精度的一年生海冰厚度估算模型是一个亟待解决的技术问题，具有重要意义。技术实现要素：针对现有的技术空白和缺点，本发明所要解决的技术问题是提供一种大范围、快速、准确的一年生海冰厚度估算方法。本发明具体采用的技术方案如下：光学遥感图像估算一年生海冰厚度的方法，包括以下步骤：s1.收集目标海域光学遥感图像和地面标定数据；所述的目标海域光学遥感图像中包含e或f波段图像，以及a、b、c和d波段图像；所述a、b、c、d、e、f波段的中心波长分别为560±20nm、650±20nm、850±10nm、1240±10nm、1640±10nm、2130±10nm；s2.对一年生海冰厚度估算模型中的参数进行标定，其中所述的一年生海冰厚度估算模型为：其中，h为海冰厚度；ρ∞,b为海冰厚度达到或超过光学厚度时b波段的反射率；ρw,b为海水在b波段的反射率；ρb、ρc和ρd分别为海冰在b、c和d波段的反射率；ρe||f为海冰在e或f波段的反射率；a、b和c为需要进行标定的参数；s3.对目标海域光学遥感图像进行预处理，获得其反射率图像；s4.提取目标海域的海冰范围，获得海冰掩膜；s5.根据s4所得海冰掩膜和s3所得反射率图像，用s2所标定的一年生海冰厚度估算模型，在海冰掩膜覆盖范围内逐像元计算海冰厚度，获得海冰厚度分布图。作为优选，所述s1中，所述的目标海域光学遥感图像可为modis影像。作为优选，所述s1中，地面标定数据包括野外晴空测量的400-2500nm范围内且光谱采样间隔不低于1nm的海冰和海水反射光谱数据及对应的冰厚数据。作为优选，所述s1中，地面标定数据的数据量在20个以上。作为优选，所述s2中，ρ∞,b取值为0.94。作为优选，所述s2中，ρw,b根据有冰像元邻近的无冰像元在b波段的反射率插值模拟得到。作为优选，所述s2中，标定方法如下：将地面标定数据中测量的反射光谱经传感器光谱响应函数卷积后，得到海冰和海水各个波段反射率数据，将其与测量的对应厚度数据代入海冰厚度估算模型，用最小二乘法计算得到该海域所对应的参数a，b和c。作为优选，所述步骤s3中，光学遥感图像进行预处理是指对光学遥感图像进行辐射定标、几何校正和大气校正，最终获得反射率图像。作为优选，所述s4中，提取海冰范围时采用归一化差分水体指数(normalizeddifferencewaterindex，ndwi)，ndwi计算方法如下：其中，ρa和ρb分别为海冰在a和b波段的反射率；若ndwi不大于预设阈值，则判定该像元为海冰；否则判定该像元为海水。作为优选，所述的ndwi的预设阈值为0.52。本发明以辐射传输模型为依据，对复杂的海冰辐射传输模型进行了适当的简化，并提出了可靠的参数标定方式。本发明可以大大减小由海冰表面反射和海冰消光系数的空间异质性所带来的海冰估算误差，提高光学遥感图像估算渤海海冰厚度的精度。附图说明图1遥感影像示意图；图2本方法估算的海冰厚度分布图；图3对比例1方法估算的海冰厚度分布图；图4对比例2方法估算的海冰厚度分布图；图5为光学遥感图像估算一年生海冰厚度的方法流程示意图。具体实施方式为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。如图5所示，一种光学遥感图像估算一年生海冰厚度的方法，包括以下步骤：s1.收集目标海域光学遥感图像和地面标定数据。目标海域光学遥感图像中包含e或f波段图像，以及a、b、c和d波段图像；所述a、b、c、d、e、f波段的中心波长分别为560±20nm、650±20nm、850±10nm、1240±10nm、1640±10nm、2130±10nm。光学遥感图像可为modis影像，也可以是其他具有这些波段的影像。地面标定数据包括野外晴空测量的400-2500nm范围内且光谱采样间隔不低于1nm的海冰和海水反射光谱数据及对应的冰厚数据。数据量应在20个以上。s2.对一年生海冰厚度估算模型中的参数进行标定，其中一年生海冰厚度估算模型为：其中，h为海冰厚度；ρ∞,b为海冰厚度达到或超过光学厚度时b波段的反射率，此处取值可为0.94；ρw,b为海水在b波段的反射率，可根据有冰像元邻近的无冰像元在b波段的反射率插值模拟得到；ρb、ρc和ρd分别为海冰在b、c和d波段的反射率；ρe||f为海冰在e或f波段的反射率；a、b和c为需要进行标定的参数。a、b和c的标定可采用如下方法：将地面标定数据中测量的反射光谱经传感器光谱响应函数卷积后，得到海冰和海水各个波段反射率数据，将其与测量的对应厚度数据代入海冰厚度估算模型，用最小二乘法计算得到该海域所对应的参数a，b和c。s3.对目标海域光学遥感图像进行辐射定标、几何校正和大气校正等预处理，获得其反射率图像；s4.提取目标海域的海冰范围，获得海冰掩膜。提取海冰范围时采用归一化差分水体指数(normalizeddifferencewaterindex，ndwi)，ndwi计算方法如下：其中，ρa和ρb分别为海冰在a和b波段的反射率；其中a波段的中心波长为560±20nm；若ndwi不大于预设阈值(可设置为0.52)，则判定该像元为海冰；否则判定该像元为海水：式中：tn为判定是否为海冰所采用的ndwi的预设阈值。该阈值可设为0.52。s5.根据s4所得海冰掩膜和s3所得反射率图像，用s2所标定的一年生海冰厚度估算模型，在海冰掩膜覆盖范围内逐像元计算海冰厚度，获得海冰厚度分布图。下面将上述方法应用于具体实施例中，以便本领域技术人员能够更好地理解本发明的效果。实施例为了验证本发明的效果，选择中国渤海作为目标海域。本实施例的实现方法如前所述，按照以下步骤实施渤海海冰厚度估算：1)收集渤海海域中分辨率成像光谱仪(moderate-resolutionimagingspectroradiometer，modis)的第1、2、4、5、7波段图像，成像日期为2010年1月29日，成像时间为北京时间11:00。所对应中心波长分别为645nm、859nm、555nm、1240nm和2130nm。选择野外晴空用asd测量渤海海冰和海水在400-2500nm范围内的反射光谱数据和对应的冰厚数据，asd光谱仪的输出光谱采样间隔为1nm。数据量共计37个，完成标定数据的收集；2)对一年生海冰厚度估算模型中的参数进行标定。一年生海冰厚度估算模型为：其中，ρ∞,1为海冰厚度达到或超过光学厚度时时modis第1波段的反射率，为0.94；ρw,1为海水在modis第1波段的反射率，根据邻近的无冰像元在modis第1波段的反射率插值模拟得到；ρ1、ρ2、ρ5和ρ7为分别为海冰在modis第1、2、5和7波段反射率；a、b和c为需要进行定标的参数。a、b和c的定标方法如下：将测量的反射光谱经modis光谱响应函数卷积后的海冰和海水各个波段反射率及其对应厚度代入海冰厚度估算模型，用最小二乘法计算得到该海域所对应的参数分别为a＝0.1182，b＝0.0132和c＝0.0075，由此完成估算模型中参数的定标。3)对光学遥感图像进行辐射定标、几何校正和大气校正等预处理，获得反射率图像，图2为对应的灰度图像；4)提取海域中的海冰范围，获得海冰掩膜。本实施例中，提取海冰像元时，采用归一化差分水体指数(normalizeddifferencewaterindex，ndwi)，ndwi的计算方法为：其中，ρ2和ρ4分别为modis的第2和4波段反射率。海冰像元提取方法为：其中，tn为阈值，阈值为0.52；获得海冰掩膜数据。5)根据4)所得海冰掩膜和3)所得反射率图像，用2)所标定的一年生海冰厚度估算模型，在海冰掩膜覆盖范围内逐像元计算海冰厚度，获得海冰厚度分布图(图2)。下面给出若干对比例，以展示本发明与其他常用的海冰厚度估算方法在效果上的差异。对比例11)收集与实施例相同的影像，即渤海海域中modis的第1、2、3、4、5和7波段图像，成像日期为2010年1月29日，成像时间为北京时间11:00。2)对光学遥感图像进行预处理，通过辐射定标和几何校正后获得大气上界反射率ρk*，其中k＝1,2,3,4,5,7；再通过大气校正获得反射率图像ρk，其中k＝1,2,3,4,5,7；3)采用双波段比值法提取海冰像元，双波段比值法公式如下：其中，和分别为modis第1和2波段的大气上界反射率。选择合适的阈值提取海冰，获得海冰掩膜数据；4)计算宽波段反照率α，其计算公式如下：α＝0.160ρ1+0.291ρ2+0.243ρ3+0.116ρ4+0.112ρ5+0.081ρ7-0.0015其中，ρ1，ρ2，ρ3，ρ4，ρ5和ρ7为分别为海冰在modis第1、2、3、4、5和7波段的反射率。5)计算海冰厚度h，计算方法为其中，αmax为海冰厚度达到或超过光学厚度时的宽波段反照率，本对比例中设置为0.7；αsea为海水的宽波段反照率，本比对例中设置为0.06；μα为宽波段衰减系数，本比对例中设置为1.209。逐像元计算后获得海冰厚度分布图(图4)。对比例21)收集与实施例相同的影像，即渤海海域中modis的第1和2波段图像，成像日期为2010年1月29日，成像时间为北京时间11:00；2)对光学遥感图像进行预处理，进行辐射定标、几何校正和大气校正，获得反射率图像ρk，其中k＝1,2,3,4,5,7；3)采用modis第一波段分区确定阈值提取海冰像元，获得海冰掩膜；4)计算海冰反照率α，计算公式如下：α＝0.423ρ1+0.577ρ2其中，ρ1和ρ2为分别为海冰在modis第1和2波段反射率。5)计算海冰厚度h，计算公式如下：其中，αmax为海冰厚度达到或超过光学厚度时的海冰反照率，本比对例中设置为0.7；αsea为海水反照率，本比对例中设置为0.1；μα为衰减系数，本比对例中设置为1.209。逐像元计算后获得海冰厚度分布图(图5)。实施例和对比例1、比对例2的部分像元计算结果比对如表1所示。其中，实施例的rmse为5.0cm，对比例1和对比例2的rmse分别为15.0和5.8cm。实施例的误差小于对比例误差，相比对比例1误差减小了10cm(66％)，相比对比例2误差减小了0.8cm(14％)。对比图3-5可以证明，相比对比例1(2-167cm)和对比例2(2-265cm)的海冰厚度范围，实施例(2-42cm)所估算的海冰厚度分布范围更加符合实际情况。综上，实施例估算海冰厚度的方法更加准确和稳定。表1部分像元计算结果比对编号实测厚度/cm实施例/cm实施例误差/cm对比例1/cm对比例1误差/cm对比例2/cm对比例2误差/cm117.520.83.335.317.821.84.3212.520.17.630.417.919.16.6320.018.1-1.925.65.69.9-10.1420.019.8-0.233.113.120.40.4515.021.16.135.620.621.76.7617.520.73.234.717.222.65.1722.021.1-0.937.115.124.92.9832.021.1-10.937.35.324.4-7.6932.020.8-11.236.04.023.1-8.91030.020.6-9.436.16.123.1-6.91120.020.70.735.115.123.03.01225.022.1-2.936.811.824.3-0.71325.021.8-3.236.211.222.8-2.21425.022.8-2.242.617.628.53.51527.522.0-5.541.614.127.2-0.31620.022.92.944.624.629.79.71720.022.12.145.825.831.911.91825.021.2-3.835.610.621.9-3.11920.022.12.140.320.325.95.92025.020.9-4.134.09.020.2-4.82124.020.3-3.733.69.622.2-1.82220.019.6-0.432.412.419.9-0.12322.020.8-1.233.811.820.6-1.4rmse——5.0—15.0—5.8以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关
技术领域：
的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。当前第1页12