InSAR微形变文化遗产影响定量刻画方法

摘要

本申请公开了一种InSAR微形变文化遗产影响定量刻画方法，包括步骤：预设待测区域，待测区域包括多个测量点，获取待测区域的影像数据和差分干涉图；预设第一目标，构建目标模型获取残余相位，迭代求解得到整体相干系数；或者，将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值，得到形变速率值的栅格图像，得到像素的空间分异性的标准差；或者，将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值，拟合得到测量点的时间维形变加速度，获取形变加速度场。本发明合理利用异常形变及时间序列方法所提供的客观实测数据，从形变的多个角度出发，整合出综合风险评价系数，提升了风险评估的可操作性。

说明书

技术领域

本申请涉及雷达干涉和遗产保护技术领域，尤其涉及一种InSAR微形变文化遗产影响定量刻画方法。

背景技术

遗产地，如文物古迹、自然奇观，是人类历史上具有重要意义和巨大价值的无以替代的财富。不幸的是，遗产地在自然灾害和人类活动的侵蚀下正在变得日益不稳定。自然灾害包括滑坡、地震、洪水、恶劣的天气以及气候突变。人类活动包括战争、资源的过度开采、城镇化建设和失控的旅游业。这些外界因素使得很多遗产地岌岌可危，监测并实施有效的保护势在必行。

传统的形变测量方法运用雷达干涉原理与数学统计模型等进行模型构建，虽能较为准确地反映文化遗产在时间、空间、综合的变形情况，但模型的精度易受卫星影像处理、人工数据筛选、数学模型选择等的精度影响。

发明内容

本申请公开了一种InSAR微形变文化遗产影响定量刻画方法包括步骤：

预设待测区域，所述待测区域包括多个测量点，获取所述待测区域的影像数据，通过所述影像数据获取差分干涉图；

预设第一目标，根据所述第一目标构建目标模型，在构建所述目标模型过程中获取所述残余相位，并将所述残余相位迭代求解得到整体相干系数；

或者，将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，通过插值方法得到所述形变速率值的栅格图像，将所述栅格图像划分为多个像素窗口，所述像素窗口包括多个阵列排布的像素，遍历所述栅格图像得到所述像素的空间分异性的标准差；

或者，将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，根据所述形变时间序列和所述形变速率值拟合得到所述测量点的时间维形变加速度，通过插值方法获取形变加速度场，优选的，获取综合风险评价系数，通过下列方法得到所述综合风险评价系数：

其中，σ为比例因子，所述比例因子控制空间、时间的贡献比例，std为所述空间分异性的标准差，acc为所述时间维形变加速度，|γ|为所述整体相干系数。

优选的，通过以下方法得到所述整体相干系数：

其中，N为所述干涉图数量，

优选的，整体相干系数|γ|的范围为0至1。

优选的，通过以下方法得到所述空间分异性的标准差：

其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，

优选的，通过以下方法得到所述变形速率均值：

其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，v

优选的，通过以下方法得到所述时间维形变加速度：

其中，

优选的，通过以下方法得到所述像素的加速度：

