復(fù)雜多變地物散射特征區(qū)域InSAR融合方法研究

夏錦 王崢輝

摘 要：【目的】PS-InSAR和SBAS-InSAR在應(yīng)對復(fù)雜多變地物散射特征沉降區(qū)域時均存在解算局限性，有必要通過對比兩者技術(shù)差異和適用性，開展融合方法研究?！痉椒ā扛鶕?jù)不同地物散射特征區(qū)域PS點(diǎn)分布的相干系數(shù)統(tǒng)計特性，尋找與強(qiáng)地物散射特征區(qū)域強(qiáng)關(guān)聯(lián)的PS點(diǎn)數(shù)據(jù)簇，再根據(jù)不同散射特征區(qū)域PS點(diǎn)分布的密度差異，采用空間聚類算法識別覆蓋城鎮(zhèn)用地的數(shù)據(jù)簇，并采用三角網(wǎng)格法對照全國土地利用現(xiàn)狀數(shù)據(jù)確定監(jiān)測數(shù)據(jù)融合的邊界，剔除邊界外低密度、低質(zhì)量的PS點(diǎn)數(shù)據(jù)簇，在邊界外采用SBAS-InSAR解算出SDFP點(diǎn)數(shù)據(jù)，得到最終結(jié)果。【結(jié)果】數(shù)據(jù)融合后，弱地物散射區(qū)域高質(zhì)量監(jiān)測點(diǎn)數(shù)據(jù)大幅增加，強(qiáng)地物散射區(qū)域內(nèi)高質(zhì)量、高密度PS點(diǎn)被保留，該部分PS點(diǎn)解算位置精度更高，可在建成區(qū)實(shí)現(xiàn)小尺度精細(xì)化監(jiān)測，最終實(shí)現(xiàn)根據(jù)區(qū)域內(nèi)的可變散射特征自動選擇匹配的干涉測量數(shù)據(jù)的結(jié)果。【結(jié)論】InSAR融合方法在應(yīng)對復(fù)雜多變地物散射特征區(qū)域時，可兼顧監(jiān)測點(diǎn)的數(shù)量和質(zhì)量。

關(guān)鍵詞：地面沉降；SBAS-InSAR；PS-InSAR；地物散射；融合方法

中圖分類號：TU196.2? ? 文獻(xiàn)標(biāo)志碼：A? ? ?文章編號：1003-5168（2024）03-0052-08

DOI：10.19968/j.cnki.hnkj.1003-5168.2024.03.011

Research on InSAR Fusion Method for Scattering Characteristic

Regions of Complex and Variable Objects

XIA Jin1? WANG Zhenghui2

（1.Nanjing Metro Construction Co.， Ltd.， Nanjing 210017，China; 2.Jiangsu Hengtong Engineering

Investigation and Testing Co.， Ltd.， Nanjing 210012，China）

Abstract： [Purposes] Both PS InSAR and SBAS InSAR have computational limitations when dealing with complex and variable terrain scattering characteristics in subsidence areas， so it is necessary to carry out research on fusion methods by comparing the technical differences and applicability of the two. [Methods] According to the statistical characteristics of the coherence coefficient of the distribution of PS points in different scattering feature areas， find the PS point data cluster that is Strongly correlated material to the strong scattering feature area. Then， according to the density difference of the distribution of PS points in different scattering feature areas， this paper uses the spatial clustering algorithm to identify the data cluster covering urban land. This paper uses the Triangle mesh method to determine the boundary of monitoring data fusion against the national land use status data， and eliminates the low density low quality PS point data cluster， and use SBAS InSAR to calculate SDFP point data outside the boundary and obtain the final result.[Findings] After data fusion， there is a significant increase in high-quality monitoring point data in weak terrain scattering areas， while high-quality and high-density PS points in strong terrain scattering areas are retained. This part of the PS points has higher accuracy in calculating their positions， and small-scale fine monitoring can be achieved in built-up areas. Ultimately， matching interferometric measurement data can be automatically selected based on the variable scattering features in the area. [Conclusions] This InSAR fusion method can balance the quantity and quality of monitoring points when dealing with complex and variable terrain scattering feature areas.

Keywords： ground subsidence; SBAS-InSAR; PS-InSAR; ground scattering; fusion method

0 引言

合成孔徑雷達(dá)差分干涉測量（Differential Interferometric Synthetic Aperture Radar， D-InSAR）技術(shù)可以監(jiān)測大范圍、全天候、亞厘米級精度的地面沉降，但易受時空失相干與大氣效應(yīng)的影響。為實(shí)現(xiàn)長時序高精度地面監(jiān)測，時序InSAR技術(shù)被提出。時序InSAR可分為基于高相干點(diǎn)目標(biāo)的單主影像時序InSAR和基于高相干面目標(biāo)的多主影像時序InSAR。前者的代表是Ferretti等［1］于2000年左右提出的PS-InSAR（Permanent Scatter InSAR）；后者的代表則是Berardino等［2］提出的SBAS-InSAR（Small Baseline Subset InSAR）。近年來，兩種方法體系不斷發(fā)展完善，為地面沉降監(jiān)測提供了技術(shù)基礎(chǔ)。

