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

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

Info

Publication number
CN112949989B
CN112949989B CN202110143545.8A CN202110143545A CN112949989B CN 112949989 B CN112949989 B CN 112949989B CN 202110143545 A CN202110143545 A CN 202110143545A CN 112949989 B CN112949989 B CN 112949989B
Authority
CN
China
Prior art keywords
deformation
acceleration
micro
region
insar
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110143545.8A
Other languages
English (en)
Other versions
CN112949989A (zh
Inventor
陈富龙
刘涵薇
周伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Aerospace Information Research Institute of CAS
Original Assignee
Aerospace Information Research Institute of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Aerospace Information Research Institute of CAS filed Critical Aerospace Information Research Institute of CAS
Priority to CN202110143545.8A priority Critical patent/CN112949989B/zh
Publication of CN112949989A publication Critical patent/CN112949989A/zh
Application granted granted Critical
Publication of CN112949989B publication Critical patent/CN112949989B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0635—Risk analysis of enterprise or organisation activities
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904—SAR modes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Remote Sensing (AREA)
  • Human Resources & Organizations (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Strategic Management (AREA)
  • Theoretical Computer Science (AREA)
  • Economics (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Quality & Reliability (AREA)
  • Operations Research (AREA)
  • Marketing (AREA)
  • Game Theory and Decision Science (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

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

Description

InSAR微形变文化遗产影响定量刻画方法
技术领域
本申请涉及雷达干涉和遗产保护技术领域，尤其涉及一种InSAR微形变文化遗产影响定量刻画方法。
背景技术
遗产地，如文物古迹、自然奇观，是人类历史上具有重要意义和巨大价值的无以替代的财富。不幸的是，遗产地在自然灾害和人类活动的侵蚀下正在变得日益不稳定。自然灾害包括滑坡、地震、洪水、恶劣的天气以及气候突变。人类活动包括战争、资源的过度开采、城镇化建设和失控的旅游业。这些外界因素使得很多遗产地岌岌可危，监测并实施有效的保护势在必行。
传统的形变测量方法运用雷达干涉原理与数学统计模型等进行模型构建，虽能较为准确地反映文化遗产在时间、空间、综合的变形情况，但模型的精度易受卫星影像处理、人工数据筛选、数学模型选择等的精度影响。
发明内容
本申请公开了一种InSAR微形变文化遗产影响定量刻画方法包括步骤：
预设待测区域，所述待测区域包括多个测量点，获取所述待测区域的影像数据，通过所述影像数据获取差分干涉图；
预设第一目标，根据所述第一目标构建目标模型，在构建所述目标模型过程中获取所述残余相位，并将所述残余相位迭代求解得到整体相干系数；
或者，将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，通过插值方法得到所述形变速率值的栅格图像，将所述栅格图像划分为多个像素窗口，所述像素窗口包括多个阵列排布的像素，遍历所述栅格图像得到所述像素的空间分异性的标准差；
或者，将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，根据所述形变时间序列和所述形变速率值拟合得到所述测量点的时间维形变加速度，通过插值方法获取形变加速度场，优选的，获取综合风险评价系数，通过下列方法得到所述综合风险评价系数：
其中，σ为比例因子，所述比例因子控制空间、时间的贡献比例，std为所述空间分异性的标准差，acc为所述时间维形变加速度，|γ|为所述整体相干系数。
优选的，通过以下方法得到所述整体相干系数：
其中，N为所述干涉图数量，wx,y,i为所述残余相位，x，y代表所述像素，为所述数据建模的值。
优选的，整体相干系数|γ|的范围为0至1。
优选的，通过以下方法得到所述空间分异性的标准差：
其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，为所述形变速率的均值，vi为形变速率值。
优选的，通过以下方法得到所述变形速率均值：
其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，vi为所述形变速率值。
优选的，通过以下方法得到所述时间维形变加速度：
其中，为所述像素的加速度。
优选的，通过以下方法得到所述像素的加速度：
其中，y(ti)为形变值，ti为各个时间相对于起始时间的差值以年为单位作归一化处理的数值，t0为所述起始时间。
优选的，根据所述形变时间序列和所述形变速率值通过最小二乘法拟合得到所述测量点的形变加速度。
优选的，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为干涉图数量，wx,y,i为残余相位，残余相位为数学模型测量值与估计值之间的差，x，y代表像素，为数据建模的值。
可选的，整体相干系数|γ|的范围为0至1。其中，整体相干系数|γ|越接近于1，即表征遗产目标微变测量的可靠性和置信度越高。
步骤S403：将待测区域根据雷达干涉形变反演得到测量点的形变时间序列和形变速率值，通过插值方法得到形变速率值的栅格图像，将栅格图像划分为多个像素窗口，像素窗口包括多个阵列排布的像素，遍历栅格图像得到像素的空间分异性的标准差。在步骤S203中，可选的，通过以下方法得到空间分异性的标准差：
其中，n为像素窗口中像素的行数或者列数，i为正整数，且i≤n，为形变速率的均值，vi为形变速率值。
可选的，通过以下方法得到变形速率均值：
其中，n为像素窗口中像素的行数或者列数，i为正整数，且i≤n，vi为形变速率值。
进一步的，空间分异性的标准差包括第一空间分异性的标准差、第二空间分异性的标准差和第三空间分异性的标准差；
第一空间分异性的标准差小于0.2mm/年，此时所监测的文化遗产为第一安全等级，第一安全等级为监测的遗产安全。
第二空间分异性的标准差范围为0.3-0.5mm/年，此时所监测的文化遗产为第二安全等级，第二安全等级为监测的遗产比较安全。
第三空间分异性的标准差大于0.5mm/年，此时所监测的文化遗产为第三安全等级，第三安全等级为监测的遗产有一定风险，具体数值、分段的方式以及数量均可以根据具体情况设置，本发明对此不作具体限定，下文不再赘述。
其中，亦可用其他指标衡量，如范围(range)，计算一定窗口邻域内的最大最小值之差以反映空间的差异性，可同样分级渲染为空间分异性专题图。但采用不同的指标、针对不同的遗产情况，所得数值范围并非一致，即分级的标准也应进行合理调整，以较为精准地反映遗产地的形变风险等级。采用邻域统计计算的方式也会在一定程度上降低图像的分辨率，故而建议采用较小的邻域窗口，如2*2像素，邻域越小，所能反映的空间分异性越精细、越敏感，分辨率也损失的更少，但计算量较大，反之，邻域越大，分异性越迟钝，分辨率损失越大，但减少了一部分计算量。但本发明不限于此，实际情况中，可以根据遗产地情况和特点选择合适的窗口和指标，下文不再赘述。
步骤S404：根据形变时间序列和形变速率值拟合得到测量点的时间维形变加速度，通过插值方法获取形变加速度场。
在步骤S404中，可选的，通过以下方法得到时间维形变加速度：
其中，为像素的加速度。
可选的，通过以下方法得到像素的加速度：
其中，y(ti)为形变值，ti为各个时间相对于起始时间的差值以年为单位作归一化处理的数值，t0为起始时间。
例如：若监测时间由2019年1月1日至2020年1月1日，共m期，则取起点时间为t0＝0，各个时间相对于起始时间的差值以年为单位作归一化处理为ti则终止时间为tm＝1。
在步骤S404中，可选的，根据所述形变时间序列和所述形变速率值通过最小二乘法拟合得到所述测量点的形变加速度。本实施例仅示意最小二乘法为例，但是不限于此，可以根据具体情况设置计算方法，下文不再赘述。
步骤S405：获取综合风险评价系数，通过下列方法得到综合风险评价系数：
其中，σ为比例因子，比例因子控制空间、时间的贡献比例，std为空间分异性的标准差，acc为时间维形变加速度，|γ|为整体相干系数。
可选的，0＜σ≤1。
可选的，σ＝0.5。本实施例仅是以σ＝0.5为例，但本发明对比例因子σ的大小不限于此，可以根据实际情况具体设置，下文不再赘述。
其中，为比例因子的范围是0＜σ≤1，据实际需要选取，优选的σ＝0.5，其控制空间、时间的贡献比例。故而CRACσ,i为范围0到1的正值，其数值越大，则表示遗产地建筑坍塌风险越大，反之则较为安全，可根据此指标计算综合风险评估专题图，并作分级渲染。同上，分级标准需根据实际情况与统计特征进行合理选择。
可以理解的是，目前雷达干涉技术应用于文化遗产微形变监测存在精细刻画短板。本专利运用三个单独指标以及一个整合性的综合指标，将测绘、遥感技术相结合，应用于遗产保护工作，实现了形变刻画以及风险评估定性评价到定量化描述的过渡；以及从建模形变反演可信度的角度、空间维的宏观展示以及时间维的长期时间序列分析出发，给予文化遗产微形变多方面的计算考量，且将各定量指标综合分析整合为统一的风险评价指标，并采取归一化的定义方法以作指示，更为清晰、直观；并且，整体相干系数表征雷达干涉形变反演的可信度；形变速率空间分异性表征了形变在空间维的异质性，指示非均匀形变，反映遗址裂缝、倾斜乃至倒塌的可能性；形变加速度为形变在时间维的变化模式，可揭示文化遗产病害现象正在进行甚至恶化；而综合风险评价系数通过整合三个定量指标，可综合表征文化遗产本体变形坍塌风险和文化遗产保护迫切程度，实现从文化遗产微变监测到风险评估的专题转化和成果落地。
实施例4：
参见图5，图5为本发明InSAR微形变文化遗产影响定量刻画方法的又一种流程图。本实施例为InSAR微形变文化遗产影响定量刻画方法的具体实施例：InSAR微形变文化遗产影响定量刻画方法包括步骤：
步骤S501：获取整体相干系数。
整体相干系数表征雷达干涉形变反演的可信度指标。
采用时序雷达干涉技术进行形变反演时，差分干涉图(去除地形、平地相位)上的相位贡献成分可包括：形变相位、DEM残余误差相位、大气延迟相位和噪声。在未知参数反演时，先提取高相干目标点(时序观测中雷达影像幅度与相位稳定)；然后利用其可靠相位观测，选用合适的形变和未知参量模型，解求地表形变线性速率、高程误差等信息。为了确保观测的相干目标点参数反演精度，可引入整体相干系数来表征观测点测量的可靠性：
其中，N为干涉图数量，wx,y,i为残余相位，x,y代表像素，为建模估计值。wx,y,i整体残余相位越小，则|γ|值就越接近于1(整体相干系数范围为0-1)；即表征未知参数建模及参数估计精度越准确、可靠。因此，整体相干系数可表征遗产目标微变测量的可靠性和置信度。
步骤S502：获取形变速率空间分异性。
形变速率空间分异性表征了形变在空间维的异质性，非均匀形变是触发遗址裂缝、倾斜乃至倒塌的最主要原因。
通过一定邻域内测量点的标准差(或其他指标：如数值范围，即最大值与最小值之差)来衡量，量测点年形变速率经过插值处理后，得到相应的年形变速率场，处理转换为栅格数据后，以计算标准差为例：设像素的窗口为n*n，即搜索该像素长宽各为n像素的邻域内所有像素值，即对应形变速率值，则形变速率均值为：
则空间分异性标准差可表示为：
计算得到所有的标准差后，即可得形变空间分异性的专题图像，根据像素间差异的大小，所得值可按照统计规律或实例经验分级渲染，数值越大则表示邻域内形变速率差异较大，空间错位性越大，所监测的遗产即形变风险越高。如某专题图示例，根据其数据特点可将标准差分为了三个等级：小于0.2mm/年，等级1，安全；0.3-0.5mm/年，等级2，比较安全；大于0.5mm/年，等级3，有一定风险。
亦可用其他指标衡量，如范围(range)，计算一定窗口邻域内的最大最小值之差以反映空间的差异性。
但采用不同的指标、针对不同的遗产情况，所得数值范围并非一致，即分级的标准也应进行合理调整，以较为精准地反映遗产地的形变风险等级。采用邻域统计计算的方式也会在一定程度上降低图像的分辨率，故而建议采用较小的邻域窗口，如2*2像素，邻域越小，所能反映的空间分异性越精细、越敏感，分辨率也损失的更少，但计算量较大，反之，邻域越大，分异性越迟钝，分辨率损失越大，但减少了一部分计算量。实际情况中，需根据遗产地情况和特点选择合适的窗口和指标。
步骤S503：获取形变加速度。
形变加速度是形变在时间维的变化模式，形变的加速暗示着病害现象正在进行甚至恶化。
相对于某一时间点的形态位置，遗产地随着时间产生的形变是呈现无特定规律的震荡变化的，这就需要从数学建模的方式拟合出误差尽可能小的变化趋势，加速度是反映形变趋势的最佳指示因子，在已知拟合所得形变年速率的情况下，可同形变值时间序列，采用最小二乘法拟合出每个量测点的加速度值，原始速度、加速度公式如下：
已知的情况下，做如下变换：
其中：y(ti)为形变值(mm)，为已知的平均形变速率(mm/year)，/>为未知值，即最小二乘方法拟合出的平均形变加速率(mm/year^2)。
例如：若监测时间由2019年1月1日至2020年1月1日，共m期，则取起点时间为t0＝0，各个时间相对于起始时间的差值以年为单位作归一化处理为ti，则终止时间为tm＝1。所有处理后的时间值与形变值、已知形变速率共同计算所有监测点m期的记为λ，即问题简化为拟合公式：
λ为m期的纵坐标向量值，t为m期的时间横坐标向量值。由最小二乘法，将每个点的m个散点坐标由以上公式拟合出误差最小的系数即为所求。
批量计算所有量测点的形变加速度后，取加速度绝对值：
可用克里金插值生成形变加速度场，分级渲染成专题图。同空间分异性一样，分级的标准需合理选择。如某示例专题图中，可根据数据统计特点分为了若干个等级，以颜色的变化指示加速度大小的不同，加速度越大的区域风险等级越高，反之越低，越稳定安全。
步骤S504：获取综合风险评价系数。
通过对以上三个定量指标的综合，从形变反演可靠性、空间维、时间维的形变情况出发，计算得到一个可综合反映遗产地形变风险程度的定量指标，可为文物保护工作者提供更为直观的遗产地保护迫切程度。由整体相干系数、形变速率空间分异性、形变加速度的数据特性和现实意义，定义CRAC如下：
σ为比例因子(范围0到1)，据实际需要选取，一般建议σ＝0.5，其控制空间、时间的贡献比例，std为表征空间分异性的标准差，acc为时间维形变加速度，|γ|为表征数据可信度的整体相干系数。对标准差及加速度作归一化处理，可避免二者量级差异带来的影响。故而CRACσ,i为范围0到1的正值，其数值越大，则表示遗产地建筑坍塌风险越大，反之则较为安全，可根据此指标计算综合风险评估专题图，并作分级渲染。同上，分级标准需根据实际情况与统计特征进行合理选择。
通过以上各实施例可知，本申请存在的有益效果是：
本发明合理利用异常形变及时间序列方法所提供的客观实测数据，从形变的多个角度出发，如可靠性、空间的错位程度、时间维的恶化趋势等，提出了多个定量刻画指标，且创造性的融合三者，整合出综合风险评价系数，提升了风险评估的可操作性，既考虑了形变的机制与文化遗产的变化动向，又采用了较为合理的数学模型与最小二乘的数据拟合方法，展示了可衡量的评价指标，提出了较为合理的遗产风险评估数学模型，为遗产保护工作者提供了科学、理性的决策参考。
上面通过附图和实施例，对于本申请的技术方案，虽然已经通过例子对本申请的一些特定实施例进行了详细说明，但是本领域的技术人员应该理解，以上例子仅是为了进行说明，而不是为了限制本发明的范围。尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。本发明的范围由所附权利要求来限定。

Claims (5)

1.一种InSAR微形变文化遗产影响定量刻画方法，其特征在于，包括步骤：
预设待测区域，所述待测区域包括多个测量点，获取所述待测区域的影像数据，通过所述影像数据获取差分干涉图；
预设第一目标，根据所述第一目标构建目标模型，在构建所述目标模型过程中获取残余相位，并将所述残余相位迭代求解得到整体相干系数；
通过以下方法得到所述整体相干系数：
其中，N为所述干涉图数量，wx,y,i为所述残余相位，x，y代表像素，为所述数据建模的值；
将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，通过插值方法得到所述形变速率值的栅格图像，将所述栅格图像划分为多个像素窗口，所述像素窗口包括多个阵列排布的像素，遍历所述栅格图像得到所述像素的空间分异性的标准差；
通过以下方法得到所述空间分异性的标准差：
其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，为所述形变速率的均值，vi为形变速率值；
通过以下方法得到所述形变速率的均值：
其中，n为所述像素窗口中所述像素的行数或者列数，i为正整数，且i≤n，vi为所述形变速率值；
将所述待测区域根据雷达干涉形变反演得到所述测量点的形变时间序列和形变速率值，根据所述形变时间序列和所述形变速率值拟合得到所述测量点的时间维形变加速度，通过插值方法获取形变加速度场；
通过以下方法得到所述时间维形变加速度；
其中，为所述像素的加速度；
获取综合风险评价系数，通过下列方法得到所述综合风险评价系数：
其中，σ为比例因子，所述比例因子控制空间、时间的贡献比例，std为所述空间分异性的标准差，acc为所述时间维形变加速度，|γ|为所述整体相干系数。
2.根据权利要求1所述的InSAR微形变文化遗产影响定量刻画方法，其特征在于，整体相干系数|γ|的范围为0至1。
3.根据权利要求1所述的InSAR微形变文化遗产影响定量刻画方法，其特征在于，通过以下方法得到所述像素的加速度：
其中，y(ti)为形变值，ti为各个时间相对于起始时间的差值以年为单位作归一化处理的数值，t0为所述起始时间。
4.根据权利要求1所述的InSAR微形变文化遗产影响定量刻画方法，其特征在于，根据所述形变时间序列和所述形变速率值通过最小二乘法拟合得到所述测量点的形变加速度。
5.根据权利要求1所述的InSAR微形变文化遗产影响定量刻画方法，其特征在于，0＜σ≤1。
CN202110143545.8A 2021-02-02 2021-02-02 InSAR微形变文化遗产影响定量刻画方法 Active CN112949989B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110143545.8A CN112949989B (zh) 2021-02-02 2021-02-02 InSAR微形变文化遗产影响定量刻画方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110143545.8A CN112949989B (zh) 2021-02-02 2021-02-02 InSAR微形变文化遗产影响定量刻画方法

Publications (2)

Publication Number Publication Date
CN112949989A CN112949989A (zh) 2021-06-11
CN112949989B true CN112949989B (zh) 2024-02-06

Family

ID=76241556

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110143545.8A Active CN112949989B (zh) 2021-02-02 2021-02-02 InSAR微形变文化遗产影响定量刻画方法

Country Status (1)

Country Link
CN (1) CN112949989B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114092837B (zh) * 2021-11-05 2022-08-26 中国科学院空天信息创新研究院 一种基于长时间尺度的遗址环境遥感监测方法及系统
CN117746239B (zh) * 2023-12-19 2024-08-06 中国科学院空天信息创新研究院 一种环境要素提取与预评估方法、系统、设备和存储介质

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2005010556A1 (en) * 2003-07-18 2005-02-03 University Of Nottingham Radar position and movement measurement for geophysical monitoring
CN102608584A (zh) * 2012-03-19 2012-07-25 中国测绘科学研究院 基于多项式反演模型的时间序列InSAR形变监测方法及装置
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法
CN104123464A (zh) * 2014-07-23 2014-10-29 中国国土资源航空物探遥感中心 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN104133996A (zh) * 2014-07-25 2014-11-05 首都师范大学 一种基于云模型和数据场的地面沉降风险等级评估方法
CN106772342A (zh) * 2017-01-11 2017-05-31 西南石油大学 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN106772377A (zh) * 2017-01-18 2017-05-31 深圳市路桥建设集团有限公司 一种基于InSAR的建筑物变形监测方法
CN109085588A (zh) * 2018-09-26 2018-12-25 云南电网有限责任公司电力科学研究院 基于Terra SAR-X高分辨率聚束模式数据电网铁塔倾斜的监测方法
CN110673145A (zh) * 2019-10-24 2020-01-10 中国地质大学（北京） 一种基于间断相干的InSAR地表形变监测方法及系统

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2005010556A1 (en) * 2003-07-18 2005-02-03 University Of Nottingham Radar position and movement measurement for geophysical monitoring
CN102608584A (zh) * 2012-03-19 2012-07-25 中国测绘科学研究院 基于多项式反演模型的时间序列InSAR形变监测方法及装置
CN104111456A (zh) * 2014-07-23 2014-10-22 中国国土资源航空物探遥感中心 一种高速铁路沿线地表形变高分辨率InSAR监测方法
CN104123464A (zh) * 2014-07-23 2014-10-29 中国国土资源航空物探遥感中心 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
CN104133996A (zh) * 2014-07-25 2014-11-05 首都师范大学 一种基于云模型和数据场的地面沉降风险等级评估方法
CN106772342A (zh) * 2017-01-11 2017-05-31 西南石油大学 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN106772377A (zh) * 2017-01-18 2017-05-31 深圳市路桥建设集团有限公司 一种基于InSAR的建筑物变形监测方法
CN109085588A (zh) * 2018-09-26 2018-12-25 云南电网有限责任公司电力科学研究院 基于Terra SAR-X高分辨率聚束模式数据电网铁塔倾斜的监测方法
CN110673145A (zh) * 2019-10-24 2020-01-10 中国地质大学（北京） 一种基于间断相干的InSAR地表形变监测方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
地表形变监测的改进相干目标法;王明洲;李陶;江利明;徐侃;吴文豪;;测绘学报(第01期);36-43 *
基于相干点目标的多基线D-InSAR技术与地表形变监测;葛大庆;王艳;郭小方;刘圣伟;范景辉;;遥感学报(第04期);574-580 *

Also Published As

Publication number Publication date
CN112949989A (zh) 2021-06-11

Similar Documents

Publication Publication Date Title
Woodget et al. Subaerial gravel size measurement using topographic data derived from a UAV‐SfM approach
CN112949989B (zh) InSAR微形变文化遗产影响定量刻画方法
CN113108764A (zh) 一种溃坝过程安全监测、预警与影响评估方法
CN102144174A (zh) Sar图像序列中的永久散射体的识别和分析
Davidson et al. Joint statistical properties of RMS height and correlation length derived from multisite 1-m roughness measurements
WO2019126972A1 (zh) 基于InSAR的形变信息提取方法、终端及存储介质
Yazici et al. Investigating persistent scatterer InSAR (PSInSAR) technique efficiency for landslides mapping: a case study in Artvin dam area, in Turkey
CN111650570B (zh) 一种地基干涉雷达三维大气校正方法及系统
CN112444188B (zh) 一种多视角InSAR海堤高精度三维形变测量方法
CN112556660A (zh) 基于卫星测高数据的海域重力异常反演方法及系统
CN113281749A (zh) 一种顾及同质性的时序InSAR高相干点选取方法
Yu et al. Structural state estimation of earthquake-damaged building structures by using UAV photogrammetry and point cloud segmentation
CN116299455A (zh) 基于PSInSAR与SqueeSAR的设施形变分析方法
CN111008355A (zh) 一种基于信任传播的气象地面要素插值方法
CN104298844B (zh) 获取点阵法测量光学遥感载荷在轨mtf测量精度的方法
Brook et al. A variational interpolation method for gridding weather radar data
CN114114257A (zh) 一种坝区形变与水位相关性检测方法和装置
Feng et al. Automatic selection of permanent scatterers-based GCPs for refinement and reflattening in InSAR DEM generation
CN110310370B (zh) 一种gps与srtm点面融合的方法
Refice et al. On the use of anisotropic covariance models in estimating atmospheric DInSAR contributions
Kim et al. The curvature interpolation method for surface reconstruction for geospatial point cloud data
CN113625241A (zh) 差异沉降监控预警方法
Kupriyanov et al. Accuracy Analysis of Remote Measurement of Thermokarst Lakes Parameters for Field Dynamics Modeling Problems
Babaee et al. Assessment of noise in InSAR timeseries using least squares variance component estimation
Luo et al. Measurement of large, discontinuous displacement from digital images

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant