孫凱，孟國杰，洪順英，董彥芳，李成龍

(1.中國地震局地震預(yù)測(cè)研究所 地震預(yù)測(cè)重點(diǎn)實(shí)驗(yàn)室，北京 100036；2.中國地震局地質(zhì)研究所 地震動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100029)

0 引言

近20多年來，合成孔徑雷達(dá)差分干涉技術(shù)(differential interferometric synthetic aperture radar， DInSAR)在測(cè)量地震形變場(chǎng)、監(jiān)測(cè)山體滑坡、城市地表沉降以及火山運(yùn)動(dòng)等領(lǐng)域發(fā)揮出重要的作用，其中應(yīng)用最有效且最廣泛的是獲取地震形變場(chǎng)[1-4]。利用DInSAR技術(shù)能夠輕松獲得受災(zāi)區(qū)域大面積的、連續(xù)空間覆蓋的高精度形變數(shù)據(jù)，有效彌補(bǔ)了常規(guī)全球定位系統(tǒng)(global positioning system， GPS)觀測(cè)中人工布設(shè)臺(tái)站數(shù)目不足、易受災(zāi)害損壞以及信息采集困難等缺陷，為研究斷層活動(dòng)性質(zhì)、地殼形變特征提供了定量化的基礎(chǔ)數(shù)據(jù)[5-7]。形變場(chǎng)數(shù)據(jù)集為反演模型提供穩(wěn)定的數(shù)據(jù)源。為了提高反演計(jì)算效率和減少信息冗余，InSAR形變數(shù)據(jù)需要進(jìn)行降采樣處理，并且在盡可能提高數(shù)據(jù)壓縮率的同時(shí)，確保反演數(shù)據(jù)源中形變信息的完整性和可靠性。

不同于常規(guī)的數(shù)據(jù)壓縮技術(shù)，形變數(shù)據(jù)降采樣就是通過最少的數(shù)據(jù)量同時(shí)體現(xiàn)原形變場(chǎng)整體區(qū)域狀況和局部重要細(xì)節(jié)特征。對(duì)于較為平滑的地表形變圖，Pritchard等[8]提出一種規(guī)則采樣法，通過計(jì)算鄰近像素值的均值，近似代表局部像元包含的形變特征，該方法計(jì)算簡(jiǎn)便，但必須在影像細(xì)節(jié)特征和數(shù)據(jù)壓縮程度二者間進(jìn)行取舍，不適用于形變特征差異大的震區(qū)形變場(chǎng)。Simons等[9]采用影像像元分解思路，有效提高數(shù)據(jù)壓縮率的同時(shí)，保留了形變場(chǎng)的細(xì)節(jié)特征，但該方法并未考慮采樣點(diǎn)數(shù)據(jù)質(zhì)量。張靜等[10]提出一種基于InSAR數(shù)據(jù)空間相關(guān)性建立協(xié)方差函數(shù)來設(shè)定四叉樹分解參數(shù)的方法，通過減少相位突變提高數(shù)據(jù)連續(xù)性，消除數(shù)據(jù)中的部分高頻噪聲。

現(xiàn)有降采樣方法缺少依據(jù)InSAR干涉特性評(píng)價(jià)形變信息的可靠性。由于雷達(dá)側(cè)視成像幾何特性和InSAR數(shù)學(xué)模型的限制，目標(biāo)區(qū)域地表反射特性隨時(shí)間變化造成的時(shí)間失相干和空間基線導(dǎo)致的空間失相干都會(huì)影響InSAR干涉質(zhì)量，干涉質(zhì)量決定了干涉相位信息的可接受程度，反映了最終形變信息的可靠性[11]。因此，相干性可作為描述影像相位噪聲水平和評(píng)價(jià)干涉圖質(zhì)量的指標(biāo)。基于此，本文提出一種顧及InSAR相干系數(shù)的概率采樣算法，估算自適應(yīng)濾波后的干涉像對(duì)相干性系數(shù)作為四叉樹分割約束參數(shù)，對(duì)待分割窗口內(nèi)相干性系數(shù)進(jìn)行自適應(yīng)核運(yùn)算得到相干表征量，并根據(jù)歸一化后的形變信息和相干性信息，合理分配形變梯度和相干表征量之間的權(quán)重，計(jì)算接受概率，獲取窗口分割指令。同時(shí)，依據(jù)數(shù)字高程模型和原始形變圖無效值(not a number，NaN)生成掩膜文件，去除采樣結(jié)果中的虛假信息，整體提高了反演數(shù)據(jù)源的質(zhì)量。本文通過獲取中國臺(tái)灣地區(qū)美濃Mw6.7級(jí)地震的形變場(chǎng)來驗(yàn)證該算法的可行性，并對(duì)采樣點(diǎn)數(shù)據(jù)質(zhì)量進(jìn)行定性與定量分析。

1 顧及相干系數(shù)的概率采樣算法

本節(jié)在回顧經(jīng)典四叉樹分割算法的基礎(chǔ)上，引入相干性系數(shù)作為額外約束參數(shù)，并結(jié)合當(dāng)前待分割窗口的像元分布狀況建立自適應(yīng)概率采樣模型。

1.1 四叉樹影像分割

InSAR形變場(chǎng)降采樣的目的是以最少的采樣點(diǎn)達(dá)到原始影像的整體最佳逼近。對(duì)于規(guī)則采樣等傳統(tǒng)降采樣方法而言，較少的采樣數(shù)目和具體的細(xì)節(jié)特征不可兼顧，因此，必須合理化參考影像內(nèi)容特征進(jìn)行自適應(yīng)采樣。作為基于遙感圖像內(nèi)容均勻性檢測(cè)的壓縮算法，自適應(yīng)四叉樹影像分割成為形變場(chǎng)降采樣的主流手段。

自適應(yīng)算法主要思想是根據(jù)給定數(shù)據(jù)特征的參數(shù)條件，決定待處理數(shù)據(jù)下一步指令的過程。四叉樹分割算法通過引入用于二維空間數(shù)據(jù)分析的四叉樹，對(duì)每次下達(dá)指令后滿足分割條件的局部影像進(jìn)行編碼。如果滿足分割參數(shù)要求，則影像的當(dāng)前窗口被分割為大小一致的4個(gè)方形，持續(xù)該步驟直至最小分割閾值。因此，為滿足后續(xù)數(shù)據(jù)處理要求，地理編碼后獲取的原始形變圖需經(jīng)裁剪、填充，使其大小滿足分割閾值的正整數(shù)倍。Jonsson[12]最先把四叉樹影像分割應(yīng)用于形變圖降采樣，將原始影像分割為4個(gè)象限，計(jì)算各象限均值，如果像元均值的均方根不滿足閾值要求，則把每個(gè)象限再次分割為4個(gè)新象限，繼續(xù)計(jì)算各象限像素均值的均方根并與閾值作比較，迭代該過程直至收斂。在實(shí)際采樣過程中，為提高算法運(yùn)算效率以及局部分割能力，一般直接將計(jì)算出的各象限最大梯度值作為參考條件，即把視線向形變量梯度值作為唯一參數(shù)，各形變圖窗口進(jìn)行自適應(yīng)分割。但該算法只關(guān)注了形變量的數(shù)值特征，并未考慮其自身數(shù)據(jù)質(zhì)量。因此，需要引入額外的描述數(shù)據(jù)質(zhì)量的特征參數(shù)，約束影像分割指令，如干涉影像的相干性系數(shù)，提高采樣點(diǎn)質(zhì)量以確保反演模型解算結(jié)果的可靠性。

1.2 相干系數(shù)計(jì)算與數(shù)據(jù)預(yù)處理

利用InSAR技術(shù)獲取高程模型和形變場(chǎng)的信息質(zhì)量取決于干涉相位質(zhì)量。干涉相位質(zhì)量由2幅單視復(fù)數(shù)(single look complex，SLC)影像對(duì)應(yīng)分辨單元內(nèi)總回波信號(hào)矢量值之間差異程度決定。地物反射特性隨時(shí)間的變化以及空間基線過長、衛(wèi)星軌道控制誤差引起的時(shí)空失相干等因素都會(huì)影響SAR影像之間的干涉質(zhì)量，一般用相干系數(shù)作為評(píng)價(jià)相干性的指標(biāo)。在實(shí)際數(shù)據(jù)處理過程中，利用式(1)所示的窗口估計(jì)方法求得歸一化相干系數(shù)近似值[13]。

(1)

SAR影像使用雷達(dá)坐標(biāo)系統(tǒng)，為便于后續(xù)數(shù)據(jù)分析，本文在計(jì)算主從SLC影像相干系數(shù)估計(jì)值后，基于軌道信息和成像參數(shù)將相干系數(shù)圖由雷達(dá)坐標(biāo)系地理編碼至WGS84坐標(biāo)系統(tǒng)。同時(shí)，如上文所述，依據(jù)四叉樹算法輸入數(shù)據(jù)大小要求和實(shí)際形變場(chǎng)分布特征，對(duì)所獲取實(shí)驗(yàn)區(qū)域的形變數(shù)據(jù)和相干質(zhì)量數(shù)據(jù)進(jìn)行裁剪，并進(jìn)行NaN值填充，以便于整體采樣。NaN值填充一般采用鄰近像素均值插值法，但對(duì)于存在NaN值較為集中分布的情況，如水體等低相干區(qū)域，均值窗口太大，會(huì)造成交界處梯度過大，導(dǎo)致局部四叉樹分割窗口中引入額外的虛假信息；而在區(qū)域面積大的情況，設(shè)置較小的均值窗口，則會(huì)降低算法效率。因此，本文提出一種適用于四叉樹采樣的擴(kuò)散填充方法，首先以四叉樹最小分割窗口大小建立初始窗口，從NaN值和有效值交界處進(jìn)行均值插值，向外擴(kuò)散插值的同時(shí)，窗口擴(kuò)大為前一個(gè)窗口的2倍，直至完全填充。NaN值分布主要有3種情況：小區(qū)域、內(nèi)部大區(qū)域和影像邊緣大區(qū)域。對(duì)于小區(qū)域的NaN值而言，利用均值計(jì)算基本不會(huì)引入誤差，因此可以不用建立掩膜文件。對(duì)于大區(qū)域的NaN值，主要是防止邊緣梯度過大，免于增加交界處有效信息區(qū)域多余采樣點(diǎn)，而對(duì)于區(qū)域內(nèi)部可生成掩膜文件，采樣結(jié)束后剔除該區(qū)域內(nèi)部的采樣點(diǎn)。圖1是對(duì)64像素×64像素的模擬數(shù)據(jù)以1×1的初始窗口得到的填充結(jié)果和掩膜文件，其中模擬數(shù)據(jù)中白色部分為NaN值區(qū)域。

注：掩膜圖像中黑色表示掩膜區(qū)域；白色表示保留區(qū)域。圖1 模擬數(shù)據(jù)擴(kuò)散填充結(jié)果和掩膜文件

1.3 基于相干系數(shù)的概率采樣模型

不同于近似規(guī)則分布的形變場(chǎng)數(shù)據(jù)，相干性系數(shù)整體關(guān)聯(lián)特征不明顯且受局部區(qū)域多因素影響，所以在設(shè)定整體相干閾值的同時(shí)，還需要設(shè)立一種描述局部窗口總相干特征的物理量，近似體現(xiàn)窗口數(shù)據(jù)質(zhì)量。因此，本文在給定相干系數(shù)閾值的同時(shí)，建立了描述窗口內(nèi)像素值接受程度的核函數(shù)K，相干系數(shù)值核運(yùn)算后獲得窗口整體相干表征量γk。該窗口相干表征量γk的表達(dá)如式(2)所示。