PS-InSAR可以從建構(gòu)筑物、橋梁、各種角散射器效應(yīng)所構(gòu)成的高相干目標(biāo)中提取大量永久散射體，廣泛應(yīng)用于城市區(qū)域地面沉降監(jiān)測。但是，失相干會使PS點(diǎn)密度降低，無法形成緊密間隔的參考網(wǎng)格，影響大氣干擾的估算和去除，測量精度難以保障［3-4］。對此，已有大量學(xué)者投入到提高PS-InSAR中PS點(diǎn)密度的研究中。針對判別PS點(diǎn)相干性的優(yōu)化算法，Sadeghi等［5］將PS-InSAR與單基線偏振相干性優(yōu)化相結(jié)合來監(jiān)測地面沉降，通過提高PS點(diǎn)的相干性增加PS點(diǎn)數(shù)量。針對識別PS點(diǎn)的優(yōu)化算法，黃長軍等［6］提出了一種基于經(jīng)驗(yàn)?zāi)B(tài)分析的PS點(diǎn)檢測方法，在PS點(diǎn)檢測過程中降低了誤判和遺漏的可能性。Perissin等［7］通過對半相干點(diǎn)進(jìn)行探測來提高PS點(diǎn)密度，但需要更多的計算時間。在地物參照點(diǎn)輔助識別PS點(diǎn)方面，Zhu等［8］通過在低相干區(qū)域布設(shè)角反射器，提高了PS點(diǎn)的空間密度。上述研究在一定程度上克服了失相干使PS點(diǎn)密度降低的問題，在低相干區(qū)域識別出更多PS點(diǎn)，但是在復(fù)雜多變的地物散射特征沉降區(qū)域，PS點(diǎn)的分布仍然存在較大密度差異，對監(jiān)測結(jié)果的局部精度存在影響［9-10］。

SBAS-InSAR通過剔除基線較大的干涉對，形成多個小基線集合，提高了時空采樣率，降低了時空失相干的概率［11-12］。SBAS-InSAR反演過程對SAR影像數(shù)據(jù)采用自適應(yīng)多視處理，從而提高影像的信噪比，降低了去相干噪聲對形變估算的影響，提高了干涉質(zhì)量，但會降低空間分辨率，影響監(jiān)測精度［13］。目前，許多學(xué)者圍繞提高SBAS-InSAR監(jiān)測結(jié)果空間分辨率進(jìn)行了研究。針對多元數(shù)據(jù)融合，F(xiàn)rancesca等［13］在處理交通網(wǎng)絡(luò)和采礦點(diǎn)沿線局部變形時，通過采用高分辨率SUAV技術(shù)，提高了監(jiān)測結(jié)果的局部分辨率。針對干涉目標(biāo)的提取算法，Mora等［14］提出了一種基于相干系數(shù)的相干目標(biāo)選擇方法，可以在時序影像數(shù)量較少的情況下提取干涉目標(biāo)。在數(shù)據(jù)視數(shù)優(yōu)化方面，Lanari等［15］提出了一種基于單視處理的算法來探測干涉目標(biāo)，同樣可以起到提高空間分辨率的效果。SBAS-InSAR通常需要通過最小成本流方法或分支切割方法進(jìn)行相位展開，市區(qū)干涉圖的干涉條紋會導(dǎo)致相位展開困難，以及提取散射體的相位信息誤差增大，在城市中的應(yīng)用還存在局限性［16］。

目前，兩種時序技術(shù)在應(yīng)對復(fù)雜多變地物散射特征沉降區(qū)域時均存在局限性，關(guān)于時序InSAR適用范圍局限性的研究還處于初始階段。近年來，許多學(xué)者提出了融合不同時序InSAR的思路，旨在增加相干點(diǎn)目標(biāo)的空間采樣率、提高分布密度，獲得更可靠的監(jiān)測結(jié)果。姜德才等［17］提出了融合兩種方法散射體的選取算法，通過構(gòu)建Delaunay三角網(wǎng)來估算形變速率，提高了目標(biāo)提取的正確率和最終結(jié)果的精度。聶運(yùn)菊等［18］通過篩選與提取分布于研究區(qū)內(nèi)的穩(wěn)定PS點(diǎn)作為SBAS-InSAR處理過程中的地面控制點(diǎn)，進(jìn)行軌道精煉與形變反演。上述研究圍繞全局目標(biāo)點(diǎn)選取的數(shù)量和質(zhì)量進(jìn)行討論，未考慮兩種方法在應(yīng)對沉降區(qū)域地物散射特征復(fù)雜多變時的解算局限性。因此，本文從兩種方法應(yīng)對沉降區(qū)域地物散射特征復(fù)雜多變的適用性出發(fā)，開展融合方法的研究。

1 InSAR技術(shù)

1.1 PS-InSAR技術(shù)

從覆蓋相同區(qū)域不同時間段的SAR影像數(shù)據(jù)中選取一幅作為公共主影像，其余影像作為從影像配準(zhǔn)到主影像上形成影像對，使所有影像對之間的空間和時間去相關(guān)效應(yīng)最小，然后去除平地相位并獲取差分干涉相位圖。其中，[φi]表示第i幅差分干涉圖的相位，具體見式（1）。

[φi=φdef _i+φtopo_i+φatm_i+φnoise_i] （1）

式中：[φdef _i]為LOS方向的形變相位；[φtopo_i]為殘余地形相位；[φatm_i]為大氣延遲相位；[φnoise_i]為其他噪聲相位。

以PS點(diǎn)為節(jié)點(diǎn)，建立不規(guī)則三角網(wǎng)，采用差分方法削弱空間距離相關(guān)的誤差（如大氣延遲），估算PS網(wǎng)絡(luò)中每條基線的參數(shù)，以及PS點(diǎn)的LOS向線性形變和高程誤差。采用振幅離差閾指數(shù)（ADI）篩選穩(wěn)定的候選點(diǎn)PSC（Permanent Scatterer Candidates），建立相鄰PS點(diǎn)相位差分模型并求解模型參數(shù)。

第[i]幅差分干涉圖的相鄰PS點(diǎn)相位差表示見式（2）。

[φidiff=4πλTiΔν+4πB⊥iΔελRsinθ+Δωi] （2）

式中：[Δν]為相鄰PS點(diǎn)之間的形變速率增量；[Δε]為高程誤差增量；[Δωi]為殘余相位誤差；[B⊥i]、[Ti]分別為干涉對[i]的時空基線；[λ]為波長；[R]為天線相位中心到目標(biāo)的距離；[θ]為雷達(dá)入射角。同時假設(shè)[Δωi<π]。

在得到各個PS點(diǎn)對間的相對形變速率和高程誤差之后，通過積分獲得各個PS點(diǎn)的形變速率和高程誤差。

1.2 SBAS-InSAR技術(shù)

獲取覆蓋研究區(qū)的多幅SAR影像，確定時空基線閾值自由組合成干涉對。假設(shè)覆蓋研究區(qū)的SAR影像數(shù)量為N，獲取時間依次是[t0]，[t1]，…，[tN?1]。確定時空基線閾值，形成M幅差分干涉圖，M滿足式（3）。

[M2≤N≤M（M?1）2] （3）

假設(shè)第[i]幅差分干涉圖的相位[φi]，見式（4）。

[φi=φdef _i+φtopo_i+φatm_i+φnoise_i] （4）

式中：[φdef _i]表示LOS方向的形變相位；[φtopo_i]為殘余地形相位；[φatm_i]為大氣延遲相位；[φnoise_i]為其他噪聲相位。

SBAS-InSAR通過逐個計算像元來獲取差分干涉圖中各個像元在時間序列上的形變。為消除原始長時序形變中的大氣延遲和相位噪聲誤差，對原始形變進(jìn)行線性回歸，反演平均形變速率，模型見式（5）。

[A?=δ?] （5）

式中：[?]為N+1個成像時刻的雷達(dá)影像上未知形變相位點(diǎn)組成的矩陣；[δ?]為M幅差分干涉圖中觀測相位值矩陣；A為一個M×N的系數(shù)矩陣，形變相位通過采用最小范數(shù)意義的最小二乘法求解。

考慮外部DEM數(shù)據(jù)導(dǎo)致的精度誤差，相位中存在殘余地形相位[ε]，可建立方程見式（6）。

[Bν+Cε=δ?] （6）

式中：[B]為M×N的矩陣，[C]為與空間基線相關(guān)M×1的系數(shù)矩陣。

利用矩陣奇異值分解方法可估算線性形變速率和高程誤差。從先前計算的原始長時序形變中減去反向線性形變，剩下非線性地面形變、大氣延遲和隨機(jī)相位噪聲。非線性形變可以通過低通濾波器分離。最后，線性形變與非線性形變之和即為長時序形變信息。

1.3 技術(shù)差異

PS-InSAR的位移監(jiān)測目標(biāo)是穩(wěn)定的點(diǎn)狀散射體，即尺寸較小、具有較好幾何特征且散射特性穩(wěn)定的物體。在對SAR影像進(jìn)行差分干涉處理時，不考慮臨界基線的限制，不進(jìn)行多普勒帶寬濾波和光譜位移濾波。PS點(diǎn)因自身算法沒有光譜位移、基線去相關(guān)或多普勒帶寬的去相關(guān)，不需要距離向和方位向的濾波。PS-InSAR對每個像素的時間序列分別進(jìn)行計算，不執(zhí)行任何相位解纏。這兩個特征使獲得的監(jiān)測結(jié)果空間分辨率最大，即相鄰像素的監(jiān)測值相對獨(dú)立。

SBAS-InSAR的位移監(jiān)測（差分相位）目標(biāo)為幾何特征不明顯的物體。反演過程通過采用自適應(yīng)多視處理，將處理后的相位作為其優(yōu)化相位，提高強(qiáng)度影像的信噪比，但會降低分辨率［12］。為減少大量失相干現(xiàn)象對監(jiān)測目標(biāo)的影響，算法流程中包含了多普勒帶寬濾波和光譜位移濾波處理。雖然在選擇PS點(diǎn)和SDFP像素點(diǎn)時運(yùn)用了相同算法，但是因?yàn)镻S點(diǎn)是對單幅主影像且沒經(jīng)過光譜位移濾波，SDFP像素點(diǎn)是對多主影像且經(jīng)過了光譜位移濾波，所以選擇出的像素對是不同的。干涉圖配準(zhǔn)后會進(jìn)行相位解纏，解纏精度與干涉相位的空間連續(xù)性有關(guān)，并且也是誤差傳播的潛在來源。

2 數(shù)據(jù)處理

2.1 研究區(qū)概況

