南京明城墙序贯PSInSAR形变估计

南京明城墙序贯PSInSAR形变估计 Sequential PSInSAR approach for the deformation monitoring of the Nanjing Ming Dynasty City Wall 1 2 first-author 王成，1998年生，男，硕士研究生，研究方向为InSAR文化遗产微变监测与评估。E-mail：wangcheng191@mails.ucas.ac.cn 王成，1998年生，男，硕士研究生，研究方向为InSAR文化遗产微变监测与评估。E-mail：wangcheng191@mails.ucas.ac.cn wangcheng191@mails.ucas.ac.cn 1 3 corresp 陈富龙，1980年生，男，研究员，研究方向为雷达遥感文化遗产智能感知与保护。E-mail：chenfl@aircas.ac.cn 陈富龙，1980年生，男，研究员，研究方向为雷达遥感文化遗产智能感知与保护。E-mail：chenfl@aircas.ac.cn chenfl@aircas.ac.cn 1 3 4 4 中国科学院空天信息创新研究院, 北京 100094 中国科学院空天信息创新研究院, 北京 100094 Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China 中国科学院大学, 北京 100049 中国科学院大学, 北京 100049 University of Chinese Academy of Sciences, Beijing 100049, China University of Chinese Academy of Sciences, Beijing 100049, China 联合国教科文组织国际自然与文化遗产空间技术中心, 北京 100094 联合国教科文组织国际自然与文化遗产空间技术中心, 北京 100094 International Centre on Space Technologies for Natural and Cultural Heritage under the Auspices of UNESC, Beijing 100094, China International Centre on Space Technologies for Natural and Cultural Heritage under the Auspices of UNESC, Beijing 100094, China 浙江省测绘科学技术研究院, 杭州 311100 浙江省测绘科学技术研究院, 杭州 311100 Zhejiang Institute of Surveying and Mapping Science and Technology, Hangzhou 311100, China Zhejiang Institute of Surveying and Mapping Science and Technology, Hangzhou 311100, China 针对永久散射体PSInSAR（Persistent Scatterer SAR Interferometry）算法流程的局限性，本文立足差分雷达干涉预处理和时序微变反演性能提升，提出了一种基于优化解空间的PSInSAR序贯处理算法。该算法能够实现新SAR影像无缝嵌入，避免时序PSInSAR整体数据处理流程溯源，实现微变参数的准实时更新。研究以南京明城墙文化遗产及其200 m缓冲区赋存环境为示范，基于2015-01—2018-02的32景降轨Cosmo-SkyMed条带影像，开展传统PSInSAR方法与优化解空间PSInSAR序贯方法形变反演与性能对比研究。结果表明，基于解空间优化的PSInSAR序贯方法简化了差分干涉流程，通过解空间搜索机制和结构改进，降低了算法时空复杂度，实现了未知参数求解约10倍效率提升。通过形变场精度交叉验证，发现两种方法形变估计结果吻合（总体误差在0—1 mm），证实了PSInSAR序贯方法的有效性与可靠性，并揭示其在遥感大数据时代文化遗产高精度、准实时微变监测中的应用潜力。 Human society has entered the big data era given the exponential growth of remote sensing data due to the emergence of higher resolution, frequent revisits, and multi-platform image acquisitions. This phenomenon raised challenges for Interferometric SAR (InSAR) and Multitemporal InSAR (MTInSAR) data processing in near real time. For instance, the traditional PSInSAR algorithm can no longer satisfy a fast response monitoring due to delay in deformation time series updating. To address the aforementioned technical limitations, we proposed a PSInSAR sequential processing algorithm characterized by the optimized searching-space to achieve the performance improvement of Differential InSAR (DInSAR) data preprocessing and MTInSAR parameter estimation. In this approach, new SAR acquisitions were seamlessly integrated into the reconstructed spatiotemporal baselines of interferograms. Then, a triple-level Delaunay network was established using the temporal coherence value (high, low, and decorrelated) on network edges. On the basis of the value inheritance from previous PSInSAR estimations, the unknown parameter estimation on PS candidates for the sequential PSInSAR was accelerated owing to the proposed searching strategy adopted to the triple-level coherence network edges. That is, the solution space was first sampled using a large searching step (for example, 10 times the measurement accuracy of unknown parameters) to determine the potential interval of the optimal solution. Then, the choice of network edge was determined on the basis of the maximum value of the temporal coherence, followed by a dedicated fine searching (with the step equivalent to the predetermined accuracy of unknown parameters) concentrating on the potential interval for the optimal inversion. Owing to the applied global-local searching strategy, the optimization of calculation efficiency and estimation accuracy can be achieved. We conducted a comparative investigation for the deformation estimation and performance assessment between the current PSInSAR and the proposed sequential PSInSAR methods using 32 scenes Cosmo SkyMed Stripmap images (in descending orbits and acquired from January 2015 to February 2018) covering the Nanjing Ming Dynasty City wall. Results indicate that a high efficiency of unknown parameter estimation (height, deformation, and thermal dilation) was obtained using the sequential PSInSAR with the adopted optimized searching-space approach, with the computation acceleration with approximately an order of 10 times. The cross comparison of the deformation velocity rates from both approaches reveals a consistent estimation as presented by the overall dispersion values ranging from 0 to 1 mm/a, which validates the feasibility and reliability of sequential PSInSAR in the deformation estimation. The driving force of detected deformation anomalies along three sections of the city wall was further exploited, providing new insights for the sustainable conservation of the heritage properties. This study implies the potential of the sequential PSInSAR method in the accurate, near real-time deformation monitoring, and preventive conservation of large-scale cultural heritage sites (i.e., Nanjing Dynasty City Wall), particularly on the emergence of big data. 遥感 PSInSAR 优化解空间 序贯 文化遗产 remote sensing PSInSAR sequential processing optimized searching-space cultural heritage 1　引言 1 文化遗产记录着历史的变迁，传承着文明的生命韵律，焕发着民族的光华，滋养着国家和民族的精神沃土，具备很高的历史价值、艺术价值、文化价值和科学价值，对其进行预防性监测与保护是时代赋予的神圣使命，也是推动文化探源和遗产地可持续发展，实现美丽中国梦的重要一环（ 陈富龙，2015 陈富龙，2015 ）。然而由于自然和人类活动的共同影响，许多遗产正经历着不可修复的破坏过程。长时间缓慢地表形变作为典型的不均匀沉降和水平位移病害，已然成为影响其景观稳定性的重要因素；因此亟需以微形变为典型定量指标，对其形变场及时空演化模式进行一盘棋、系统把控。 目前传统形变测量方法如水准测量、全球定位系统GPS（Global Positioning System）和位移传感器等实施成本高，效率低下，且无法实现科学有效的长期观测（ Chen等，2020 Chen等，2020 ）。差分雷达干涉DInSAR（Differential SAR Interferometry）因其特有的宏观性、高效性、无损性（ 邓飚 等，2010 邓飚 等，2010 ）等优势，通过提取SAR复数据中干涉相位信息来获取地表三维信息和变化信息，被广泛应用到地表形变监测研究中；然而因受到时间、空间去相干和对流层/电离层大气扰动等影响，该方法较难获得毫米级时间序列形变信息（ Gabriel等，1989 Gabriel等，1989 ）。永久散射体雷达干涉PSInSAR（Persistent Scatterer SAR Interferometry）能有效克服上述局限性，其使用观测时间段内所获取的SAR影像集，通过将影像按照一定的组合模式进行配对形成干涉图，以多时相SAR影像上受各种去相干现象影响较小的稳定点为载体，对干涉相位进行信号分离和提取，能够获取目标地物毫米级形变场与时序信息（ Ferretti等，2001 Ferretti等，2001 ）。在地表大尺度形变监测中，该方法较为成熟，已成功应用于火山活动、地震灾害、城市地表沉降等所引起的地表形变等。 在传统PSInSAR整体算法流程中，影像集差分干涉预处理和模型参数解算（包括线性速度、热动力学膨胀和结构高度）是微形变反演的关键环节，同时也是最耗时的步骤。前者在影像集中筛选出高置信度的候选点并提供可靠的主辅影像差分相位信息，后者通过对候选点组网与网络空间搜索解算形变信息。然而，传统算法为单一时间跨度的批量处理，不具备时间维拓展性。即当新影像数据引入后，由于影像集组合策略和网络结构的改变，差分干涉预处理与网络模型解算只得重新处理，增加处理时间和计算量。另外，算法中模型参数解算的方法主要采用基于最大时序相关系数的解空间搜索法。该方法本质上为空间谱估计算法，即利用拟牛顿法的原理寻找目标函数的最优值（徐华平和王碧君，2014）；在缺乏先验知识的情况下，为了保证反演参数的精度，往往需要使用较大解空间进行全局遍历搜索，算法计算时间随解空间规模成倍数增加。尤其是在SAR大数据背景下，高重访率的卫星带来的海量观测数据，传统处理方法已然不能满足地表准实时监测的需求（ Li等，2018 Li等，2018 ）。为了在传统PSInSAR技术上拓宽时间维适用性，时序微变处理应向着序贯处理发展（ Ansari等，2017 Ansari等，2017 ）。本文提出了基于优化解空间的PSInSAR序贯微变反演方法，在保证未知参数反演精度基础上，能实现示范区高效形变反演。 2　研究方法及数据 1 2.1　研究区域 2 江苏省南京市是长江下游地区第二大城市，濒临长江，毗邻东海，盛行东南季风，气候常年温暖湿润，降水量充足，年平均降雨量达1090 mm，月平均气温从1月的2.7℃到7月的28.1℃。受益于怡人的气候条件和重要的地理位置，南京曾作为历代王朝都城长达500余年，是中国首批国家历史文化名城之一，也是中华文明的重要发祥地，积累了明孝陵和明城墙等弥足珍贵的文化遗产。 明城墙建于1366年—1393年，当时有宫城、皇城、京城和外郭城四重城墙；而现在通称的“南京城墙”一般特指其中的京城城墙（ 沈家祺和邵晓峰，2020 沈家祺和邵晓峰，2020 ）。该城墙地处南京腹地，东傍钟山，西纳石城，北控后湖，南拥外秦淮河，得山川之利，控江湖之势。一改旧时方、矩城墙旧制，形成了一个不规则的环绕面积约 <math id="M1"><mn mathvariant="normal">60</mn><mtext> </mtext><mi mathvariant="normal">k</mi><msup><mrow><mi mathvariant="normal">m</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703949&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703947&type= 10.07533360 3.04800010 的城池，当前城墙完整保存长度为25.1 km，是世界范围内最长、规模最大、保存原真性最好的古城墙，在中国的城市规划史上具有很高的研究价值（ 周成海和丁勇，2006 周成海和丁勇，2006 ）。本研究选取了南京明城墙月牙湖公园至神策门、狮子山至定淮门、西水关至东水关3段具有代表性的现存墙体本体及200 m缓冲区内墙体赋存环境为研究区（ 图1 图1 ）。综合考虑候选点构网密度、参数解算精度和形变累计等测度，扩大缓冲区范围，即选取了墙体2 km SAR影像感兴趣区进行数据处理，以提高PS点精度与算法鲁棒性（形变控制点位于中华门附近，如 图1 图1 所示）。 2.2　数据 2 本研究中使用的多源时空数据包括：（1）2015-01—2018-02的32景降轨Cosmo-SkyMed条带影像，影像中心入射角为34.0°，空间分辨率优于3 m。（2）美国地质调查局（USGS）的SRTM1高程数据，空间分辨率为1弧度秒（约30 m）。高程数据将用来模拟地形相位生成差分干涉图，并提供形变场地理编码经纬度转换，即将SAR影像坐标系转换为通用横轴墨卡托坐标系。（3）SAR影像成像日期温度及天气数据（ http：//lishi.tianqi.com/nanjing/index.html http://lishi.tianqi.com/nanjing/index.html ［2021-03-02］）；Cosmo-SkyMed卫星每次固定在当地时间18：02左右飞过研究区上空，获取温度数据可用于地物热幅值建模与参数解算。SAR数据各成像日期南京温度基本信息如 表1 表1 所示。 编号 成像日期 传感器 天气 成像温度/℃ 1 2015-01-17 CSK2 晴 3 2 2015-02-22 CSK4 晴转多云 8 3 2015-03-22 CSK2 晴转阴 15 4 2015-04-11 CSK4 多云转晴 14 5 2015-05-13 CSK4 多云 23 6 2015-06-14 CSK4 晴 29 7 2015-07-16 CSK4 多云 25 8 2015-08-17 CSK4 晴 23 9 2015-09-18 CSK4 多云转阵雨 21 10 2015-10-16 CSK2 晴 20 11 2015-11-17 CSK2 阴转小雨 14 12 2015-12-19 CSK2 阵雨 5 13 2016-01-20 CSK2 晴转多云 2 14 2016-03-12 CSK4 晴转多云 9 15 2016-04-13 CSK4 阵雨 13 16 2016-05-15 CSK4 阴转小雨 13 17 2016-06-16 CSK4 多云转阵雨 24 18 2016-07-14 CSK2 小雨转多云 26 19 2016-08-15 CSK2 晴 28 20 2016-09-20 CSK4 晴 18 21 2016-11-07 CSK4 阴转小雨 14 22 2016-12-09 CSK4 晴 8 23 2017-02-27 CSK4 晴 9 24 2017-04-16 CSK4 多云转阵雨 21 25 2017-05-18 CSK4 晴 22 26 2017-06-03 CSK4 多云转晴 25 27 2017-07-21 CSK4 晴 31 28 2017-09-19 CSK2 多云 23 29 2017-10-05 CSK2 阵雨转阴 15 30 2017-11-10 CSK4 晴转多云 15 31 2017-12-12 CSK4 晴 5 32 2018-02-14 CSK4 晴转多云 13 3　研究方法与数据处理 1 本研究拟采用传统PSInSAR和优化解空间PSInSAR序贯形变监测方法，结合多源时空数据进行示范地南京明城墙形变监测，并通过两者交叉验证，评估序贯方法的处理性能。 3.1　PSInSAR 2 PSInSAR方法面向覆盖相同区域的多景时序SAR影像，基于时空分布选取适宜的主影像，并以主影像为基准配准影像集合，进行差分干涉处理，提取在长时空基线下高相干、稳定PS点，并利用精细差分干涉相位模型，反演计算PS点高置信度的高程、形变等未知参数。 具体而言，对于 <math id="M2"><mi>N</mi><mo>+</mo><mn mathvariant="normal">1</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703969&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703968&type= 8.89000034 2.70933342 景原始单视复影像，经预处理后得到差分干涉相位集，基于振幅离散度与平均强度提取PS候选点，对于第 <math id="M3"><mi>i</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703974&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703971&type= 1.10066664 2.62466669 幅干涉图，各候选点的真实差分干涉相位模型为（ Chen等，2017 Chen等，2017 ） <math display="block" id="M4"><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mi>W</mi><mfenced separators="|"><mrow><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">e</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">A</mi><mi mathvariant="normal">P</mi><mi mathvariant="normal">S</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub><mo>+</mo><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">c</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703978&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703976&type= 64.85466003 7.95866632 <math display="block" id="M5"><mi>i</mi><mo>=</mo><mn mathvariant="normal">1,2</mn><mo>,</mo><mo>⋯</mo><mo>,</mo><mi>N</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703981&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703979&type= 18.45733452 3.21733332 式中， <math id="M6"><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703986&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703984&type= 2.96333337 4.40266657 为差分干涉相位， <math id="M7"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">e</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703990&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703988&type= 4.74133301 4.40266657 为高程相位， <math id="M8"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703994&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703992&type= 4.82600021 4.40266657 为形变相位， <math id="M9"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">a</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">m</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703998&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703996&type= 5.24933338 4.40266657 为大气延迟相位差， <math id="M10"><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">c</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704002&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23703999&type= 5.16466665 4.40266657 为去相干噪声相位， <math id="M11"><mi>W</mi><mo stretchy="false">(</mo><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704006&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704004&type= 6.09600019 3.47133350 表示相位缠绕。考虑到材料固有热力学膨胀，引入温度基线和热幅系数，高程相位和形变相位可表示为 <math display="block" id="M12"><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">e</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">e</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><msub><mrow><mi>ξ</mi></mrow><mrow><mi>i</mi></mrow></msub><mi>s</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704009&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704007&type= 21.67466736 6.26533365 <math display="block" id="M13"><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mfenced separators="|"><mrow><msub><mrow><mi>η</mi></mrow><mrow><mi>i</mi></mrow></msub><mi>υ</mi><mo>+</mo><msub><mrow><mi>ζ</mi></mrow><mrow><mi>i</mi></mrow></msub><mi>k</mi></mrow></mfenced><mo>+</mo><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi><mo>_</mo><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mo>_</mo><mi mathvariant="normal">n</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">e</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704012&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704011&type= 58.25066376 6.26533365 式中， <math id="M14"><msub><mrow><mi>η</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>/</mo><mi>λ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704016&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704014&type= 14.64733315 4.48733330 ， <math id="M15"><mtext> </mtext><msub><mrow><mi>t</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704020&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704018&type= 2.53999996 3.80999994 为时间基线，表示各影像与主影像之间的观测时间差，单位为年。 <math id="M16"><mi>υ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704025&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704022&type= 1.86266661 2.62466669 表示目标点沿斜距向的线性形变速率，单位为 <math id="M17"><mi mathvariant="normal">m</mi><mi mathvariant="normal">m</mi><mo>/</mo><mi mathvariant="normal">a</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704028&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704027&type= 8.04333401 2.70933342 。 <math id="M18"><msub><mrow><mi>ζ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><msub><mrow><mi>T</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>/</mo><mi>λ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704032&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704030&type= 15.15533447 4.40266657 ， <math id="M19"><msub><mrow><mi>T</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704037&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704035&type= 2.79399991 3.72533321 为温度基线，表示各影像与主影像之间气温差，单位为 <math id="M20"><mo>℃</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704041&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704038&type= 3.30200005 2.87866688 。 <math id="M21"><mi>k</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704044&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704043&type= 1.69333339 2.62466669 表示热幅系数，表示地物由温度变化导致的形变幅度，单位为 <math id="M22"><mi mathvariant="normal">m</mi><mi mathvariant="normal">m</mi><mo>/</mo><mo>℃</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704048&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704047&type= 9.65200043 2.87866688 ， <math id="M23"><msub><mrow><mfenced separators="|"><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">f</mi></mrow></msub></mrow></mfenced></mrow><mrow><mi>i</mi><mo>_</mo><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mo>_</mo><mi mathvariant="normal">n</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">e</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704052&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704050&type= 18.54199982 6.26533365 表示残余的非线性形变相位。 由于邻近PS点大气相位和残余的非线性形变相位在空间域是低频的，采用Delaunay三角网法布设候选点网络，并对网络中联通点进行二次差分，以此削弱大气相位与残余非线性形变相位的影响。二次差分相位表示为 <math display="block" id="M24"><mi mathvariant="normal">Δ</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mfenced separators="|"><mrow><msub><mrow><mi>ξ</mi></mrow><mrow><mi>i</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>s</mi><mo>+</mo><msub><mrow><mi>η</mi></mrow><mrow><mi>i</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>υ</mi><mo>+</mo><msub><mrow><mi>ζ</mi></mrow><mrow><mi>i</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>k</mi></mrow></mfenced><mo>+</mo><mi mathvariant="normal">Δ</mi><msub><mrow><mi>φ</mi></mrow><mrow><mi mathvariant="normal">i</mi><mo>_</mo><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704056&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704054&type= 61.12933731 5.16466665 式中， <math id="M25"><msub><mrow><mi>ξ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">2</mn><msub><mrow><mi>b</mi></mrow><mrow><mo>⊥</mo><mi>i</mi></mrow></msub><mo>/</mo><mi>λ</mi><mi>r</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mtext> </mtext><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704060&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704058&type= 25.39999962 4.40266657 ，表示高程向空间频率， <math id="M26"><mo>∆</mo><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi><mo>_</mo><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704064&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704062&type= 8.80533314 4.40266657 表示非线性形变相位、大气相位和噪声相位贡献和，且通常 <math id="M27"><mfenced open="|" close="|" separators="|"><mrow><mo>∆</mo><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi><mo>_</mo><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi></mrow></msub></mrow></mfenced><mo><</mo><mi mathvariant="normal">π</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704069&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704066&type= 18.79599953 5.24933338 。采用时序相关系数（NTM）作为模型匹配测度（ Ma和Lin，2016 Ma和Lin，2016 ）来解求参数和评价PS点质量： <math display="block" id="M28"><mi mathvariant="bold-italic">γ</mi><mfenced separators="|"><mrow><mi mathvariant="normal">Δ</mi><mi mathvariant="bold-italic">s</mi><mo>,</mo><mi mathvariant="normal">Δ</mi><mi mathvariant="bold-italic">υ</mi><mo>,</mo><mi mathvariant="normal">Δ</mi><mi mathvariant="bold-italic">k</mi></mrow></mfenced><mo>=</mo><munder><mrow><mi mathvariant="normal">M</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi></mrow><mrow><mi>γ</mi></mrow></munder><mfenced separators="|"><mrow><mfrac><mrow><mfenced open="|" close="|" separators="|"><mrow><mi>a</mi><msup><mrow><mfenced separators="|"><mrow><mi mathvariant="normal">Δ</mi><mi>s</mi><mo>,</mo><mi mathvariant="normal">Δ</mi><mi>υ</mi><mo>,</mo><mi mathvariant="normal">Δ</mi><mi>k</mi></mrow></mfenced></mrow><mrow><mi>H</mi></mrow></msup><mi>y</mi></mrow></mfenced></mrow><mrow><mi>N</mi></mrow></mfrac></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704072&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704070&type= 62.31467056 13.63133430 <math display="block" id="M29"><mi mathvariant="bold-italic">a</mi><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><mi mathvariant="normal">e</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">p</mi><mtext> </mtext><mfenced separators="|"><mrow><mi>j</mi><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mfenced separators="|"><mrow><msub><mrow><mi>ξ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi mathvariant="normal">Δ</mi><mi>s</mi><mo>+</mo><msub><mrow><mi>η</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi mathvariant="normal">Δ</mi><mi>υ</mi><mo>+</mo><msub><mrow><mi>ζ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mi mathvariant="normal">Δ</mi><mi>k</mi></mrow></mfenced></mrow></mfenced></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi mathvariant="normal">e</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">p</mi><mtext> </mtext><mfenced separators="|"><mrow><mi>j</mi><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mfenced separators="|"><mrow><msub><mrow><mi>ξ</mi></mrow><mrow><mi>N</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>s</mi><mo>+</mo><msub><mrow><mi>η</mi></mrow><mrow><mi>N</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>υ</mi><mo>+</mo><msub><mrow><mi>ζ</mi></mrow><mrow><mi>N</mi></mrow></msub><mi mathvariant="normal">Δ</mi><mi>k</mi></mrow></mfenced></mrow></mfenced></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704076&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704074&type= 59.60533142 20.23533440 式中， <math id="M30"><mi mathvariant="bold-italic">γ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704080&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704078&type= 2.20133328 3.38666677 表示时序相干系数， <math id="M31"><msup><mrow><mfenced separators="|"><mrow><mo>•</mo></mrow></mfenced></mrow><mrow><mi>H</mi></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704085&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704082&type= 5.84200001 4.40266657 表示共轭转置， <math id="M32"><mi mathvariant="bold-italic">y</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704089&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704087&type= 1.77800000 3.47133350 为二次差分复数序列。通过在离散空间中计算时序相关系数，并设定相干系数最小阈值 <math id="M33"><msub><mrow><mi>T</mi></mrow><mrow><mi>γ</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704093&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704091&type= 3.30200005 3.72533321 以淘汰低置信度的边，并优化网络结构。考虑解空间搜索的离散性和噪声点的存在，采用稳健的加权最小二乘法估计（M-estimator）获得连续域上的稳健参数增量（ Ma和Lin，2016 Ma和Lin，2016 ）；并进一步使用网络平差参量集成获取实验区域内PS点相对参考点的形变场和时序信息。 3.2　解空间优化序贯算法思路 2 算法立足差分干涉数据预处理与网络模型解算两个环节来提高计算效率，其主要处理步骤如下所示（ 图2 图2 ）。 输入： <math id="M34"><mi mathvariant="normal">新</mi><mi mathvariant="normal">增</mi><mi mathvariant="normal">影</mi><mi mathvariant="normal">像</mi><mi mathvariant="normal">序</mi><mi mathvariant="normal">列</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704103&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704101&type= 22.26733208 3.55599999 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704103&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704103&type=middle 120.01499939 141.98599243 以及前序PSInSAR估算参数 <math id="M35"><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704123&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704121&type= 3.55599999 3.72533321 <math display="block" id="M36"><mtext> </mtext><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>=</mo><msub><mrow><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mover accent="true"><mi>s</mi><mo>̂</mo></mover></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mrow><mover accent="true"><mi>s</mi><mo>̂</mo></mover></mrow><mrow><mi>Q</mi></mrow></msub></mtd></mtr></mtable><mtext>      </mtext><mtable><mtr><mtd><msub><mrow><mover accent="true"><mi>υ</mi><mo>̂</mo></mover></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mrow><mover accent="true"><mi>υ</mi><mo>̂</mo></mover></mrow><mrow><mi>Q</mi></mrow></msub></mtd></mtr></mtable><mtext>      </mtext><mtable><mtr><mtd><msub><mrow><mover accent="true"><mi>k</mi><mo>̂</mo></mover></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mrow><mover accent="true"><mi>k</mi><mo>̂</mo></mover></mrow><mrow><mi>Q</mi></mrow></msub></mtd></mtr></mtable></mrow></mfenced></mrow><mrow><mi>Q</mi><mo>×</mo><mn mathvariant="normal">3</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704111&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704109&type= 37.59199905 16.67933464 Q （PS点数）维向量 <math id="M37"><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704115&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704114&type= 3.55599999 3.72533321 存储的前序微变数值。 输出： <math id="M38"><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704120&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704118&type= 3.55599999 3.72533321 流程： （1）分析新获取影像时空基线分布，与主影像配准、干涉与差分干涉。 （2）网络重组与外拓：重新提取PS候选点并构建全局Delaunay三角网。检索PS候选点与前序PS点的重叠度，确定内部边与外延边的网络组成：其中内部边仅含前序PS点，外延边为前序PS点和新提取PS候选点的组合。然后，对所有PS点依据网络拓扑进行相位二次差分。 （3）参数继承：外延边中所有未知候选点的初始参数参考最邻近PS点，内部边点初始参数等于 <math id="M39"><msub><mrow><mi mathvariant="bold-italic">X</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704123&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704121&type= 3.55599999 3.72533321 。然后进一步将点参数依照网络拓扑关系转换为边上的高程修正增量 <math id="M40"><mo>∆</mo><mi>s</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704126&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704124&type= 3.38666677 2.62466669 、线性形变速率增量 <math id="M41"><mo>∆</mo><mi>ν</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704131&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704128&type= 3.80999994 2.62466669 和热幅增量 <math id="M42"><mo>∆</mo><mi>k</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704135&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704133&type= 3.72533321 2.62466669 。 （4）解空间采样与划分：将所得网络参数纳入模型匹配测度计算初始时序相干系数，设置相应阈值将网络边分为高相关边（例如，时间相干系数>0.7）、欠相关边（时间相干系数位于［0.4，0.7］）和失相关边（时间相干系数<0.4）；基于相干系数空间谱分布，对各相关边设置相应解空间。 （5）优化解空间搜索：1）可行域搜索：对解空间进行采样分析，探测最优解分布子区间；2）子空间最优解搜索：基于上述子空间，进行精细遍历搜索，获取全局最优解。 （6）M-estimator与岭回归：将离散解空间网络转化为连续域内稳健解，然后基于预设的参照点，通过带权岭回归分析方法对网络最优参数进行求解，获取各节点相对参考点的最小二乘最优解。 3.3　数据处理 2 本研究首先以28景影像（2015-01—2017-09）作为初始数据，基于卫星时空分布与地面气候状况选取2016-08-15期影像为参考主影像，以尽量减少时空去相干和大气相位系统误差。使用基于幅度的互相关方法对SAR数据集进行了联合配准，然后进行差分雷达干涉与时序微变反演。由于常规基于幅度离差等方法在提取PS候选点过程中未能综合考虑线性遗产的几何特征及其SAR成像姿态参数，导致局部线性遗产沿线PS候选点分布过于稀疏（例如墙体走向与卫星飞行方向垂直）；因此在案例中将墙体上所有散射点合并至候选PS点集，以实现对线性墙体微变的高覆盖提取。 为测试优化解空间序贯方法在明城墙文化遗产景观微变监测的能力及潜在瓶颈，模拟新增4景影像（影像时空基线，如 图3 图3 所示），分别使用传统PSInSAR与PSInSAR序贯算法进行形变估算。 对于传统PsInSAR算法，本文不再赘述，将着眼于解空间优化序贯算法的具体实施方案：新增4景影像数据后（ 图3 图3 ），影像集时空基线分布未发生实质性改变。因此无需重复处理前序影像序列，只需对新增影像进行配准、裁剪和差分干涉等预处理。 事实上，预处理流程中影像集配准耗时最长，当初始影像序列数量级远大于新增影像集时，主影像在序贯前、后两个时间跨度下可表征为连贯性。新影像集的无缝嵌入可有效提高算法计算效率。 网络模型解算阶段，序贯算法通过未知参数搜索机制和空间结构优化来提高参数解求效率。在初始目标点参数已知条件下，新影像引入时，首先采取与传统PSInSAR同样的策略选取PS候选点，并基于候选点进行Delaunay全局构网。接下来，保留三角网络中所有包含前序PS点的网络边，依据前序已知解参数及最邻近插值法求取外延边所有候选点参数，并基于拓扑关系转化为网络边上的高程增量 <math id="M43"><mtext> </mtext><mo>∆</mo><mi>s</mi><mtext> </mtext></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704142&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704141&type= 4.91066647 2.79399991 、形变增量（ <math id="M44"><mo>∆</mo><mi>v</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704144&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704143&type= 3.55599999 2.62466669 、 <math id="M45"><mo>∆</mo><mi>k</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704148&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704146&type= 3.72533321 2.62466669 ），进而确定由内部边和外延边共同组成Delaunay网络未知参数增量（高程、形变和热幅系数）先验值。考虑到：（1）形变（线性形变、材料热膨胀）是持续变化的微小量，即无明显地质运动（地震、滑坡）的情况下，其在时空维度可表征为高相关；（2）因城市区域地物高程存在突变，导致插值点高程值不准确；加之多散射体信号混叠，PS候选点高程初始值可与实际数值发生偏差。鉴于此，本研究首先计算初始网络边时序相干系数，并设置0.7与0.4经验双阈值，将网络边划分为高相关边，低相关边与失相关边，如 图4 图4 所示；在本研究中，有47%的网络边仍保持高相关，且仅有少部分（21%）网络边失相关，三者数量之比为47∶32∶21。若采用时序批处理方式对整体网络进行解算将不可避免的增加无效计算量。然后基于时序相干系数确认网络边空间搜索机制，采用策略及对应原理如下：时序相干系数搜索空间谱本质上是多个正弦波的相互叠加，即在复数空间内表征为距原点不同长度的离散点分布。在单散射体主导的空间谱内，相干系数可存在多个波峰（但仅存单个主峰），因此参数搜索空间可与主峰波谱宽度对应（ 图5 图5 ）。对于高相关边，其网络边参数位于空间谱波峰邻近区域，只需在初始参数上设置上下微小浮动区间（一般不超过原有全局解空间的20%）；对于低相关边，其可能位于主峰或旁瓣区域，需相应扩大浮动区间（经验值为原有全局解空间的50%），其中 <math id="M46"><mo>∆</mo><mi>s</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704151&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704149&type= 3.38666677 2.62466669 应上调增长幅度，以减少高程突变的影响；对于失相关边，高程将采用原始解空间，浮动空间进一步增大（经验值为原有全局解空间的70%）。对于不包含PS点的网络边，综合考虑其对整体网络结构的影响以及其小数量，则可忽略其网络边集合。 1,2;3,4;5; 搜索机制确立后，研究采用全局—局部由粗到精分支限界搜索策略，以避免整个解空间的枚举：首先以较大搜索步长（例如，未知参数测量精度的10倍）对解空间进行采样分析，确定最优解潜在分布区间（注：搜索采样步长不宜过大，以避免落入局部最优区间）；同时根据区间内最大时间相干系数确定网络边的取舍，即当 <math id="M47"><msubsup><mrow><mi>T</mi></mrow><mrow><mi>γ</mi></mrow><mrow><mi>'</mi></mrow></msubsup><mo><</mo><msub><mrow><mi>T</mi></mrow><mrow><mi>γ</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704201&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704199&type= 11.09133339 4.06400013 时，直接舍弃该网络边。其次，基于潜在分布区间，使用精细搜索步长（预定未知参数测量精度），对最优解空间进行遍历搜素，实现未知参数最优解的获取。综上所述，算法权衡了大范围全局搜索和小范围局部寻优能力，提高了最优解估算收敛速度。 4　结果与分析 1 传统PSInSAR和PSInSAR序贯方法在目标区域内分别获取到37219，37096个PS点，研究区形变速率结果如 图6 图6 — 图8 图8 所示（GF-2为背景影像），形变速率为负表示目标沿远离雷达的方向运动；反之，形变速率为正表示目标沿朝向雷达的方向运动。其中西水关至东水关段研究区由于墙体走势与卫星飞行轨道近似垂直，在雷达图像上表现为弱散射体，大部分PS点分布在城墙城门与建筑物区域，墙体本体沿线PS点覆盖较为稀疏。从结果中可以看出，大部分PS整体平均沉降速率位于-2—+2 mm/a，揭示在成像时段内，墙体及其赋存环境不存在显著形变。 1,2;3,4;5,6; 1,2;3,4; 1,2;3,4; GPS或精确水准测量是交叉验证InSAR测量精度的科学方法（ Chen等，2019 Chen等，2019 ）。然而，由于缺乏大地测量数据，无法利用地面测量资料对形变探测结果进行定量化验证。考虑到传统PSInSAR在城市地表沉降估算的精度和可靠性，本研究采用两种算法的交叉比对来间接测评PSInSAR序贯算法的性能。研究将两种算法获取的线性形变速率组成散点图（ 图9 图9 ），可以发现两种算法的形变速率结算结果具有高度一致性，相关系数为0.9966，置信度大于99.58%。另外，选取狮子山至定淮门段3个典型点（仪风门、挹江门、华严岗门）的时序形变信息组成折线图（ 图10 图10 ），可知两种算法的时序形变解算结果保持一致。 1,2;3;4; 研究进一步随机选用两条剖线（ 图11 图11 — 图12 图12 ），展示了两种算法形变反演在选定的水平和纵向剖面中亦获得了一致的位移测量。 1,2;3; 1,2;3; 为了评估序贯方法的形变测量精度，以传统PSInSAR形变量为基准，定义偏差 <math id="M48"><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">f</mi><mi mathvariant="normal">f</mi></mrow></msub><mtext> </mtext></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704372&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704370&type= 6.01133299 3.80999994 ，其计算公式如下： <math display="block" id="M49"><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">d</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">f</mi><mi mathvariant="normal">f</mi></mrow></msub><mo>=</mo><mfenced open="|" close="|" separators="|"><mrow><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">P</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">A</mi><mi mathvariant="normal">R</mi></mrow></msub><mo>-</mo><msub><mrow><mi>σ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">q</mi><mi mathvariant="normal">u</mi><mi mathvariant="normal">e</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">t</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">l</mi><mtext> </mtext><mi mathvariant="normal">P</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">n</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">A</mi><mi mathvariant="normal">R</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704378&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=23704376&type= 45.63533020 5.24933338 图13 图13 为两种方法所有PS点的偏差分量值，横轴坐标表示PS点的序号。图中最大偏差分量值为3 mm/a，大部分PS点偏差小于0.5 mm/a，总体误差在0—1 mm/a。 如前文所述，PSInSAR序贯算法旨在数据差分干涉处理和时序形变反演两大环节计算效率的提升。为评价PSInSAR序贯方法的计算性能，同时排除计算机硬件水平对方法效率的影响，新增影像数据后，开展两者方法平均耗时测评（ 图14 图14 ）；结果表明序贯方法能总体实现约10倍效率的提升。 综上所述，本研究提出的序贯方法具有如下优点：（1）相对传统时序方法，序贯方法考虑了影像集时空基线分布，可实现新影像无缝嵌入，避免了数据集整体流程重复处理。（2）序贯方法充分利用了前序PS点参数，并对整体网络进行筛选；专注于已有PS点的形变变化，舍弃了不包含PS点的网络边，进而简化了整体网络结构；同时在网络参数求解中，序贯算法对保留网络采用不同的搜索机制，可减少网络边的冗余搜索。（3）优化解空间搜索兼顾全局搜索与局部寻优能力，指导解空间向可行域搜索，有效提高了搜索空间计算效率。 5　讨论 1 如 图6 图6 — 图8 图8 ，研究区内形变量呈高斯正态分布，区域整体沉降趋于稳定；其中最大沉降区位于集庆门南部内环西线段，沉降速率值为-6 mm/a。此外，区域内零星分布有少量形变异常值，局部墙体需要进行重点关注。为探讨沉降驱动机制，本研究聚焦于明城墙集庆门与东水关城门形变靶区，将目标线性沉降速率与材料热动力学运动异质性（热膨胀系数）分离（ 图15 图15 和 图18 图18 ），旨在揭示地面变形和结构不稳定性两者对沉降的贡献。 1,2; 1,2; 1,2; 明城墙集庆门靶区的线性形变和热膨胀参数较为显著，且呈现一定正相关（ 图15 图15 ）。考虑到该区域城墙无植被覆盖，植物根系生物破坏可忽略。城市市政路网修建改造是触发形变的直接驱动力：即临时性土方开挖与回填导致城墙局部松动。特别是2013年内环西线的修建（集庆门隧道深基坑开挖）引发城墙下方土层应力状态改变，实地观察到裂缝与水泥填补脱落现象（ 图16 图16 ）。 另外，线性形变与热膨胀数值空间局部正相关也可源自PSInSAR算法：（1）模型解算算法本身的局限性，当线性形变与热膨胀形变分量两者在时空维度存在耦合时，最大时序相干系数法可能无法分离两者贡献。（2）解空间为有限三维离散型分布，存在1/2搜索步长系统误差；经网络传递后，该误差可累积放大。（3）PSInSAR算法基于单永久散射体点SPS（Single Permanent Scatter），没有考虑叠掩双永久散射体点DPS（Double Permanent Scatters）现象（ 图17 图17 ），在DPS复合出现时，可导致微变参数的错误解算。 图18 图18 为明城墙东水关城门靶区（线性、热膨胀）沉降分布图，由于该段墙体内部多采用黄土夹碎石夯实（ 杨新华，2006 杨新华，2006 ），自然雨水及长期风蚀作用导致墙体明显退化。此外，植物根系的季节性生长对岩石层产生不同程度的挤压破坏（ 图19 图19 ）。 1,2; 另一方面，材料热胀冷缩在空间上可表现为一定差异性。即由于城墙底部采用条石砌筑，上部（顶部雉堞和女墙）采用城砖砌筑；这种高度向材料和结构的差异，在温度应力作用下可发生墙体不规则臌胀与变形。 6　结论 1 遥感大数据的到来，对文化遗产地形变近实时监测及其时空演化模式诠析提出更高要求。针对传统PSInSAR 技术在形变时序更新面临的技术瓶颈，本文提出了一种基于优化解空间的PSInSAR序贯处理算法。算法立足于差分干涉数据预处理与网络模型解算两个环节，充分借鉴前序形变反演先验知识，通过目标点参数的有效继承、网络重组与外拓、解空间采样与划分、全局/局部搜索并行等途径改善了解空间收敛过程；同时应用M-estimator与岭回归提高了反演结果的稳健性。 研究以南京明城墙文化遗产为研究对象，使用PSInSAR和序贯PSInSAR算法开展形变监测和性能比对。交叉验证结果证实了序贯处理方法的有效性，即该方法在确保形变估算精度的基础上，有效提高了参数解算效率（例如，本试验提升效率约10倍）；揭示了序贯PSInSAR方法在未来InSAR大数据近实时处理和形变时序更新的应用潜力。 最后，研究结合南京明城墙案例，对墙体异常形变进行了驱动力分析。研究表明，周期性材料热胀冷缩、趋势性线性形变及其他贡献因子的分离，有助形变模型参数精确反演，以及参数错误解算成因的探究；从而为PSInSAR建筑遗产形变监测和算法模型优化指明了方向。 参考文献（References） Ansari H , Zan F D and Bamler R . 2017 . Sequential estimator: toward efficient InSAR time series analysis . IEEE Transactions on Geoscience and Remote Sensing , 55 ( 10 ): 5637 - 5652 [ DOI: 10.1109/TGRS.2017.2711037 http://dx.doi.org/10.1109/TGRS.2017.2711037 ] Sequential estimator: toward efficient InSAR time series analysis Chen F L . 2015 . Principle, application and prospects of synthetic aperture radar (SAR) remote sensing in archaeology . Remote Sensing Technology and Application , 30 ( 5 ): 835 - 841 Principle, application and prospects of synthetic aperture radar (SAR) remote sensing in archaeology 陈富龙 . 2015 . 雷达遥感考古机理、应用与展望 . 遥感技术与应用 , 30 ( 5 ): 835- 841 [ DOI: 10.11873/j.issn.1004-0323.2015.5.0835 http://dx.doi.org/10.11873/j.issn.1004-0323.2015.5.0835 ] 雷达遥感考古机理、应用与展望 Chen F L , Wu Y H , Zhang Y M , Parcharidis I , Ma P F , Xiao R Y , Xu J , Zhou W , Tang P P and Foumelis M . 2017 . Surface motion and structural instability monitoring of Ming Dynasty City walls by two-step Tomo-PSInSAR approach in Nanjing City, China . Remote Sensing , 9 ( 4 ): 371 [ DOI: 10.3390/rs9040371 http://dx.doi.org/10.3390/rs9040371 ] Surface motion and structural instability monitoring of Ming Dynasty City walls by two-step Tomo-PSInSAR approach in Nanjing City, China Chen F L , Zhou W , Chen C F and Ma P F . 2019 . Extended D-TomoSAR displacement monitoring for Nanjing (China) city built structure using high-resolution TerraSAR/TanDEM-X and Cosmo SkyMed SAR data . Remote Sensing , 11 ( 22 ): 2623 [ DOI: 10.3390/rs11222623 http://dx.doi.org/10.3390/rs11222623 ] Extended D-TomoSAR displacement monitoring for Nanjing (China) city built structure using high-resolution TerraSAR/TanDEM-X and Cosmo SkyMed SAR data Chen F L , Zhou W , Xu H , Parcharidis I , Lin H and Fang C Y . 2020 . Space technology facilitates the preventive monitoring and preservation of the Great Wall of the Ming Dynasty: a comparative study of the Qingtongxia and Zhangjiakou sections in China . IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing , 13 : 5719 - 5729 [ DOI: 10.1109/JSTARS.2020.3023297 http://dx.doi.org/10.1109/JSTARS.2020.3023297 ] Space technology facilitates the preventive monitoring and preservation of the Great Wall of the Ming Dynasty: a comparative study of the Qingtongxia and Zhangjiakou sections in China Deng B , Guo H D , Wang C L and Nie Y P . 2010 . Applications of remote sensing technique in archaeology: a review . Journal of Remote Sensing , 14 ( 1 ): 187 - 206 Applications of remote sensing technique in archaeology: a review 邓飚 , 郭华东 , 王长林 , 聂跃平 . 2010 . 遥感技术在考古中的应用综述 . 遥感学报 , 14 ( 1 ): 187- 206 [ DOI: 10.11834/jrs.20100115 http://dx.doi.org/10.11834/jrs.20100115 ] 遥感技术在考古中的应用综述 Ferretti A , Prati C and Rocca F . 2001 . Permanent scatterers in SAR interferometry . IEEE Transactions on Geoscience and Remote Sensing , 39 ( 1 ): 8 - 20 [ DOI: 10.1109/36.898661 http://dx.doi.org/10.1109/36.898661 ] Permanent scatterers in SAR interferometry Gabriel A K , Goldstein R M and Zebker H A . 1989 . Mapping small elevation changes over large areas: differential radar interferometry . Journal of Geophysical Research : Solid Earth , 94 ( B7 ): 9183 - 9191 [ DOI: 10.1029/JB094iB07p09183 http://dx.doi.org/10.1029/JB094iB07p09183 ] Mapping small elevation changes over large areas: differential radar interferometry Li S , Du S Y , Yu A X , Dong Z , Zhang X L and Zhu X X . 2018 . A PS- InSAR sequential processing method based on network reconstruction // 2018 International Conference on Microwave and Millimeter Wave Technology (ICMMT) . Chengdu, China : IEEE : 1 - 3 [ DOI: 10.1109/ICMMT.2018.8563850 http://dx.doi.org/10.1109/ICMMT.2018.8563850 ] A PS- InSAR sequential processing method based on network reconstruction Ma P F and Lin H . 2016 . Robust detection of single and double persistent scatterers in urban built environments . IEEE Transactions on Geoscience and Remote Sensing , 54 ( 4 ): 2124 - 2139 [ DOI: 10.1109/TGRS.2015.2496193 http://dx.doi.org/10.1109/TGRS.2015.2496193 ] Robust detection of single and double persistent scatterers in urban built environments Shen J Q and Shao X F . 2020 . Research on the design of guide system under the cultural background——Taking the outer wall of Ming Dynasty in Nanjing as an example . Art and Literature for the Masses , ( 2 ): 141 - 142 Research on the design of guide system under the cultural background——Taking the outer wall of Ming Dynasty in Nanjing as an example 沈家祺 , 邵晓峰 . 2020 . 文化底蕴下的导视系统设计研究——以南京明城墙外郭为例. 大众文艺, (2): 141 - 142 [DOI: 10.3969/j.issn.1007-5828.2020.02.088] 文化底蕴下的导视系统设计研究——以南京明城墙外郭为例. 大众文艺, (2): Xu H P and Wang B J . 2014 . Parameter estimation method of PS-DInSAR surface deformation measurement based on optimized solution space search method . China, CN 104091064 Parameter estimation method of PS-DInSAR surface deformation measurement based on optimized solution space search method A (徐华平 , 王碧君 . 2014 . 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法. 中国, CN104091064A) 基于优化解空间搜索法的PS-DInSAR地表形变测量参数估计方法. 中国, CN104091064A) Yang X H . 2006 . City wall of the Ming Dynasty in Nanjing . Nanjing : Nanjing University Press : 55 - 594 City wall of the Ming Dynasty in Nanjing 杨新华 . 2006 . 南京明城墙 . 南京 : 南京大学出版社 : 55 - 594 南京明城墙 Zhou C H and Ding Y . 2006 . Thoughts on the history and existing condition of Ming dynasty’s rampart in Nanjing . Shanxi Architecture , 32 ( 24 ): 43 - 44 Thoughts on the history and existing condition of Ming dynasty’s rampart in Nanjing 周成海 , 丁勇 . 2006 . 南京明城墙的历史与现状思考 . 山西建筑 , 32 ( 24 ): 43- 44 [ DOI: 10.3969/j.issn.1009-6825.2006.24.027 http://dx.doi.org/10.3969/j.issn.1009-6825.2006.24.027 ] 南京明城墙的历史与现状思考 10.11834/jrs.20211098 1007-4619（2021）12-2381-15 王成，陈富龙，周伟，余华芬，吴迪.2021.南京明城墙序贯PSInSAR形变估计.遥感学报，25（12）： 2381-2395 王成，陈富龙，周伟，余华芬，吴迪.2021.南京明城墙序贯PSInSAR形变估计.遥感学报，25（12）： 2381-2395 Wang C，Chen F L，Zhou W，Yu H F and Wu D. 2021. Sequential PSInSAR approach for the deformation monitoring of the Nanjing Ming Dynasty City Wall. National Remote Sensing Bulletin， 25（12）：2381-2395 Wang C，Chen F L，Zhou W，Yu H F and Wu D. 2021. Sequential PSInSAR approach for the deformation monitoring of the Nanjing Ming Dynasty City Wall. National Remote Sensing Bulletin， 25（12）：2381-2395 国家自然科学基金（编号：41771489）；国家重点研发计划国际合作重点专项（编号：2017YFE0134400） Supported by National Natural Science Foundation of China (No. 41771489); National Key Research and Development Program of China (No.2017YFE0134400) 2021-03-02 2021-12-07 2022-01-11 2021-09-20 遥感学报 12 25