(2)

式中：2n表示當(dāng)前窗口大??；γr，c和Kr，c分別代表窗口內(nèi)r行c列相干系數(shù)值及其核元素值。

在當(dāng)前待分割窗口內(nèi)，對(duì)于滿足整體相干閾值γg的像素，其對(duì)應(yīng)核元素設(shè)為1，不滿足閾值要求設(shè)為k。以圖2中4×4矩陣為例，假定整體相干性閾值是0.5，創(chuàng)建對(duì)應(yīng)的核。

圖2 核創(chuàng)建示意圖

k值的大小決定了低相干像元在整體相干表征量中的權(quán)重，在相同大小窗口內(nèi)，k值越大，低質(zhì)量點(diǎn)位相干系數(shù)權(quán)重越高，γk的值越小，接受程度越低。圖3表示在大小為4×4窗口內(nèi)，低相干數(shù)目m分別是1、5、8、11和15時(shí)，低相干像元對(duì)整體表征量的影響程度關(guān)于k值的變化曲線。其中，4個(gè)*號(hào)點(diǎn)代表高、低質(zhì)量2種像元點(diǎn)對(duì)整體表征量影響程度相等時(shí)k值的取值。當(dāng)2種像元點(diǎn)數(shù)目相同時(shí)，k值取1，影響程度相同，隨著低相干點(diǎn)位數(shù)目的增加，k值左移。在實(shí)際情況下，k取值一般在1到2之間，以保證在提高低相干點(diǎn)位對(duì)采樣結(jié)果影響程度的同時(shí)，均衡高相干質(zhì)量點(diǎn)位權(quán)重，以免導(dǎo)致采樣點(diǎn)過少丟失重要細(xì)節(jié)信息。

圖3 低相干像元數(shù)目和k值對(duì)整體表征量影響曲線

本文算法中，形變梯度值和整體相干表征量同時(shí)作為四叉樹分割的約束參數(shù)，而整體相干表征量可在歸一化的相干系數(shù)值基礎(chǔ)上計(jì)算得到，并且核運(yùn)算結(jié)果并未改變其定量本質(zhì)，因此，只須對(duì)形變量做歸一化處理。經(jīng)歸一化處理后，更利于概率函數(shù)中變量權(quán)重因子的分配和調(diào)節(jié)，擴(kuò)大算法優(yōu)化空間。相較于常規(guī)的全局相干系數(shù)閾值降采樣方法而言，本文新算法提高整體閾值γg的同時(shí)，設(shè)定了一個(gè)較低的極限閾值γ′，一般設(shè)為0.1至0.15之間，以確保合理采樣較低相干點(diǎn)中的關(guān)鍵細(xì)節(jié)形變信息。當(dāng)待分割窗口形變梯度滿足閾值要求，而相干質(zhì)量小于整體閾值時(shí)，局部窗口的下一步指令由梯度值和相干表征量綜合決定。基于sigmoid函數(shù)單調(diào)、映射區(qū)間處于0到1之間、關(guān)于中心對(duì)稱等特點(diǎn)，經(jīng)適當(dāng)變換創(chuàng)建接受分割概率函數(shù)如式(3)所示。

(3)

式中：τ是根據(jù)sigmoid函數(shù)在自變量大于4后因變量變化趨于平緩，近似為1的性質(zhì)，作為壓縮常量，一般取4到8之間即可，默認(rèn)值為6；α和β分別是歸一化后的形變量和相干表征量的權(quán)重系數(shù)；d′和γ′是形變梯度閾值和相干系數(shù)極限閾值；x=α·d+β·γ，d和γ是當(dāng)前窗口的最大形變梯度和整體相干表征量。系數(shù)α和β的比值越大，則采樣更偏重于低質(zhì)量區(qū)域中的形變信息；反之，更偏重當(dāng)前窗口的相干性，可根據(jù)實(shí)際采樣要求進(jìn)行設(shè)置。分割概率曲線如圖4所示。