宿州西水源地位于宿州城區(qū)所在的埇橋區(qū)西南方向，占地面積135 km2。區(qū)域地質(zhì)為新生界松散覆蓋層，覆蓋層主要由粉砂、細(xì)砂、黏土質(zhì)砂、黏土、砂質(zhì)黏土、鈣質(zhì)黏土等組成，自上而下可分成4個含水層（組）和3個隔水層（組）。以合徐高速為界，水源地東側(cè)建筑密度較大，城鎮(zhèn)化水平較高，存在大量建構(gòu)筑物和硬化場地；西側(cè)多為河渠、水田、林地和旱地，植被覆蓋表面受自然因素影響變化較大，且建構(gòu)筑物較少。

宿州西水源地合徐高速公路以西區(qū)域設(shè)有地下水開采區(qū)，供水層位主要是第二含水層及第三含水層上部地下水，有地下水開采井共110口，水井密度高達(dá)6～8眼/km2，其中97%為中深層地下水，其余為淺層地下水。由于井位布局不合理，地下水超采嚴(yán)重，地下水補(bǔ)給有限，已形成較大的地下水水位降落漏斗，造成區(qū)域性地面沉降，跨越城市建成區(qū)和非建成區(qū)，若沉降量持續(xù)增加，可能進(jìn)一步導(dǎo)致地面高程損失，引起市政、道路設(shè)施破壞，建筑物結(jié)構(gòu)破壞，抗震能力降低等問題。

鑒于宿州西水源地內(nèi)既存在建構(gòu)筑物密度較高的建成區(qū)，又存在以河渠、水田、林地和旱地為主的非建成區(qū)，且由于地下水降落漏斗擴(kuò)大導(dǎo)致的地面沉降跨越建成區(qū)和非建成區(qū)，選取宿州西水源地包括群井在內(nèi)的20 km×20 km矩形區(qū)域作為研究區(qū)，如圖1所示。

2.2 研究數(shù)據(jù)概況

Sentinel-1A衛(wèi)星是歐洲航天局2014年發(fā)射的地球觀測衛(wèi)星，該衛(wèi)星載有C波段合成孔徑雷達(dá)，采用漸進(jìn)掃描成像模式生成全球影像，重訪觀測周期為12 d。衛(wèi)星影像采用IW成像模式，250 km幅寬和VV極化方式。本文選取了覆蓋宿州西水源地研究區(qū)的38幅升軌數(shù)據(jù)。

因?yàn)楦缮鎴D中的殘差條紋包括有誤差的軌道信息，會帶來系統(tǒng)誤差，需要在InSAR數(shù)據(jù)預(yù)處理流程中利用衛(wèi)星軌道數(shù)據(jù)來校正。該試驗(yàn)采用的數(shù)據(jù)時間跨度較大，對監(jiān)測結(jié)果的精度要求較高，因此應(yīng)使用精密定軌星歷數(shù)據(jù)。精密定軌星歷數(shù)據(jù)會每天生成一個文件，該數(shù)據(jù)發(fā)布時間會延遲21 d，定位精度優(yōu)于5 cm。

2.3 數(shù)據(jù)處理結(jié)果

在38幅SAR影像中，PS-InSAR處理流程選取2019-01-22期影像作為主影像，將主影像與所有從影像依次配準(zhǔn)形成37個干涉對，然后經(jīng)過差分干涉處理、PSC的識別和提取，最終反演出形變。SBAS-InSAR處理流程將所有SLC影像配準(zhǔn)重采樣后通過設(shè)置時空基線閾值，生成272個干涉對，然后經(jīng)過差分干涉處理、軌道精煉與重去平，最終反演出形變。PS-InSAR和SBAS-InSAR的處理結(jié)果如圖2和圖3所示。由圖可知，研究區(qū)內(nèi)PS-InSAR結(jié)果顯示年平均形變速率介于-35.4～10.4 mm·a-1之間；SBAS-InSAR結(jié)果顯示年平均形變速率介于-43.2～4.4 mm·a-1之間。

在城市建成區(qū)，由于人工建筑等強(qiáng)散射體分布較密集，PS點(diǎn)分布密度高，可形成緊密間隔的參考網(wǎng)格，監(jiān)測結(jié)果可達(dá)到理想狀態(tài)下的精度。SBAS-InSAR通常需要通過最小成本流方法進(jìn)行相位展開，城市建成區(qū)干涉圖的復(fù)雜干涉條紋會導(dǎo)致相位展開困難，提取散射體的相位信息誤差增大。在城市非建成區(qū)，失相干使PS點(diǎn)密度降低，無法形成緊密間隔的參考網(wǎng)格，影響大氣干擾的估算和去除，測量精度難以保障。相比之下，SBAS-InSAR能較好地克服時空失相干，在該類型區(qū)域獲得了更多監(jiān)測數(shù)據(jù)，降低了零星監(jiān)測點(diǎn)的偶然誤差。

在應(yīng)對復(fù)雜多變地物散射特征沉降區(qū)域時，PS-InSAR在弱地物散射區(qū)域的監(jiān)測精度易衰減，SBAS-InSAR在弱地物散射區(qū)域克服失相干能力更強(qiáng)，數(shù)據(jù)密度更高，但在強(qiáng)地物散射區(qū)域的監(jiān)測精度不及理想狀態(tài)下PS-InSAR的監(jiān)測精度。單一時序InSAR難以兼顧監(jiān)測點(diǎn)的數(shù)量和質(zhì)量。對此，本文對融合兩種時序InSAR技術(shù)的方法進(jìn)行研究，旨在根據(jù)區(qū)域內(nèi)的可變散射特征自動選擇匹配的干涉測量數(shù)據(jù)。

