一种基于雷达后向散射系数时间序列的石漠化判断方法与流程

本发明属于石漠化治理领域，尤其是一种基于雷达后向散射系数时间序列的石漠化判断方法。

背景技术

石漠化地区往往多云多雨，而传统光学遥感图像虽然具备了丰富的光谱信息，但极易受到时间、天气等多重因素的影响，这就对石漠化地区数据的获取制造了相当大的困难。雷达影像具有全天时全天候不受天气条件约束的成像特点，影像纹理信息丰富，地物轮廓清晰，但本身只具有单波段信息，两种方法均不利于对象地区是否属于石漠化类型的判断。所以，利用多时相雷达影像提取石漠化信息，可以克服许多不利因素，有效弥补各自的缺点。

技术实现要素：

本发明要解决的技术问题是：提供一种基于雷达后向散射系数时间序列的石漠化判断方法，以解决现有技术传统光学遥感图像极易受到时间、天气等多重因素的影响，雷达影像只具有单波段信息，两种方法均不利于对象地区是否属于石漠化类型的判断的问题。

本发明采用的技术方案为：一种基于雷达后向散射系数时间序列的石漠化判断方法，包括如下步骤：

步骤一：选取对象一年内的雷达数据，每个月份至少取一份雷达数据，且雷达数据的极化方式保持一致，以确保某地区石漠化信息提取结果的精确性；

步骤二：为了对雷达数据的原始振幅信息进行定量表达，对步骤一中的雷达数据进行辐射定标处理，辐射定标处理采用的公式为：

pd为后向散射强度，pt为发射功率，pn为附着功率，是接收天线关于在仰角θel和方位角θaz下的功率，g^e为系统的电流增益，gp为雷达处理器的常数，r为雷达影像在传播过程中造成的损耗，ls为系统损耗，la为大气损耗，a为雷达影像的散射面积，σ⁰为雷达亮度系数；

步骤三：选出步骤二中辐射定标处理的雷达数据所在区域的同分辨率数字高程模型数据，并对该数据进行地形辐射校正，其校正公式为：

式中：σ⁰为校正后成像单元的后向散射系数，λ⁰为初始成像单元的后向散射系数，β1为距离向坡度角，β2为方位向坡度角，δ为影像入射角，θ为影像反射角；

步骤四：采用精致lee算子方法对步骤三中校正后的雷达数据进行影像斑点噪声的滤波处理；

步骤五：将步骤四中经滤波处理后的雷达数据按照相同的像元大小分割成m行×n列像元构成的数据，则m×n的值为步骤一中雷达数据的份数，然后按照时间顺序合成一幅影像，使得不同时相影像中的像元值一一对应，影像中的x轴为经度，y轴为纬度，z轴为时间轴，那么每个像元的后向散射系数则构成了后向散射系数的时间序列，雷达数据后向散射系数时间序列表达式为：

t＝{[(x1,y1),z1],[(x2,y2),z2],…[(xi,yi),zi]…,[(xm×n,ym×n),zm×n]}

式中：(xi，yi)为第i个像元的经纬坐标值，zi是基于第i个像元构建的后向散射系数时间序列；

步骤六：利用dtw算法对步骤五中构建的时间序列进行石漠化信息提取，得到dtw值，记为d(i，j)。

步骤七：用目视解译法选取步骤五中雷达数据所在地区的典型石漠化区域边缘位置的混合像元，计算混合像元在与雷达数据各个时间节点相对应的后向散射系数值的平均值，这样就构建了一条混合像元的时间序列曲线，最后计算混合像元时间序列和石漠化像元标准时间序列的dtw值，以该值作为雷达数据所在区域像元类型判断的阈值，记为d1(i，j)，若d(i，j)＜d1(i，j)，则该像元属于石漠化类型，若d(i，j)≥d1(i，j)，则该像元不属于石漠化类型。