圖4 分割概率曲線

該算法從整體到局部，結(jié)合形變信息、相干信息對(duì)形變場(chǎng)進(jìn)行四叉樹分割，理論上可以有效保留細(xì)節(jié)形變特征和增加高質(zhì)量像元的采樣率，但一定程度上不能高效判別形變場(chǎng)附近山體疊掩引入的虛假形變信息。因此，考慮利用該地區(qū)的數(shù)字高程模型，設(shè)置合理的高程閾值并優(yōu)化上文中的掩膜文件，從而保證采樣點(diǎn)的可靠性。圖5為算法主要部分的流程圖。

圖5 基于相干系數(shù)的概率采樣算法流程圖

2 應(yīng)用實(shí)例

為了驗(yàn)證本文算法的實(shí)用性和可靠性，選取2016年2月6日中國臺(tái)灣地區(qū)美濃Mw6.7級(jí)地震作為實(shí)例分析，獲取2016年1月9日和2016年2月14日的一對(duì)歐空局Sentinel-1A數(shù)據(jù)進(jìn)行干涉處理，并引入30 m分辨率的SRTM DEM作為參考地形數(shù)據(jù)，使用GMTSAR軟件獲取升軌形變場(chǎng)。裁剪其中大小為4 096像素×4 096像素的形變區(qū)域進(jìn)行算法實(shí)驗(yàn)和結(jié)果評(píng)價(jià)。

為保證同時(shí)顧及形變場(chǎng)的形變特征和質(zhì)量特性，設(shè)置核元素k值為2，顧及形變特征的同時(shí)更注重像素值質(zhì)量。極限閾值γ′設(shè)為1.5，根據(jù)待分割影像窗口內(nèi)部像元質(zhì)量分布情況創(chuàng)建相應(yīng)的核，計(jì)算出相干表征量γk。由于實(shí)驗(yàn)區(qū)整體相干質(zhì)量較低，為提高采樣結(jié)果整體質(zhì)量，整體閾值γg設(shè)為0.7，采樣高質(zhì)量點(diǎn)位的同時(shí)，合理分配較可信點(diǎn)位的接受概率。

圖6給出了本文提出的優(yōu)化算法的采樣結(jié)果和常規(guī)四叉樹采樣結(jié)果的對(duì)比，為保證實(shí)驗(yàn)結(jié)果分析的合理性和公平性，2種采樣算法的形變特征分割指令都是基于形變場(chǎng)歸一化后，以0.25作為形變梯度閾值獲得。其中，原始形變圖中的R1、R2矩形框，L1、L2線段分別作為局部采樣結(jié)果對(duì)比分析區(qū)域和遠(yuǎn)距離剖線分析區(qū)域。為合理呈現(xiàn)整體對(duì)比結(jié)果，本文在圖6中給出算法最終采樣結(jié)果的同時(shí)，統(tǒng)計(jì)掩膜處理前2種采樣算法所獲取的點(diǎn)位在不同質(zhì)量區(qū)間的數(shù)目，繪制直方圖進(jìn)行直觀對(duì)比。圖6(e)是根據(jù)擴(kuò)散填充法生成的NaN值區(qū)域以及高程大于1 000 m的山地區(qū)域生成的掩膜文件。對(duì)比常規(guī)四叉樹和優(yōu)化算法的采樣點(diǎn)分布圖(圖6(c)和圖6(f))以及數(shù)字高程模型，形變場(chǎng)右上角和右下角采樣的虛假形變信息可以被很好地去除；同時(shí)，整體數(shù)據(jù)量由16 777 216個(gè)分別減少至5 815個(gè)(圖6(c))和2 853個(gè)數(shù)據(jù)點(diǎn)(圖6(f))。相對(duì)于經(jīng)典四叉樹采樣，優(yōu)化算法數(shù)據(jù)點(diǎn)數(shù)目下降了50.94%；掩膜之前數(shù)據(jù)量為4 285個(gè)，同比下降26.31%。在一定程度上，優(yōu)化算法在保持形變特征的同時(shí)，壓縮了采樣點(diǎn)數(shù)量，從數(shù)據(jù)源本身提高了后續(xù)反演的計(jì)算效率。并且，由優(yōu)化算法影像分割示意圖(圖6(d))和最終的優(yōu)化采樣結(jié)果(圖6(f))可知，本文算法能夠較好地保留形變場(chǎng)的整體形變信息。根據(jù)在不同相干系數(shù)區(qū)間采樣點(diǎn)數(shù)目的統(tǒng)計(jì)直方圖(圖7)可知，相干質(zhì)量越低的區(qū)間，采樣點(diǎn)數(shù)目大量減少，例如0～0.3相干區(qū)間采樣點(diǎn)數(shù)目降低了43.66%，而對(duì)于相干性系數(shù)大于0.5的采樣點(diǎn)幾乎都被接受。