3 融合方法研究

3.1 相干系數(shù)統(tǒng)計特性研究

由PS-InSAR監(jiān)測結(jié)果可知，PS-InSAR可在人工建筑表面解算出更多的PS點(diǎn)，且相干性高于河渠、水田、林地等區(qū)域解算出的PS點(diǎn)。城市建成區(qū)以人工建筑為主，城市非建成區(qū)以植被為主，在建成區(qū)和非建成區(qū)高相干干涉目標(biāo)的空間分布密度和相干性存在差異。基于此特征確定數(shù)據(jù)融合的邊界，采用基于密度的聚類算法識別建成區(qū)范圍內(nèi)的高相干PS點(diǎn)，根據(jù)建成區(qū)和非建成區(qū)內(nèi)PS點(diǎn)相干性強(qiáng)弱的差異，提高篩選建成區(qū)內(nèi)PS點(diǎn)的準(zhǔn)確性。

以中國科學(xué)院2020年完成的GlobeLand30 V2020 全國土地利用現(xiàn)狀數(shù)據(jù)為參照，對PS點(diǎn)的相干系數(shù)統(tǒng)計特性進(jìn)行研究，在ArcGIS中載入PS點(diǎn)數(shù)據(jù)，得到PS點(diǎn)在不同類型用地上的密度特征，如圖4所示。研究區(qū)主要有城鎮(zhèn)用地、農(nóng)田、農(nóng)村居民點(diǎn)和旱地四種用地類型。PS點(diǎn)主要分布于城鎮(zhèn)用地和農(nóng)村居民點(diǎn)，水田和旱地區(qū)域分布較少。將研究區(qū)劃分為城鎮(zhèn)用地和非城鎮(zhèn)用地兩部分，對所有PS點(diǎn)按照相干系數(shù)值進(jìn)行統(tǒng)計分析。相干系數(shù)以0.01為間隔，計算各區(qū)間內(nèi)非城鎮(zhèn)用地和城鎮(zhèn)用地區(qū)域PS點(diǎn)數(shù)量的比值，統(tǒng)計結(jié)果如圖5所示。

結(jié)果顯示，隨著相干系數(shù)增加，非城鎮(zhèn)用地和城鎮(zhèn)用地區(qū)域內(nèi)PS點(diǎn)數(shù)量比例呈下降趨勢，即非城鎮(zhèn)用地內(nèi)PS點(diǎn)數(shù)量的下降速率大于城鎮(zhèn)用地。當(dāng)相干系數(shù)大于0.9時，非城鎮(zhèn)用地和城鎮(zhèn)用地區(qū)域內(nèi)PS點(diǎn)的數(shù)量比例驟降，由0.31下降至0.21，此時非城鎮(zhèn)用地和城鎮(zhèn)用地區(qū)域內(nèi)PS點(diǎn)的密度差異發(fā)生明顯變化，原因在于河渠、水田、林地和旱地等非建成區(qū)內(nèi)相干系數(shù)大于0.9的永久散射體數(shù)量較少，故將篩選閾值設(shè)為0.9。以相干系數(shù)0.9為閾值對所有PS點(diǎn)進(jìn)行篩選，篩選結(jié)果顯示，高相干PS點(diǎn)的分布滿足以下密度特征：在城市建成區(qū)集中且點(diǎn)距小，在城市非建成區(qū)較分散且點(diǎn)距大。

3.2 基于PS點(diǎn)密度的聚類分析

空間聚類是一種研究空間內(nèi)對象之間距離相似性的聚類方法?？臻g對象屬性特征相似或空間位置關(guān)聯(lián)度較高的數(shù)據(jù)集被定義為“簇”，不同簇之間空間對象屬性特征或空間位置關(guān)聯(lián)度差異較大。DBSCAN （Density-Based Spatial Clustering of Applications with Noise）是一種基于密度的空間聚類算法，該算法考慮空間內(nèi)數(shù)據(jù)對象間遠(yuǎn)近及疏密關(guān)系的結(jié)構(gòu)特征，通過全局參數(shù)[Eps]和[MinPts]來定義密度標(biāo)準(zhǔn)查找聚類子集，根據(jù)空間對象的密度尋找被低密度區(qū)域分離的高密度區(qū)域。傳統(tǒng)DBSCAN算法需要人為確定聚類參數(shù)，聚類精度取決于參數(shù)選擇的是否合理。本研究從數(shù)據(jù)集本身的統(tǒng)計特性出發(fā)來確定參數(shù)，采用基于非參數(shù)核密度估計理論的Kernel-DBSCAN算法對PS點(diǎn)數(shù)據(jù)集進(jìn)行聚類分析。該算法利用核函數(shù)描述數(shù)據(jù)對象的分布特征，在識別數(shù)據(jù)對象中的噪聲點(diǎn)時準(zhǔn)確性較高。該算法利用核函數(shù)的實(shí)施過程如下。

假設(shè)數(shù)據(jù)集包含[x1]，[ x2]，…，[xi]個對象，概率密度函數(shù)為[f]，公式見式（7）。

