一种基于分形维数和极化分解的月表平坦区域选取方法

1.本发明涉及雷达图像处理领域，尤其是一种基于分形维数和极化分解的月表平坦区域选取方法。


背景技术：

2.合成孔径雷达（synthetic aperture radar, sar）是一种主动式微波成像系统，具有全天时全天候以及穿透成像的工作能力，被广泛的用于地球遥感与行星探测。极化合成孔径雷达提供了除目标散射回波的幅度和相位特性外的极化散射信息，为推断行星表面的物理特性提供了一种有效的途径。极化目标分解通过对表征目标极化散射特征的矩阵进行分解，提取出具有物理意义的各种极化参数，是极化sar图像解译的一个重要分支。yamaguchi四元素分解是一种非相干极化分解方法，在地面目标的散射机制识别与分类中应用广泛。
3.分形理论是描述不规则几何形态的有力数学工具，其中分形维数（d）是基于分形的方法测量表面粗糙度的一个重要参数。对于二维表面而言，其值在2.0和3.0之间。光滑的二维曲面的分形维数值为2.0，随着表面粗糙度的增加，分形维数值增加，对于非常粗糙的曲面，分形维数值可达到3.0。因此，分形维数可以用来描述不同的土地覆盖特征，并且分形维数和地质表面的纹理特征（或粗糙度）之间存在一一对应的关系。
4.平坦地区的选取在月球探测中具有重要需求。例如，在月球极区永久阴影区的水冰探测中，需要寻找较为平坦的区域以排除表面粗糙度的干扰，并且为未来永久阴影区着陆器的着陆和采样提供支持。此外，月球表面的平坦地区还可作为分布式目标用于未来月球极化sar数据的定标。目前已有的寻找月球平坦区域的方法有：通过高分辨率光学图像（月球轨道光学相机或月球车载相机）反演得到月表粗糙度信息；由月球轨道激光高度计数据反演得到月球表面的数字高程模型（dem），进而得到平坦地区。然而光学图像无法提供月球极区永久阴影区的信息，并且月球车影像覆盖范围小，只覆盖月球表面极少数感兴趣的区域。此外，受限于固有的激光斑点间距，已有的月球表面数字高程模型（dem）数据在一些感兴趣区域分辨率较低。


技术实现要素：

5.为了提取出月球表面感兴趣地区的平坦区域，本发明提出了一种联合分形维数（空间几何特征）和极化分解技术（极化散射特征）的极化sar图像平坦区域选取方法。可以联合利用分形维度和极化分解技术，通过识别月球表面的纹理特征与极化散射特征来进行平坦区域的选取。本发明提出的方法能够在不依赖额外dem数据的情况下，利用极化sar数据实现自然地质表面的平坦地区选取，为将来全极化sar数据在月球的应用研究提供支撑。
6.该方法可分为两个主要步骤：粗选取与精选取。首先是粗选取阶段：利用雷达后向散射总功率图像计算分形维数，然后提取出分形维数小于分形维数均值的区域。精细选取：对极化sar图像进行yamaguchi四元素分解，在粗选取出的候选区域内，进一步筛选出单次
散射功率占比最大的区域。具体的，本发明的技术方案为：一种基于分形维数和极化分解的月表平坦区域选取方法，包括如下步骤：步骤1、获取全极化sar数据；步骤2、根据雷达后向散射总功率计算分形维数d，得到分形维数图；步骤3、对分形维数图进行二值化，大于分形维数均值的像素值置为白色，小于均值的置为黑色，得到二值化后的分形维数图；步骤4、对二值化后的分形维数图，筛选出所有满足不包含白色像素的连续n
×
n像素区域，作为粗选取出的候选区域，n为像素区域宽度；步骤5、利用yamaguchi四元素分解方法对sar回波进行极化分解，获取单次散射分量、二次散射分量、体散射分量以及螺旋散射分量的大小；步骤6、针对每一个候选区域，筛选出单次散射分量占比最大的n
×
n像素区域，即为精选取出的平坦地区，其中n 《 n。
7.有益效果：本发明针对月球轨道sar获取的全极化回波数据，首先对后向散射总功率图像计算分形维数，然后筛选出小于分形维数均值的连续n
×
n像素区域作为候选区域；接着利用yamaguchi四元素分解方法对雷达回波进行极化分解；随后针对每一个候选区域，筛选出单次散射分量占比最大的连续n
×
n像素区域，作为最终选取的平坦地区。本发明具有三个优点：（1）针对月球两极永久阴影区等光学传感器无法提供观测信息的区域，提供了一种利用极化sar数据进行平坦地区选取的方法；（2）综合利用地质表面的纹理特征和极化散射特征，显著提升了平坦地区选取的准确性；（3）为未来极化sar数据在月球的应用研究提供支撑。此外，本发明利用高分三号卫星在沙漠地区的全极化sar数据以及srtm1 dem数据对所提方法进行了验证实验，结果表明，在不依赖额外dem数据的情况下，本发明提出的方法能够准确提取所选区域的地质特征，有效实现平坦地区的识别。
8.附图说明
9.图1为本发明的方法流程图；图2为三棱柱法表面积计算示意图；图3为局部三棱柱法不同栅格步长示意图；（a）步长为1；（b）步长为2；（c）步长为4；（d）步长为8；图4（a）后向散射总功率图；图4（b）分形维数图，分形维数d大于其均值的像素置为白色，小于均值的像素置为黑色；图5为搜索到的符合条件的n
×
n候选区域（白色框标识）；（a）白色实线框标识出分形维数图中搜索到的200
×
200像素区域；（b）分形维数放大图；（c）后向散射总功率放大图；图6为最终筛选出的40
×
40像素（320m
×
320m）平坦区域（白色虚线框），白色实线
矩形框为根据分形维数筛选出的粗选取区域。
10.图7为整个sar图像场景内，全部40
×
40像素窗口的高程标准差与最大高程偏差的散点密度图。
具体实施方式
11.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅为本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域的普通技术人员在不付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。
12.为达到上述目的，本发明的技术方案如图1所示，一种基于分形维数和极化分解的月表平坦区域选取方法，包括如下步骤：步骤1、获取全极化sar数据；步骤2、根据雷达后向散射总功率计算分形维数d，得到分形维数图；步骤3、对分形维数图进行二值化，大于分形维数均值的像素值置为白色，小于均值的置为黑色，得到二值化后的分形维数图；步骤4、对二值化后的分形维数图，筛选出所有满足不包含白色像素的连续n
×
n像素区域，作为粗选取出的候选区域，n为像素区域宽度；步骤5、利用yamaguchi四元素分解方法对sar回波进行极化分解，获取单次散射分量、二次散射分量、体散射分量以及螺旋散射分量的大小；步骤6、针对每一个候选区域，筛选出单次散射分量占比最大的n
×
n像素区域，即为精选取出的平坦地区，其中n《n。
13.根据本发明的实施例，所述的步骤2，根据雷达后向散射总功率计算分形维数，得到分形维数图；具体如下：首先采用局部三棱柱法计算雷达后向散射总功率图像的分形维数。该方法将图像上每个像素的灰度值看作“高度”值，然后利用一个9
×
9的窗口在图像上移动，在每一个窗口中心的位置上，计算不同栅格步长下的表面积。表面积的计算通过引入中心像素点来进行，中心像素点将栅格正方形分割成4个三角形，而中心像素点的像素值为4个角点像素的平均值，分别计算这4个三角形的面积之和即为栅格的表面积，如图2所示。a、b、c、d为栅格像素点，a’，b’，c’，d’为相应的灰度值，p为正方形中心点，p点像素值为4个角点像素值的均值。
14.改变栅格步长（步长分别为1，2，4，8），重复上述计算，从而得到栅格表面积与步长之间的对应关系，如图3所示，为局部三棱柱法不同栅格步长示意图。从左至右，（a）（b）（c）（d）步长依次为1,2,4,8。在对数坐标系下（以2为底），对得到的表面积与步长进行最小二乘拟合，获得拟合直线的斜率k，然后利用下式即可求出窗口中心位置上的分形维数d：d=2-k（1）之后将窗口移动至下一像素，重复上述过程，直至整个图像计算完毕。
15.根据本发明的一个实施例，下面以高分三号在沙漠地区的全极化sar数据的后向散射总功率图像为例来计算分形维数。所选sar数据的尺寸为6197
×
8035像素，数据编号为：gf3_kas_qpsi_007178_e83.8_n39.9_20171221_l1a_ahv_l10002865872。该图像的分辨率为8m，中心经度为83
°
47'21"e，纬度为39
°
56'36"n。
16.步骤3、对分形维数图进行二值化，大于分形维数均值的像素值置为白色，小于均值的置为黑色，得到二值化后的分形维数图。
17.对计算出的分形维数d计算均值，将大于均值的像素置为255（白色），小于均值的像素置为0（黑色）。分形维数越大，表征地形起伏越大，地表越粗糙；分形维数越小，表示地形起伏相对平坦，表面越光滑。
18.图4（a）展示了其后向散射总功率图，图4（b）为对计算得到的分形维数d图像进行了二值化，其中大于分形维数d平均值的置为255（白色），小于均值的像素置为0（黑色）。从图4（a）、图4（b）中可以看出，分形维数图的纹理特征基本可以反映总功率图的纹理特征，在丘陵地区等起伏较大的区域较亮，在相对平坦的地区较暗。因此，可以根据后向散射总功率计算的分形维数来进行地表平坦地区的“粗选取”。
19.步骤4、对二值化后的分形维数图，筛选出所有满足不包含白色像素的连续n
×
n像素区域，作为粗选取出的候选区域，n为像素区域宽度；具体包括：根据分形维数图，选取不包含白色像素的最大连续区域n
×
n，作为平坦地区的候选区域。具体实现方法为采用一个n
×
n的滑窗在图像上滑动，检查每次滑动后窗口内部像素值的和是否恒小于1，若满足即记录下此时滑窗的位置信息，若不满足则继续向下滑动，直至全部搜索完毕。
20.假设n=200，对图4的分形维数图进行逐像素点筛查，可获得多个满足条件的n
×
n像素候选区域，如图5所示，从图5中可以看出，搜索到的区域（白色实线矩形框）在后向散射系数图中(c)显示为雷达后向散射的均匀区域，无明显地形起伏，这表明根据后向散射系数的分形维数d进行平坦区域的粗选取是可靠的。
21.步骤5、利用yamaguchi四元素分解方法对sar回波进行极化分解，获取单次散射分量、二次散射分量、体散射分量以及螺旋散射分量的大小；散射分量的计算由yamaguchi四元素分解方法实现，该分解方法将地物目标的极化相干矩阵t分解为四种散射机制，分别为：单次散射、二次散射、体散射和螺旋散射分量，如下式所示：（2）其中，为泡利基散射矢量，上标表示共轭转置，s为目标散射系数，s
hh
、s
vv
、s
hv
的下标h和v分别表示水平和垂直极化，fs，fd，fv和fc为待确定的展开系数。单次散射由布拉格表面散射模型表征，其相干矩阵具有如下形式：（3）二次散射由二面角散射模型表征，其相干矩阵为：
（4）其中，α和β为待确定系数，|α|《1，|β|《1。
22.体散射模型的相干矩阵tv将根据的大小，选择服从不同概率密度函数的偶极子模型。随机取向的偶极子模型的相干矩阵如下所示：（5）服从正弦分布的偶极子模型的相干矩阵如下所示：（6）当-2db《p《2db时，体散射模型为t
v1
；当p《-2db时，体散射模型为t
v2
；当p》2db时，体散射模型为t
v3
。
23.对于具有旋转不变性的螺旋散射分量，其相干矩阵可以表示为：（7）对上述gf3数据进行yamaguchi四元素极化分解，j为虚数单位。
24.步骤6、针对每一个候选区域，筛选出单次散射分量占比最大的n
×
n像素区域，即为精选取出的平坦地区，其中n《n。
25.平坦地区主要以单次散射分量占主导，因此，精选取就是在粗选取出的候选区域内，根据每个像素的各散射分量功率的占比进行搜索，寻找单次散射功率占比最大的连续区域（n
×
n像素），即为最终所选择的平坦区域。
26.对候选区域中的像素按照各散射分量的占比进行n
×
n像素区域的搜索。具体实现方法为采用一个n
×
n的滑窗在图像上滑动，每滑动一次就计算滑窗内部所有像素的单次散射功率总和，直至滑动完整个候选区域。最后记录下单次散射功率占比最大的滑窗位置，即为最终所选取的平坦区域。
27.考虑到gf3全极化影像分辨率为8m
×
8m，在此我们设定n=40，最终得到如下所示选取区域（白色虚线矩形框，320m
×
320m），该区域拥有候选区域内（白色实现矩形框）最大的
单次散射功率占比，达到88.7%。如图6，最终筛选出的40
×
40像素（320m
×
320m）平坦区域（白色虚线框2）。白色实线矩形框1为根据分形维数筛选出的粗选取区域。如图7为整个sar图像场景内，全部40
×
40像素窗口的高程标准差与最大高程偏差的散点密度图。
28.为了验证所提方法的有效性，在此通过30m分辨率的srtm1dem数据分析了该区域的高程信息。首先将dem数据和sar图像进行地理编码，然后计算整个场景内，所有40
×
40像素窗口内的最大高程偏差和高程标准差，得到其散点密度分布图如图7所示。利用本发明所提方法选取出的平坦区域（白色虚线框所示区域）的最大高程偏差为3.39m，高程标准差为0.72m。因此，本发明所提方法可以在不依赖额外dem数据的情况下，有效实现自然地质表面的平坦地区选取。
29.尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，且应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。