步骤四中，精致lee算子方法的算法原理是：若计算结果显示雷达数据像元值大小较为一致，那么该地区属于匀质区域，此时可以利用envi中的低通滤波方法进行斑点噪声滤波，若计算结果显示雷达数据像元值大小差异较大，则采用envi中的高通滤波方法进行斑点噪声滤波，一般我们认为如果像元值的方差小于5则可以认为差异较小；像元值的方差大于等于5的话则可以认为差异较大。

步骤五中，相同的像元大小为100-200。

步骤六中，dtw算法是指dynamictimewarping，动态时间调整算法，该算法提取石漠化信息的步骤如下：首先，提取典型石漠化像元的后向散射系数时间序列构建石漠化像元标准时间序列，接着计算雷达数据所在区域内像元的时间序列与石漠化地物标准时间序列之间的dtw距离，其算法公式为：

d(i,j)＝min{d(i-1,j-1),d(i-1,j),d(i,j-1),d(i-1,j)}，

式中：i＝2，…，m；j＝2，3，…，n；d(i，j)为雷达数据所在区域像元的时间序列曲线与石漠化像元标准时间序列距离的最小累加值，即为dtw值。

本发明与现有技术相比的优点在于：

1、本方法对雷达数据进行辐射定标处理、地形辐射校正处理和影像斑点噪声的滤波处理，克服了单时相雷达数据由于受到斑点噪声影像而导致关键时间节点信息缺失的缺陷；

2、与光学遥感提取石漠化信息相比，本方法选取雷达数据作为原始数据克服了光学遥感数据提取石漠化信息时极易受到时间、天气等多重因素影响的缺陷，提高了石漠化信息判断的精度。

附图说明

图1实施例1的石漠化类型判断结果图。

具体实施方式

下面结合附图对本发明方法的实施方式做详细说明：

实施例1

一种基于雷达后向散射系数时间序列的石漠化判断方法，包括如下步骤：

步骤一：选取对象一年内的雷达数据，每个月份至少取一份雷达数据，且雷达数据的极化方式保持一致，以确保某地区石漠化信息提取结果的精确性；

选择贵州省毕节市七星关区内经纬跨度为104°2′00″至104°49′26〞e、26°27′58〞至26°50′08〞n的地区(以下简称对象区域)作为发明方案实施地点。雷达数据选择c波段sar雷达数据，数据的获取时间分别为2016年1月9日、2月5日、3月3日、4月16日、5月20日、6月21日、7月17日、8月19日、9月21日、10月23日、11月18日以及12月19日，数据极化方式均为vv。

步骤二：为了对雷达数据的原始振幅信息进行定量表达，对步骤一中的雷达数据进行辐射定标处理，辐射定标处理采用的公式为：

pd为后向散射强度，pt为发射功率，pn为附着功率，是接收天线关于在仰角θel和方位角θaz下的功率，g^e为系统的电流增益，gp为雷达处理器的常数，r为雷达影像在传播过程中造成的损耗，ls为系统损耗，la为大气损耗，a为雷达影像的散射面积，σ⁰为雷达亮度系数；

打开envi软件，载入步骤一中的12幅对象区域的雷达数据，选择radiometriccorrection中的applygainandoffset工具，在弹出的gainandoffsetvalues对话框中输入辐射定标参数，从而获取12幅雷达数据的后向散射系数值。

步骤三：选出步骤二中辐射定标处理的雷达数据所在区域的同分辨率数字高程模型数据，并对该数据进行地形辐射校正，其校正公式为：

式中：σ⁰为校正后成像单元的后向散射系数，λ⁰为初始成像单元的后向散射系数，β1为距离向坡度角，β2为方位向坡度角，δ为影像入射角，θ为影像反射角；

在envi软件中载入对象区域的aster数字高程模型作为校正模板。在envi软件中选择工具箱中的empiricallinecomputefactorsandcorrect工具，对弹出的对话框中选择校正模板文件，接着输入校正参数完成地形辐射校正。

