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

(1.中國地震局地震預測研究所 地震預測重點實驗室，北京 100036；2.中國地震局地質(zhì)研究所 地震動力學國家重點實驗室，北京 100029)

0 引言

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

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

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

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

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

1.1 四叉樹影像分割

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

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

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

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

(1)

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

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

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

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

(2)

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

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

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

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

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

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

(3)

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

圖4 分割概率曲線

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

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

2 應用實例

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

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

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

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

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

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

注：剖線方向為從左至右。圖9 局部L1、L2形變采樣結(jié)果

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

3 結(jié)束語

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

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