圖7 掩膜前不同相干系數(shù)區(qū)間采樣點(diǎn)數(shù)目統(tǒng)計(jì)圖

圖8是R1、R2區(qū)域基于經(jīng)典四叉樹和本文算法的影像分割過程示意圖，以及相應(yīng)區(qū)域的相干性系數(shù)圖。由圖8第1行(R1)可知，在局部相干質(zhì)量較差時(shí)，本文算法在保證形變特征充分采樣的前提下，有效減少了低質(zhì)量點(diǎn)位的采樣數(shù)目；而對(duì)比第2行(R2)的結(jié)果，局部影像左下區(qū)域和右下區(qū)域相干性系數(shù)較高，采樣點(diǎn)數(shù)目并未減少。因此，本文算法可以在相干質(zhì)量分布具有明顯差異的局部區(qū)域進(jìn)行高效率自適應(yīng)采樣，從而提高整體采樣的可靠性。同時(shí)，本文給出2條主形變剖面的形變量變化曲線(圖9)，無論是整體變化特征還是局部細(xì)節(jié)特征都可以被完整地保留，體現(xiàn)了本文算法的可靠性和采樣優(yōu)勢(shì)。如表1所示，計(jì)算原始形變場(chǎng)、常規(guī)四叉樹及本文算法采樣結(jié)果的均方根值(root mean square，RMS)作為有效值評(píng)估參考，計(jì)算結(jié)果表明3種數(shù)據(jù)RMS大致相當(dāng)，說明有效信息采樣率很高。本文算法運(yùn)算時(shí)間比常規(guī)算法長，但在可接受范圍內(nèi)。

圖8 局部R1和R2區(qū)域采樣結(jié)果

注：剖線方向?yàn)閺淖笾劣?。圖9 局部L1、L2形變采樣結(jié)果

表1 采樣算法結(jié)果對(duì)比

3 結(jié)束語

本文提出了顧及InSAR相干系數(shù)的自適應(yīng)概率采樣算法，并給出數(shù)據(jù)降采樣處理的過程模型，即“影像擴(kuò)散填充—計(jì)算形變梯度和相干表征量—計(jì)算局部自適應(yīng)接受分割概率—掩膜”。通過對(duì)2016年中國臺(tái)灣地區(qū)美濃Mw6.7級(jí)同震形變場(chǎng)進(jìn)行采樣實(shí)驗(yàn)，充分證明了該算法的可靠性和有效性。

依據(jù)真實(shí)數(shù)據(jù)的定性定量分析，將形變場(chǎng)相干性系數(shù)作為影像自適應(yīng)分割的約束參數(shù),能夠同時(shí)從整體和局部上兼顧細(xì)節(jié)形變特征和數(shù)據(jù)質(zhì)量特征，彌補(bǔ)了目前形變場(chǎng)降采樣模型缺少考慮InSAR干涉質(zhì)量的缺陷，從數(shù)據(jù)源本身保證了后續(xù)反演流程的可靠性，有利于提高構(gòu)造形變模型反演的準(zhǔn)確性。