[f?x=1n∑ni=1K?（x?xi）=1n?∑ni=1Kx?xi?] （7）

式中：h為帶寬；[K（x）]為核函數(shù)，同時[K（x）]滿足以下條件，見式（8）、式（9）。

[Kx≥0，Kxdx=1] （8）

[xKxdx=0，x2Kxdx>0] （9）

其中，帶寬會影響擬合結(jié)果的精度。一般地，在聚類分析中選擇均平方積分誤差函數(shù)[MISE?]來優(yōu)化帶寬，公式見式（10）。

[MISE?=Efx?fx2dx] （10）

在弱假設(shè)下，則有式（11）至式（14）。

[MISE?=AMISE?+o（1n?+?4）] （11）

[AMISE?=RKn?+14m2K2?4R（f''）] （12）

[RK=Kx2dx>0] （13）

[m2K=x2Kxdx>0] （14）

最小化[MISE?]等價于最小化[AMISE?]，求偏導(dǎo)并令導(dǎo)數(shù)等于0，則有式（15）、式（16）。

[??? AMISE?=?RKn?2+m2K2?3Rf''=0] （15）

[?AMISE=RK15m2K25R（f''）15n15] （16）

根據(jù)[Eps]的預(yù)估值，采用數(shù)學(xué)期望的方法，根據(jù)距離矩陣[DISTn×n]，在給定數(shù)據(jù)集中計算[MinPts]的值，公式見式（17）。

[MinPts=1n∑ni=1Pi] （17）

式中：[Pi]為對象[i]的[Eps]鄰域內(nèi)包含的樣本個數(shù)。

采用Kernel-DBSCAN算法對PS點(diǎn)數(shù)據(jù)集進(jìn)行聚類分析，以期獲得覆蓋城鎮(zhèn)用地的數(shù)據(jù)簇。試驗(yàn)首先計算所有PS點(diǎn)之間的空間距離數(shù)據(jù)集，再分別計算[Eps]和[MinPts]，此過程基于MATLAB編程語言實(shí)現(xiàn)。根據(jù)PS點(diǎn)間的空間距離數(shù)據(jù)集繪制其核密度估計曲線，結(jié)果如圖6所示。由圖6可知，橫軸表示PS點(diǎn)間的距離，縱軸表示介于該距離附近PS點(diǎn)的密度程度，經(jīng)過了無量綱標(biāo)準(zhǔn)化處理，利用核密度估計曲線刻畫PS點(diǎn)間不同距離的密度特征。

分析圖6可知，核密度曲線只出現(xiàn)了一次峰值，峰值點(diǎn)對應(yīng)的距離為4 286.9 m，這表示點(diǎn)與點(diǎn)之間的距離在4 286.9 m左右的數(shù)量最多，數(shù)據(jù)集存在唯一的中心簇。由圖2可知，研究區(qū)內(nèi)城鎮(zhèn)化發(fā)展集中于東北象限，與聚類結(jié)果一致，因此可認(rèn)為當(dāng)聚類半徑達(dá)到4 286.9 m時效果最佳。密度曲線的峰值既可以表征單個簇在距離上的分布特征，也可以表征簇與簇之間在距離上的分布特征。當(dāng)存在多城鎮(zhèn)化中心時，核密度曲線將相應(yīng)出現(xiàn)多個峰值，第一個峰值表示單個簇空間內(nèi)對象之間的緊密程度，第二個峰值表示不同簇空間內(nèi)對象之間的緊密程度。此時，選擇第一個峰值對應(yīng)的距離值作為[Eps]效果最佳。根據(jù)[Eps]的預(yù)估值，采用數(shù)學(xué)期望的方法，根據(jù)距離矩陣計算[MinPts]的值。

最終，Kernel-DBSCAN算法選取[Eps]=4 286.9、[MinPts]=1 983作為輸入?yún)?shù)，對PS點(diǎn)數(shù)據(jù)集進(jìn)行聚類分析，結(jié)果如圖7所示。聚類結(jié)果顯示，高相干PS點(diǎn)數(shù)據(jù)集根據(jù)密度特征被劃分為2類簇，城鎮(zhèn)用地內(nèi)存在2類簇，城鎮(zhèn)用地邊界基本無噪聲。

3.3 數(shù)據(jù)融合與結(jié)果分析

以GlobeLand30 V2020全國土地利用現(xiàn)狀數(shù)據(jù)為參照，Kernel-DBSCAN算法能夠識別覆蓋城鎮(zhèn)用地的數(shù)據(jù)簇，即城鎮(zhèn)用地數(shù)據(jù)簇，僅將其分裂為兩部分，有效合并同一類型的簇。在MATLAB中對城鎮(zhèn)用地數(shù)據(jù)簇的邊界進(jìn)行提取。首先，將所有PS點(diǎn)作為節(jié)點(diǎn)創(chuàng)建三角網(wǎng)格，然后對所有三角網(wǎng)格進(jìn)行以下計算，如果三角網(wǎng)格的邊不存在公共邊則將其定義為“邊緣”，遍歷所有三角網(wǎng)格得到最終的“邊緣”集合，將其作為城鎮(zhèn)用地數(shù)據(jù)簇的初始邊界，結(jié)果如圖8所示。

