本发明涉及矿区地形探测技术领域，特别涉及一种基于多级匹配窗口约束下的矿区大梯度形变提取方法。
背景技术
合成孔径雷达差分干涉测量技术(d-insar，differentialinterferometricsyntheticapertureradar)已被国际上认为一种先进的遥感技术，能够为矿区提供大范围空间连续覆盖的地表形变信息，这不仅有助于定位和确定由开采沉陷引起的形变范围，而且有助于描述其开采沉陷机理。然而，受矿区地理要素时间失相干严重影响，去相干噪声的存在则难以反映其实际形变特征；其次，矿区开采沉降中心区域通常表现为沉陷速度快、大梯度且非线性沉降特征，这可能使得相邻像素之间缠绕相位差的绝对值超出其最大可监测到的形变梯度，利用现有的相位解缠算法很难恢复到真实形变梯度的正确解。
为了解决d-insar技术受矿区失相干严重和大梯度开采沉陷等影响，像素跟踪方法(pixeltrackingmethod，又称偏移跟踪方法)，由于其提取到的距离向(range)形变包含与差分干涉测量的沿雷达视线(los)方向近似相同的形变信息，通过利用sar影像强度信息可以广泛地用于矿区像元级大梯度形变提取的研究，这已被证明是传统d-insar技术在受低相干和大梯度开采沉陷监测失败情况下进行像元级大梯度形变提取的最佳替代方法。一般来说，基于归一化互相关(ncc)算法下的像素跟踪方法估算出形变信息精度高度依赖于所选择的匹配窗口大小。窗口足够小，可以提取高质量的形变信息，而误匹配的几率会由此显著增加；相对应于足够大的窗口，散斑噪声得到抑制，进一步提高了像素偏移位置匹配的置信度，但形变监测的准确性和可靠性会由此显著降低。在很大程度上，由于不同匹配窗口大小内信号方差(包括亮度或对比度)和噪声水平的显著差异性和不确定性，单凭研究者根据矿区开采沉陷规模、sar影像分辨率大小以及噪声水平，采用单一匹配窗口大小的像素跟踪方法很难寻求到理想(最优)大梯度形变估计。此外，针对上述问题，有研究着重提出关于自适应匹配窗口选择的方法，然而，这通常会造成在匹配窗口选择、计算效率方面产生冗余，且对计算机处理性能会提出更高要求，很难在实际中得以推广应用。
技术实现要素：
本发明旨在至少在一定程度上解决上述技术中存在的关键问题之一。为此，本发明的目的在于提出一种基于多级匹配窗口约束下的矿区大梯度形变提取方法，能够提高矿区大梯度形变提取过程中的计算效率，并能够提高所提取的形变信息的质量，具有简单易行、可靠性高等优点。
为达到上述目的，本发明提出了一种基于多级匹配窗口约束下的矿区大梯度形变提取方法，包括以下步骤：获取覆盖矿区开采区域大梯度形变发生前、后的两景sar影像，并对参考影像和搜索影像进行配准处理；采用多级不同大小的匹配窗口分别提取到所述矿区开采区域的初始大梯度形变信息；去除所述初始大梯度形变信息中与开采形变无关的偏移分量，以分别得到不同大小的匹配窗口下与开采形变相关的大梯度形变信息；根据所述大梯度形变信息计算同一参考像素在不同大小的匹配窗口下的像素偏移量估计的均方根误差，并根据所述均方根误差获取在不同大小的匹配窗口下的高质量大梯度形变信息；对获取不同大小的匹配窗口下的高质量大梯度形变信息进行垂向投影变换和地理编码处理，以得到地理坐标系下所述矿区开采区域的大梯度形变场空间分布信息。由此，能够提高矿区大梯度形变提取过程中的计算效率，并能够提高所提取的形变信息的质量，具有简单易行、可靠性高等优点。
另外，根据本发明上述实施例提出的基于多级匹配窗口约束下的矿区大梯度形变提取方法还可以具有如下附加的技术特征：
根据本发明的一个实施例，所述根据均方根误差获取在不同大小的匹配窗口下的高质量大梯度形变信息，具体包括：判断所述均方根误差是否小于预设阈值；选取均方根误差小于所述预设阈值的参考像素，且认为该参考像素在不同大小的匹配窗口下得到偏移估计具有高质量匹配特征，并取具有高质量匹配的参考像素的像素偏移估计均值作为高质量大梯度形变估计。
进一步地，所述根据所述均方根误差获取在不同大小的匹配窗口下的高质量大梯度形变信息，还包括：选取均方根误差不小于所述预设阈值的参考像素作为误匹配的参考像素；对于所述误匹配的参考像素，在不同大小的匹配窗口下依据snr寻找最优偏移量匹配，以获取高质量大梯度形变估计。
进一步地，在不同大小的匹配窗口下依据snr寻找最优偏移量匹配，以获取高质量大梯度形变估计，具体包括：获取具有误匹配偏移估计的参考像素在不同大小的匹配窗口下的snr；根据所述snr定义优化系数，并选取优化系数极大值对应的匹配窗口下的像素偏移量估计作为高质量大梯度形变估计。
或者，在不同大小的匹配窗口下依据snr寻找最优偏移量匹配，以获取高质量大梯度形变估计，具体包括：获取具有误匹配偏移估计的参考像素在不同大小的匹配窗口下的snr；根据所述snr对所述误匹配参考像素的像素偏移量估计进行修正，并取误匹配参考像素在修正后的像素偏移量估计均值作为高质量大梯度形变估计。
根据本发明的一个实施例，所述采用多级不同大小的匹配窗口分别提取到所述矿区开采区域的初始大梯度形变信息，具体包括：采用多级不同大小的匹配窗口，以常规像素追踪法分别获取方位向和距离向的二维亚像素级形变信息。
根据本发明的一个实施例，所述去除所述初始大梯度形变信息中与开采形变无关的偏移分量，具体包括：对获取的初始大梯度形变信息选取非矿区开采区域且具有较高置信度像素匹配的稳定像素点偏移估计，采用基于最小二乘意义上的二次多项式拟合函数进行全局偏移量估计，以消除因星载sar轨道及姿态变化造成sar影像成像几何变化而引起的全局性偏移；采用基于sar影像对基线信息和外部dem数据分别对距离向和方位向进行地形偏移量改正，以消除与地形起伏有关的局部性偏移。
附图说明
图1为根据本发明实施例的基于多级匹配窗口约束下的矿区大梯度形变提取方法的流程图；
图2为根据本发明一个具体实施例的基于多级匹配窗口约束下的矿区大梯度形变提取方法的流程图。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
下面结合附图来描述本发明实施例的基于多级匹配窗口约束下的矿区大梯度形变提取方法。
如图1所示，本发明实施例的基于多级匹配窗口约束下的矿区大梯度形变提取方法，包括以下步骤：
s1，获取覆盖矿区开采区域大梯度形变发生前、后的两景sar影像，并对参考影像和搜索影像进行配准处理。
可通过星载sar获取矿区开采区域的影像，并根据开采沉陷区域选择合适的sar影像，在选择sar影像时要充分考虑sar传感器的特性(这里主要指波长、空间分辨率、卫星轨道和重复周期)、成像质量、时间和空间基线以及地面特征等。在本发明的一个实施例中，选取位于浅层部开采且导致大梯度沉降的在开采矿区，采用覆盖研究矿区的两景1.33×2.06m分辨率的c波段radarsat-2卫星影像，范围可为6400×5760像元。
参考影像和搜索影像分别为矿区开采形变前、后的sar影像，不同于传统基于sar影像强度信息的像素偏移跟踪算法仅依据配准多项式进行像素级重采样，在本发明的一个实施例中，顾及了由地形起伏变化引起的偏移误差，即利用基于参考影像和与参考影像成像几何相一致的30m分辨率srtmdem数据确定出的配准查找表来连接两景sar影像成像几何，实现参考影像和搜索影像的像素级粗配准及精配准，从而在辅影像重采样到主影像过程中尽可能地消除由地形起伏变化引起的配准误差。
s2，采用多个不同大小的匹配窗口矿区开采区域的初始大梯度形变信息。
具体地，可采用多级不同大小的匹配窗口，对配准得到的sar影像对以常规像素追踪(pixeltracking)法分别获取矿区开采区域在方位向和距离向的二维亚像素级形变信息。此外，在利用不同匹配窗口提取初始大梯度形变信息时，可对同一参考像素位置开始进行像素偏移量估计，从而建立其相同的像素偏移参考几何，有助于后续得到同一参考像素在不同匹配窗口下的最优匹配像素位置的确定。
其中，不同大小的多级匹配窗口可依据矿区开采沉陷规模、sar影像分辨率大小以及噪声水平进行选择。在本发明的一个具体实施例中，可采用大小为32×32、40×40、48×48、56×56、64×64的匹配窗口获取在不同匹配窗口下的形变信息，其中，64×64为理想匹配窗口，即利用该窗口得到的误匹配噪声数目最少。
s3，去除初始大梯度形变信息中与开采形变无关的偏移分量，以得到不同大小的匹配窗口下与开采形变相关的大梯度形变信息。
在本发明的一个实施例中，可对获取的初始大梯度形变信息选取非矿区开采区域且具有较高置信度像素匹配的稳定像素点偏移估计，采用基于最小二乘意义上的二次多项式拟合函数进行全局偏移量估计，以消除因星载sar轨道及姿态变化造成sar影像成像几何变化而引起的全局性偏移。并且，可采用基于sar影像对基线信息和外部dem数据分别对距离向和方位向进行地形偏移量改正，以消除与地形起伏有关的局部性偏移。
具体地，初始大梯度形变信息中的初始偏移量可由以下若干分量组成：
ototal＝odisp+oorbit+oatti+otopo+onoise
其中，ototal、odisp、oorbit、oatti、otopo、onoise分别为初始偏移量、与开采形变相关的像素偏移量估计、因星载sar轨道位置引起的全局性偏移、因星载sar姿态变化引起的全局性偏移、与地形起伏变化有关的局部性偏移、局部性误匹配偏移估计。
首先，可对获取的初始大梯度形变信息选取非矿区开采区域且具有较高置信度像素匹配的稳定像素点偏移估计，采用基于最小二乘意义上的二次多项式拟合函数oglobal＝a0+a1x+a2y+a3xy+a4x2+a5y2进行全局偏移量估计，消除因星载sar轨道及姿态变化造成sar影像成像几何变化而引起的全局性偏移oglobal；其次，考虑到更高分辨率且长基线相结合的sar影像对，地形起伏变化可能造成超过0.5像素偏移，针对与地形起伏有关的局部性偏移可以进一步采用sar影像对基线信息和外部dem数据利用公式
s4，根据大梯度形变信息计算同一参考像素在不同大小的匹配窗口下的像素偏移量估计均方根误差，并根据均方根误差获取在不同大小的匹配窗口下的高质量大梯度形变信息。
在得到不同大小的匹配窗口下与开采形变相关的像素偏移量估计后，可获取同一参考像素在不同大小的匹配窗口下的像素偏移量估计的离散程度，即计算同一参考像素在不同大小的匹配窗口下的像素偏移量估计的均方根误差
在求得均方根误差后，可判断每个参考像素对应的偏移估计均方根误差是否小于预设阈值。在本发明的一个具体实施例中，预设阈值δothreshold可取0.5m。
在作出上述判断后，可选取均方根误差小于预设阈值的参考像素且认为该像素在不同大小的匹配窗口下得到偏移估计具有高质量匹配特征，并取高质量匹配的参考像素在不同匹配窗口下的像素偏移量估计均值
与此同时，可筛选出均方根误差不小于预设阈值的参考像素，并认为其在不同大小的匹配窗口下存有误匹配偏移估计的可能；对于存在误匹配偏移估计的参考像素，可在不同大小的匹配窗口下依据snr寻找最优偏移量匹配，以获取高质量大梯度形变估计。具体地，可获取误匹配的参考像素在不同大小的匹配窗口下的snr，并根据snr定义优化系数，选取优化系数极大值对应的匹配窗口下的像素偏移量估计作为高质量大梯度形变估计；或者，可获取误匹配的参考像素在不同大小的匹配窗口下的snr，并根据snr对误匹配的参考像素的像素偏移量估计进行修正，以及取修正后误匹配参考像素的像素偏移量估计均值作为高质量大梯度形变估计。
更具体地，在本发明的一个实施例中，对于误匹配的参考像素，为了计算出在不同匹配窗口下最大归一化互相关值在其局部互相关表面(这里为8×8窗口)中的峰值水平，可定义snr值：
针对较小匹配窗口存在误匹配偏移估计的情况下，获取高质量大梯度形变估计有如下两种实现方式：
方式一：获取含有误匹配偏移量估计的像素点在不同匹配窗口下的snr值：{snr1,snr2,...,snrn}，n为匹配窗口的个数，采取逐级匹配窗口下的分析比较方式，通过定义优化系数
方式二：选取snr极大值对应的匹配窗口，并找到最大互相关值对应于匹配窗口中的所在位置(xl,yl)，然后用于引导含有误匹配噪声的小匹配窗口内进行再匹配计算，使其限定在一个较小范围内(xl+δx,yl+δy)(通常为8×8窗口)寻找最大互相关值，从而减弱小匹配窗口内出现误匹配偏移量估计的可能。然后，利用多级匹配窗口约束进行匹配位置引导，依次完成小窗口内所有出现误匹配偏移量估计的修正。最后，与上述均方根误差小于预设阈值的参考像素类似，取修正后参考像素的像素偏移量估计均值
s5，对获取的不同大小的匹配窗口下的高质量大梯度形变信息进行垂向投影变换和地理编码处理，以得到地理坐标系下矿区开采区域的大梯度形变场空间分布信息。
在本发明的一个具体实施例中，如图2所示，基于多级匹配窗口约束下的矿区大梯度形变提取方法包括以下步骤：
s101，获取参考影像。即获取矿区大梯度形变发生前的sar影像作为参考影像。
s102，获取搜索影像。即获取矿区大梯度形变发生后的sar影像作为搜索影像。
s103，获取外部dem数据。
s104，进行配准处理。具体可基于外部dem及sar成像几何实现影像对精确配准。
s105，采用常规像素跟踪方法提取到在匹配窗口w1～wn下的初始大梯度形变信息ot1～otn。步骤s105所包括的子步骤如图2所示。
s106，进行全局偏移量的估计与去除。即对初始大梯度形变信息进行全局偏移量估计与去除，消除因星载sar轨道及姿态变化造成sar影像成像几何变化而引起的全局性偏移。
s107，进行残余地形偏移量的估计与去除。
s108，在保留误匹配噪声的前提下，得到每个匹配窗口下的像素偏移量估计{o(1)，o(2)，…，o(n-1)，o(n)}。该步骤后执行步骤s109和s110。
s109，判断是否有rmse＜δo。δo为预设阈值。如果是，则执行步骤s113；如果否，则执行步骤s111。
s110，得到每个匹配窗口下的snr值{snr(1)，snr(2)，…，snr(n-1)，snr(n)}。该步骤后执行步骤s111。
s111，进行最佳偏移量估计。
s112，进行噪声剔除和再匹配。即对误匹配的参考像素的像素偏移量估计进行修正。该步骤后返回步骤s108。
s113，得到高质量大梯度形变估计。
s114，对高质量大梯度形变估计进行垂向投影变换。
s115，对高质量大梯度形变估计进行地理编码。
综上所述，根据本发明实施例的基于多级匹配窗口约束下的矿区大梯度形变提取方法，通过获取覆盖矿区开采区域大梯度形变发生前、后的两景sar影像，并对参考影像和搜索影像进行配准处理；采用多级不同大小的匹配窗口分别提取到所述矿区开采区域的初始大梯度形变信息；去除所述初始大梯度形变信息中与开采形变无关的偏移分量，以分别得到在不同大小的匹配窗口下与开采形变相关的大梯度形变信息；根据所述大梯度形变信息计算同一参考像素在不同大小的匹配窗口下的像素偏移量估计的均方根误差，并根据所述均方根误差获取在不同大小的匹配窗口下的高质量大梯度形变信息；对获取不同大小的匹配窗口下的高质量大梯度形变信息进行垂向投影变换和地理编码处理，以得到地理坐标系下所述矿区开采区域的大梯度形变场空间分布信息。由此，能够提高矿区大梯度形变提取过程中的计算效率，并能够提高所提取的形变信息的质量，具有简单易行、可靠性高等优点。
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”、“顺时针”、“逆时针”、“轴向”、“径向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
在本发明中，除非另有明确的规定和限定，第一特征在第二特征“上”或“下”可以是第一和第二特征直接接触，或第一和第二特征通过中间媒介间接接触。而且，第一特征在第二特征“之上”、“上方”和“上面”可是第一特征在第二特征正上方或斜上方，或仅仅表示第一特征水平高度高于第二特征。第一特征在第二特征“之下”、“下方”和“下面”可以是第一特征在第二特征正下方或斜下方，或仅仅表示第一特征水平高度小于第二特征。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。