其中，y(t

优选的，根据所述形变时间序列和所述形变速率值通过最小二乘法拟合得到所述测量点的形变加速度。

优选的，0＜σ≤1。

与现有技术相比，本发明提供的InSAR微形变文化遗产影响定量刻画方法，达到如下有益效果：

本发明合理利用异常形变及时间序列方法所提供的客观实测数据，从形变的多个角度出发，如可靠性、空间的错位程度、时间维的恶化趋势等，提出了多个定量刻画指标，且创造性的融合三者，整合出综合风险评价系数，提升了风险评估的可操作性，既考虑了形变的机制与文化遗产的变化动向，又采用了较为合理的数学模型与最小二乘的数据拟合方法，展示了可衡量的评价指标，提出了较为合理的遗产风险评估数学模型，为遗产保护工作者提供了科学、理性的决策参考。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1为本发明InSAR微形变文化遗产影响定量刻画方法的一种流程图；

图2为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图；

图3为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图；

图4为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图；

图5为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。应注意到，所描述的实施例实际上仅仅是本发明一部分实施例，而不是全部的实施例，且实际上仅是说明性的，决不作为对本发明及其应用或使用的任何限制。本申请的保护范围当视所附权利要求所界定者为准。

实施例1：

参见图1至图3所示，图1为本发明InSAR微形变文化遗产影响定量刻画方法的一种流程图，图2为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图，图3为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图。本实施例提供的InSAR微形变文化遗产影响定量刻画方法，包括以下几种情况：

第一种：继续结合图1所示，InSAR微形变文化遗产影响定量刻画方法，包括步骤：

步骤S101：预设待测区域，待测区域包括多个测量点，获取待测区域的影像数据，通过影像数据获取差分干涉图。其中，影像数据指的是待测区域的SAR影像，SAR"SyntheticAperture Radar，”合成孔径雷达。

步骤S102：预设第一目标，根据第一目标构建目标模型，在构建目标模型过程中获取残余相位，并将残余相位迭代求解得到整体相干系数。

在步骤S102中，本方法中的“构建目标模型”是指通过已有的数据和参数，根据各个参数存在的数学关系，建立数学公式模型，以逼近模拟计算得到所需要的未知参数，此过程为建模。具体到本方法中，则是根据形变特点、建立形变速率、高程误差、轨道误差、大气相位等参数的方程，进行相关参数的迭代求解，在此过程中，每一期迭代，均会由逐像素的测量值和模型估计值之差产生一个残余相位，通过多次迭代求解，使得残余项减小并收敛，使得模型估计更为精确。

第二种：继续结合图2所示，InSAR微形变文化遗产影响定量刻画方法，包括步骤

步骤:S201：预设待测区域，待测区域包括多个测量点，获取待测区域的影像数据，通过影像数据获取差分干涉图。其中，影像数据指的是待测区域的SAR影像，SAR"SyntheticAperture Radar，”合成孔径雷达。

步骤S202：将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值，通过插值方法得到形变速率值的栅格图像，将栅格图像划分为多个像素窗口，像素窗口包括多个阵列排布的像素，遍历栅格图像得到像素的空间分异性的标准差。

第三种：继续结合图3所示，InSAR微形变文化遗产影响定量刻画方法，包括步骤

步骤S301：预设待测区域，待测区域包括多个测量点，获取待测区域的影像数据，通过影像数据获取差分干涉图。

其中，影像数据指的是待测区域的SAR影像，SAR"Synthetic Aperture Radar”合成孔径雷达。

步骤S302：将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值；

步骤S303：根据形变时间序列和形变速率值拟合得到测量点的时间维形变加速度，通过插值方法获取形变加速度场。

可以理解的是，上述任意一种InSAR微形变文化遗产影响定量刻画方法，均可以用于对微形变文化遗产影响定量刻画。其中，第一种InSAR微形变文化遗产影响定量刻画方法，可以利用得到的整体相干系数来定量刻画，整体相干系数值越接近于1(整体相干系数范围为0-1)，即表征未知参数建模及参数估计精度越准确、可靠。因此，整体相干系数可表征遗产目标微变测量的可靠性和置信度。第二种，InSAR微形变文化遗产影响定量刻画方法，可以利用空间分异性标准差来定量刻画，空间分异性标准差针对形变的空间角度，计算得到所有的空间分异性标准差后，即可得形变空间分异性的专题图像，根据像素间差异的大小，所得值可按照统计规律或实例经验分级渲染，数值越大则表示邻域内形变速率差异较大，空间错位性越大，所监测的遗产即形变风险越高。第三种InSAR微形变文化遗产影响定量刻画方法，可以利用得到的形变加速度场来定量刻画，可用克里金插值生成形变加速度场，分级渲染成专题图。本发明对上述InSAR微形变文化遗产影响定量刻画方法，不做限定，可以根据实际情况设置，下文不再赘述。

实施例2：

参见图4所示，图4为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图，

步骤S401，预设待测区域，所述待测区域包括多个测量点，获取所述待测区域的影像数据，通过所述影像数据获取差分干涉图；

步骤S402，预设第一目标，根据所述第一目标构建目标模型，在构建所述目标模型过程中获取所述残余相位，并将所述残余相位迭代求解得到整体相干系数；

步骤S403，将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，通过插值方法得到所述形变速率值的栅格图像，将所述栅格图像划分为多个像素窗口，所述像素窗口包括多个阵列排布的像素，遍历所述栅格图像得到所述像素的空间分异性的标准差；

步骤S404，根据所述形变时间序列和所述形变速率值拟合得到所述测量点的时间维形变加速度，通过插值方法获取形变加速度场。

步骤S405，获取综合风险评价系数，通过下列方法得到综合风险评价系数：

其中，σ为比例因子，比例因子控制空间、时间的贡献比例，std为空间分异性的标准差，acc为时间维形变加速度，|γ|为整体相干系数。

可以理解的是，目前雷达干涉技术应用于文化遗产微形变监测存在精细刻画短板。本专利运用三个单独指标以及一个整合性的综合指标，将测绘、遥感技术相结合，应用于遗产保护工作，实现了形变刻画以及风险评估定性评价到定量化描述的过渡；以及从建模形变反演可信度的角度、空间维的宏观展示以及时间维的长期时间序列分析出发，给予文化遗产微形变多方面的计算考量，且将各定量指标综合分析整合为统一的风险评价指标，并采取归一化的定义方法以作指示，更为清晰、直观；并且，整体相干系数表征雷达干涉形变反演的可信度；形变速率空间分异性表征了形变在空间维的异质性，指示非均匀形变，反映遗址裂缝、倾斜乃至倒塌的可能性；形变加速度为形变在时间维的变化模式，可揭示文化遗产病害现象正在进行甚至恶化；而综合风险评价系数通过整合三个定量指标，可综合表征文化遗产本体变形坍塌风险和文化遗产保护迫切程度，实现从文化遗产微变监测到风险评估的专题转化和成果落地。

其中，InSAR指的是干涉雷达指采用干涉测量技术的合成孔径雷达，是新近发展起来的空间对地观测技术，是传统的SAR遥感技术与射电天文干涉技术相结合的产物。它利用雷达向目标区域发射微波，然后接收目标反射的回波，得到同一目标区域成像的SAR复图像对，若复图像对之间存在相干条件，SAR复图像对共轭相乘可以得到干涉图，根据干涉图的相位值，得出两次成像中微波的路程差，从而计算出目标地区的地形、地貌以及表面的微小变化，可用于数字高程模型建立、地壳形变探测等。

实施例3：

继续结合图4所示，本实施例提供的InSAR微形变文化遗产影响定量刻画方法，包括步骤：

步骤S401：预设待测区域，待测区域包括多个测量点，获取待测区域的影像数据，通过影像数据获取差分干涉图。

步骤S402：预设第一目标，根据第一目标构建目标模型，在构建目标模型过程中获取残余相位，并将残余相位迭代求解得到整体相干系数。

在步骤S402中，可选的，通过以下方法得到整体相干系数：

其中，N为干涉图数量，w

可选的，整体相干系数|γ|的范围为0至1。其中，整体相干系数|γ|越接近于1，即表征遗产目标微变测量的可靠性和置信度越高。

步骤S403：将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值，通过插值方法得到形变速率值的栅格图像，将栅格图像划分为多个像素窗口，像素窗口包括多个阵列排布的像素，遍历栅格图像得到像素的空间分异性的标准差。在步骤S203中，可选的，通过以下方法得到空间分异性的标准差：

其中，n为像素窗口中像素的行数或者列数，i为正整数，且i≤n，

可选的，通过以下方法得到变形速率均值：

其中，n为像素窗口中像素的行数或者列数，i为正整数，且i≤n，v

进一步的，空间分异性的标准差包括第一空间分异性的标准差、第二空间分异性的标准差和第三空间分异性的标准差；

第一空间分异性的标准差小于0.2mm/年，此时所监测的文化遗产为第一安全等级，第一安全等级为监测的遗产安全。

第二空间分异性的标准差范围为0.3-0.5mm/年，此时所监测的文化遗产为第二安全等级，第二安全等级为监测的遗产比较安全。

第三空间分异性的标准差大于0.5mm/年，此时所监测的文化遗产为第三安全等级，第三安全等级为监测的遗产有一定风险，具体数值、分段的方式以及数量均可以根据具体情况设置，本发明对此不作具体限定，下文不再赘述。

其中，亦可用其他指标衡量，如范围(range)，计算一定窗口邻域内的最大最小值之差以反映空间的差异性，可同样分级渲染为空间分异性专题图。但采用不同的指标、针对不同的遗产情况，所得数值范围并非一致，即分级的标准也应进行合理调整，以较为精准地反映遗产地的形变风险等级。采用邻域统计计算的方式也会在一定程度上降低图像的分辨率，故而建议采用较小的邻域窗口，如2*2像素，邻域越小，所能反映的空间分异性越精细、越敏感，分辨率也损失的更少，但计算量较大，反之，邻域越大，分异性越迟钝，分辨率损失越大，但减少了一部分计算量。但本发明不限于此，实际情况中，可以根据遗产地情况和特点选择合适的窗口和指标，下文不再赘述。

步骤S404：根据形变时间序列和形变速率值拟合得到测量点的时间维形变加速度，通过插值方法获取形变加速度场。

在步骤S404中，可选的，通过以下方法得到时间维形变加速度：

其中，

可选的，通过以下方法得到像素的加速度：

其中，y(t

例如：若监测时间由2019年1月1日至2020年1月1日，共m期，则取起点时间为t

在步骤S404中，可选的，根据所述形变时间序列和所述形变速率值通过最小二乘法拟合得到所述测量点的形变加速度。本实施例仅示意最小二乘法为例，但是不限于此，可以根据具体情况设置计算方法，下文不再赘述。

步骤S405：获取综合风险评价系数，通过下列方法得到综合风险评价系数：

其中，σ为比例因子，比例因子控制空间、时间的贡献比例，std为空间分异性的标准差，acc为时间维形变加速度，|γ|为整体相干系数。

可选的，0＜σ≤1。

可选的，σ＝0.5。本实施例仅是以σ＝0.5为例，但本发明对比例因子σ的大小不限于此，可以根据实际情况具体设置，下文不再赘述。

其中，为比例因子的范围是0＜σ≤1，据实际需要选取，优选的σ＝0.5，其控制空间、时间的贡献比例。故而CRAC

可以理解的是，目前雷达干涉技术应用于文化遗产微形变监测存在精细刻画短板。本专利运用三个单独指标以及一个整合性的综合指标，将测绘、遥感技术相结合，应用于遗产保护工作，实现了形变刻画以及风险评估定性评价到定量化描述的过渡；以及从建模形变反演可信度的角度、空间维的宏观展示以及时间维的长期时间序列分析出发，给予文化遗产微形变多方面的计算考量，且将各定量指标综合分析整合为统一的风险评价指标，并采取归一化的定义方法以作指示，更为清晰、直观；并且，整体相干系数表征雷达干涉形变反演的可信度；形变速率空间分异性表征了形变在空间维的异质性，指示非均匀形变，反映遗址裂缝、倾斜乃至倒塌的可能性；形变加速度为形变在时间维的变化模式，可揭示文化遗产病害现象正在进行甚至恶化；而综合风险评价系数通过整合三个定量指标，可综合表征文化遗产本体变形坍塌风险和文化遗产保护迫切程度，实现从文化遗产微变监测到风险评估的专题转化和成果落地。

实施例4：

参见图5，图5为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图。本实施例为InSAR微形变文化遗产影响定量刻画方法的具体实施例：InSAR微形变文化遗产影响定量刻画方法包括步骤：

步骤S501：获取整体相干系数。

整体相干系数表征雷达干涉形变反演的可信度指标。

采用时序雷达干涉技术进行形变反演时，差分干涉图(去除地形、平地相位)上的相位贡献成分可包括：形变相位、DEM残余误差相位、大气延迟相位和噪声。在未知参数反演时，先提取高相干目标点(时序观测中雷达影像幅度与相位稳定)；然后利用其可靠相位观测，选用合适的形变和未知参量模型，解求地表形变线性速率、高程误差等信息。为了确保观测的相干目标点参数反演精度，可引入整体相干系数来表征观测点测量的可靠性：

其中，N为干涉图数量，w

步骤S502：获取形变速率空间分异性。

形变速率空间分异性表征了形变在空间维的异质性，非均匀形变是触发遗址裂缝、倾斜乃至倒塌的最主要原因。

通过一定邻域内测量点的标准差(或其他指标：如数值范围，即最大值与最小值之差)来衡量，量测点年形变速率经过插值处理后，得到相应的年形变速率场，处理转换为栅格数据后，以计算标准差为例：设像素的窗口为n*n，即搜索该像素长宽各为n像素的邻域内所有像素值，即对应形变速率值，则形变速率均值为：

则空间分异性标准差可表示为：

计算得到所有的标准差后，即可得形变空间分异性的专题图像，根据像素间差异的大小，所得值可按照统计规律或实例经验分级渲染，数值越大则表示邻域内形变速率差异较大，空间错位性越大，所监测的遗产即形变风险越高。如某专题图示例，根据其数据特点可将标准差分为了三个等级：小于0.2mm/年，等级1，安全；0.3-0.5mm/年，等级2，比较安全；大于0.5mm/年，等级3，有一定风险。

亦可用其他指标衡量，如范围(range)，计算一定窗口邻域内的最大最小值之差以反映空间的差异性。

但采用不同的指标、针对不同的遗产情况，所得数值范围并非一致，即分级的标准也应进行合理调整，以较为精准地反映遗产地的形变风险等级。采用邻域统计计算的方式也会在一定程度上降低图像的分辨率，故而建议采用较小的邻域窗口，如2*2像素，邻域越小，所能反映的空间分异性越精细、越敏感，分辨率也损失的更少，但计算量较大，反之，邻域越大，分异性越迟钝，分辨率损失越大，但减少了一部分计算量。实际情况中，需根据遗产地情况和特点选择合适的窗口和指标。

步骤S503：获取形变加速度。

形变加速度是形变在时间维的变化模式，形变的加速暗示着病害现象正在进行甚至恶化。

相对于某一时间点的形态位置，遗产地随着时间产生的形变是呈现无特定规律的震荡变化的，这就需要从数学建模的方式拟合出误差尽可能小的变化趋势，加速度是反映形变趋势的最佳指示因子，在已知拟合所得形变年速率的情况下，可同形变值时间序列，采用最小二乘法拟合出每个量测点的加速度值，原始速度、加速度公式如下：

已知

其中：y(t

例如：若监测时间由2019年1月1日至2020年1月1日，共m期，则取起点时间为t

λ为m期的纵坐标向量值，t为m期的时间横坐标向量值。由最小二乘法，将每个点的m个散点坐标由以上公式拟合出误差最小的系数

批量计算所有量测点的形变加速度后，取加速度绝对值：

可用克里金插值生成形变加速度场，分级渲染成专题图。同空间分异性一样，分级的标准需合理选择。如某示例专题图中，可根据数据统计特点分为了若干个等级，以颜色的变化指示加速度大小的不同，加速度越大的区域风险等级越高，反之越低，越稳定安全。

步骤S504：获取综合风险评价系数。

通过对以上三个定量指标的综合，从形变反演可靠性、空间维、时间维的形变情况出发，计算得到一个可综合反映遗产地形变风险程度的定量指标，可为文物保护工作者提供更为直观的遗产地保护迫切程度。由整体相干系数、形变速率空间分异性、形变加速度的数据特性和现实意义，定义CRAC如下：

σ为比例因子(范围0到1)，据实际需要选取，一般建议σ＝0.5，其控制空间、时间的贡献比例，std为表征空间分异性的标准差，acc为时间维形变加速度，|γ|为表征数据可信度的整体相干系数。对标准差及加速度作归一化处理，可避免二者量级差异带来的影响。故而CRAC

通过以上各实施例可知，本申请存在的有益效果是：

本发明合理利用异常形变及时间序列方法所提供的客观实测数据，从形变的多个角度出发，如可靠性、空间的错位程度、时间维的恶化趋势等，提出了多个定量刻画指标，且创造性的融合三者，整合出综合风险评价系数，提升了风险评估的可操作性，既考虑了形变的机制与文化遗产的变化动向，又采用了较为合理的数学模型与最小二乘的数据拟合方法，展示了可衡量的评价指标，提出了较为合理的遗产风险评估数学模型，为遗产保护工作者提供了科学、理性的决策参考。

上面通过附图和实施例，对于本申请的技术方案，虽然已经通过例子对本申请的一些特定实施例进行了详细说明，但是本领域的技术人员应该理解，以上例子仅是为了进行说明，而不是为了限制本发明的范围。尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。本发明的范围由所附权利要求来限定。