步骤四：采用精致lee算子方法对步骤三中校正后的雷达数据进行影像斑点噪声的滤波处理；

在matlab软件中运行精致lee算子代码对对象区域雷达数据进行数据斑点噪声的滤波处理。计算结果像元值的方差小于5，显示对象区域雷达数据像元值大小较为一致，此时可以利用envi软件中的低通滤波方法进行斑点噪声滤波。通过对对象区域雷达数据运用idl中的“smooth”函数可以完成低通滤波处理。

步骤五：将步骤四中经滤波处理后的雷达数据按照相同的像元大小分割成m行×n列像元构成的数据，则m×n的值为步骤一中雷达数据的份数，然后按照时间顺序合成一幅影像，使得不同时相影像中的像元值一一对应，影像中的x轴为经度，y轴为纬度，z轴为时间轴，那么每个像元的后向散射系数则构成了后向散射系数的时间序列，雷达数据后向散射系数时间序列表达式为：

t＝{[(x1,y1),z1],[(x2,y2),z2],…[(xi,yi),zi]…,[(xm×n,ym×n),zm×n]}

式中：(xi，yi)为第i个像元的经纬坐标值，zi是基于第i个像元构建的后向散射系数时间序列；

打开envi软件中的layerstacking工具，把每一个时序当做一个通道进行layerstacking，完成构建基于对象区域雷达数据后向散射系数的时间序列处理。

步骤六：利用dtw算法对步骤五中构建的时间序列进行石漠化信息提取，得到dtw值，记为d(i，j)；

步骤七：用目视解译法选取步骤五中雷达数据所在地区的典型石漠化区域边缘位置的一定数目像元，本实施例选取50个，记为混合像元，计算混合像元在与雷达数据各个时间节点相对应的后向散射系数值的平均值，这样就构建了一条混合像元的时间序列曲线，最后计算混合像元时间序列和石漠化像元标准时间序列的dtw值，以该值作为雷达数据所在区域像元类型判断的阈值，记为d1(i，j)，若d(i，j)＜d1(i，j)，则该像元属于石漠化类型，若d(i，j)≥d1(i，j)，则该像元不属于石漠化类型。

计算对象区域混合像元时间序列和对象区域石漠化像元标准时间序列的dtw值，经计算，d1(i，j)的值为47。若对象区域内像元的时间序列和石漠化像元标准时间序列的dtw值小于47，则该像元属于石漠化类型；若像元的时间序列和石漠化像元标准时间序列的dtw值大于等于47，则该像元不属于石漠化类型。结果显示：对象区域内石漠化面积为3457.62公顷，占对象区域比例为55.31％；非石漠化面积为2793.17公顷，占对象区域比例为44.69％，如图1所示。

步骤四中，精致lee算子方法的算法原理是：若计算结果显示雷达数据像元值大小较为一致，那么该地区属于匀质区域，此时可以利用envi中的低通滤波方法进行斑点噪声滤波，若计算结果显示雷达数据像元值大小差异较大，则采用envi中的高通滤波方法进行斑点噪声滤波，一般我们认为如果像元值的方差小于5则可以认为差异较小；像元值的方差大于等于5的话则可以认为差异较大。

步骤五中，相同的像元大小为100-200，本实施例取为150。

步骤六中，dtw算法是指dynamictimewarping，动态时间调整算法，该算法提取石漠化信息的步骤如下：首先，提取典型石漠化像元的后向散射系数时间序列构建石漠化像元标准时间序列，接着计算雷达数据所在区域内像元的时间序列与石漠化地物标准时间序列之间的dtw距离，其算法公式为：

d(i,j)＝min{d(i-1,j-1),d(i-1,j),d(i,j-1),d(i-1,j)}，

式中：i＝2，…，m；j＝2，3，…，n；d(i，j)为雷达数据所在区域像元的时间序列曲线与石漠化像元标准时间序列距离的最小累加值，即为dtw值。