對三角網(wǎng)格進(jìn)行以下計算，假設(shè)三角網(wǎng)格邊界上的最大邊長為D，如果三角形位于邊界上的邊大于D，則刪除該三角形。遍歷三角網(wǎng)格的邊界，如果三角形位于邊界上的邊在當(dāng)前迭代中小于最大邊長，則停止遍歷。對三角網(wǎng)格邊界的最大邊長進(jìn)

行計算，依次剔除最大邊并計算邊界修正后對應(yīng)的城鎮(zhèn)用地數(shù)據(jù)簇覆蓋面積。計算城鎮(zhèn)用地數(shù)據(jù)簇區(qū)域與土地利用現(xiàn)狀數(shù)據(jù)中城鎮(zhèn)用地區(qū)域的重疊率[ω]，見式（18），將其作為衡量邊界修正后的精度指標(biāo)。

[ω=S1∩S2S1∪S2] （18）

式中：S1為土地利用現(xiàn)狀數(shù)據(jù)城鎮(zhèn)用地區(qū)域面積；S2為城鎮(zhèn)用地數(shù)據(jù)簇覆蓋面積；[ω]為重疊率。

邊界修正精度評價見表1。試驗(yàn)結(jié)果表明如果不經(jīng)過邊界修正，城鎮(zhèn)用地數(shù)據(jù)簇覆蓋面積與土地利用現(xiàn)狀數(shù)據(jù)城鎮(zhèn)用地區(qū)域面積的重疊率為83.31%。當(dāng)剔除最大邊長4 604.77 m時，重疊率達(dá)到最大值84.62%，此時可認(rèn)為邊界精度達(dá)到最優(yōu)。將此時城鎮(zhèn)用地數(shù)據(jù)簇覆蓋區(qū)域作為城市建成區(qū)，其他區(qū)域作為非建成區(qū)，修正后的邊界作為數(shù)據(jù)融合的邊界。

在ArcGIS中利用融合邊界對所有PS點(diǎn)進(jìn)行篩選，剔除非建成區(qū)低質(zhì)量、低密度的PS點(diǎn)，并采用SDFP像素點(diǎn)數(shù)據(jù)，最終得到融合方法的監(jiān)測結(jié)果，如圖9所示。對研究區(qū)內(nèi)的PS點(diǎn)和SDFP像素點(diǎn)數(shù)量進(jìn)行統(tǒng)計。在建成區(qū)一共保留了110 308個PS點(diǎn)，PS點(diǎn)密度由原來的375個/ km2提高至858個/ km2，提高了128.8%；在非建成區(qū)等弱散射特征區(qū)域探測到的PS點(diǎn)總量僅為47 054個，密度為348個/ km2，SBAS-InSAR監(jiān)測結(jié)果顯示該區(qū)域內(nèi)存在170 469個SDFP像素點(diǎn)，密度為1 261個/ km2，采用SDFP像素點(diǎn)后，監(jiān)測點(diǎn)密度提高262.3%，在監(jiān)測公路、橋梁附近的區(qū)域性地面沉降時監(jiān)測點(diǎn)更為連續(xù)。PS-InSAR監(jiān)測結(jié)果在宿永路和宿渦路之間的區(qū)域存在因?yàn)闀r空失相干導(dǎo)致數(shù)據(jù)缺失的現(xiàn)象，SBAS-InSAR在該區(qū)域監(jiān)測到了沉降速率介于30.00～43.17 mm·a-1的區(qū)域。

綜上所述，數(shù)據(jù)融合后，弱地物散射區(qū)域低質(zhì)量、低密度的PS點(diǎn)被剔除。采用SDFP像素點(diǎn)大幅提高了監(jiān)測點(diǎn)密度，在強(qiáng)地物散射區(qū)域的高質(zhì)量、高密度PS點(diǎn)被保留，該部分PS點(diǎn)解算位置精度更高，可以在建成區(qū)實(shí)現(xiàn)小尺度精細(xì)化監(jiān)測。

4 結(jié)論

本研究提出了一種在復(fù)雜地物散射特征區(qū)域內(nèi)融合PS-InSAR和SBAS-InSAR干涉測量數(shù)據(jù)的地面沉降監(jiān)測新方法。在PS-InSAR解算中獲取PS點(diǎn)，篩選出高相干值的PS點(diǎn)，采用DBSCAN算法對PS點(diǎn)坐標(biāo)分布聚類分析。PS點(diǎn)分布密集區(qū)對應(yīng)建成區(qū)，PS點(diǎn)稀疏區(qū)對應(yīng)非建成區(qū)，兩者邊界線參照全國土地利用現(xiàn)狀數(shù)據(jù)劃定并優(yōu)化，最終確定在建成區(qū)采用PS-InSAR監(jiān)測數(shù)據(jù)，在非建成區(qū)采用SBAS-InSAR監(jiān)測數(shù)據(jù)。圍繞宿州西水源地融合方法的應(yīng)用情況，得出以下結(jié)論。

①根據(jù)建成區(qū)和非建成區(qū)內(nèi)PS點(diǎn)相干性強(qiáng)弱的差異來提高篩選建成區(qū)內(nèi)PS點(diǎn)的準(zhǔn)確性。弱地物散射區(qū)域內(nèi)相干系數(shù)大于0.9的永久散射體數(shù)量較少，若相干系數(shù)大于0.9，PS點(diǎn)出現(xiàn)在人工建筑上的概率大幅提升，可作為高相干PS點(diǎn)的篩選閾值。

②K-DBSCAN算法可以高效識別城鎮(zhèn)用地內(nèi)的數(shù)據(jù)簇和噪聲，利用三角網(wǎng)格法對照全國土地利用現(xiàn)狀數(shù)據(jù)可進(jìn)一步提高融合邊界的精度。由于宿州西水源地研究區(qū)具有唯一的城市化中心，因此對于發(fā)生地面沉降的多城市化中心研究區(qū)K-DBSCAN算法的聚類效果有待進(jìn)一步研究。

③融合方法將弱地物散射區(qū)域低質(zhì)量、低密度的PS點(diǎn)被剔除，采用SDFP像素點(diǎn)大幅提高了監(jiān)測點(diǎn)密度，在強(qiáng)地物散射區(qū)域的高質(zhì)量、高密度PS點(diǎn)被保留，該部分PS點(diǎn)解算位置精度更高，可以在建成區(qū)實(shí)現(xiàn)小尺度精細(xì)化監(jiān)測，兼顧了復(fù)雜多變地物散射特征區(qū)域監(jiān)測點(diǎn)的數(shù)量和質(zhì)量。

參考文獻(xiàn)：

［1］FERRETTI? A，PRATI C， ROCCA F. Permanent scatterers in SAR interferometry［J］. IEEE Transactions on Geoscience and Remote Sensing， 2001， 39（1）： 8-20.

［2］BERARDINO P，F(xiàn)ORNARO G， LANARI R， et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms［J］. IEEE Transactions On Geoscience and Remote Sensing， 2002，40（11）：2375-2383.

［3］GAMA F F， PARADELLA W R， MURA J C， et al. Advanced DINSAR analysis on dam stability monitoring： a case study in the Germano mining complex （Mariana， Brazil） with SBAS and PSI techniques［J］. Remote Sensing Applications， 2019，16（1）：100267.

［4］YU M， HANSSEN R F. Deformation parameter estimation in low coherence areas using a multisatellite InSAR approach［J］. IEEE Transactions on Geoence and Remote Sensing， 2013，53（8）：249-252.

［5］SADEGHI Z， ZOEJ M， MULLER J P. Combination of persistent scatterer interferometry and single-baseline polarimetric coherence optimisation to estimate deformation rates with application to tehran basin［J］. Pfg–Journal of Photogrammetry， Remote Sensing and Geoinformation Science， 2017，85（5）：327-340.

［6］黃長軍， 胡紀(jì)元， 楊亞夫. 一種基于經(jīng)驗(yàn)?zāi)B(tài)分解的永久散射體探測方法［J］. 光學(xué)學(xué)報， 2019，39（5）：371-381.

［7］PERISSIN D ， WANG T. Repeat-pass SAR interferometry with partially coherent targets［J］. IEEE Transactions on Geoscience and Remote Sensing， 2012，50（1）：271-280.

［8］ZHU W，ZHANG Q，DING X L， et al. Landslide monitoring by combining of CR-InSAR and GPS techniques［J］. Advances in Space Research， 2014，53（3）：430-439.

［9］QU F， ZHANG Q， ZHAO C，et al.Simultaneous estimation of building height and ground deformation over Xi'an City， China using multi-temporal InSAR method［J］.IEEE， 2016.

［10］姚佳明， 姚鑫， 吳作啟， 等. 基于InSAR三維分解技術(shù)的貴州省貞豐某煤礦地下采空區(qū)反演［J］. 工程地質(zhì)學(xué)報， 2020，28（4）：852-866.

［11］OGUSHI F， MATSUOKA M， DEFILIPPI M，et al. Improvement of persistent scatterer interferometry to detect large non-linear displacements with the 2π ambiguity by a non-parametric approach［J］.Remote Sensing （Basel， Switzerland）， 2019，11（21）：2467.

［12］高二濤， 范冬林， 付波霖， 等. 基于PS-InSAR和SBAS技術(shù)監(jiān)測南京市地面沉降［J］. 大地測量與地球動力學(xué)， 2019，39（2）：158-163.

［13］FRANCESCA C，VANESSA B，ALEXANDER D， et al. Mapping ground instability in areas of geotechnical infrastructure using satellite InSAR and small UAV surveying： a case study in Northern Ireland［J］. Geosciences， 2017，7（3）：1-25.

［14］MORA O， MALLORQUI J J， BROQUETAS A. Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images［J］. IEEE Transactions On Geoscience & Remote Sensing， 2003，41（10）：2243-2253.

［15］LANARI R， MORA O， MANUNTA M， et al. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms［J］. IEEE Transactions on Geoscience & Remote Sensing， 2004，42（7）：1377-1386.

［16］林琿， 馬培峰， 王偉璽. 監(jiān)測城市基礎(chǔ)設(shè)施健康的星載MT-InSAR方法介紹［J］. 測繪學(xué)報， 2017，46（10）：1421-1433.

［17］姜德才， 張繼賢， 張永紅， 等. 百年煤城地表沉降融合PS/SBAS InSAR監(jiān)測：以徐州市為例［J］. 測繪通報， 2017（1）：58-64.

［18］聶運(yùn)菊， 熊佳誠， 程朋根， 等. 結(jié)合PS特征點(diǎn)的SBAS地表形變監(jiān)測［J］. 測繪通報， 2020（4）：91-95.