基于快速分布式目標(biāo)探測的時(shí)序雷達(dá)干涉測量方法:以Lost Hills油藏區(qū)為例

蔣彌， 丁曉利， 何秀鳳， 李志偉， 史國強(qiáng)

1 河海大學(xué)地球科學(xué)與工程學(xué)院， 南京　210098 2 香港理工大學(xué)土地測量與地理資訊學(xué)系， 香港九龍 3 中南大學(xué)地球科學(xué)與信息物理工程學(xué)院， 長沙　410083

?

基于快速分布式目標(biāo)探測的時(shí)序雷達(dá)干涉測量方法:以Lost Hills油藏區(qū)為例

蔣彌1， 丁曉利2， 何秀鳳1， 李志偉3， 史國強(qiáng)1

1 河海大學(xué)地球科學(xué)與工程學(xué)院， 南京210098 2 香港理工大學(xué)土地測量與地理資訊學(xué)系， 香港九龍 3 中南大學(xué)地球科學(xué)與信息物理工程學(xué)院， 長沙410083

針對當(dāng)前分布式目標(biāo)雷達(dá)干涉測量運(yùn)算效率低、選點(diǎn)困難等問題，本文提出了一種建立在快速同質(zhì)點(diǎn)選取下的干涉數(shù)據(jù)處理框架.相比之前的時(shí)序數(shù)據(jù)處理方法，新方法具有選點(diǎn)快速、自適應(yīng)性強(qiáng)的特點(diǎn)，能在保留影像分辨率基礎(chǔ)之上增加空間點(diǎn)密度.另外，在統(tǒng)計(jì)推斷的基礎(chǔ)上，提出基于無偏空間相干性估計(jì)的分布式目標(biāo)選擇方法，進(jìn)而彌補(bǔ)了傳統(tǒng)經(jīng)驗(yàn)閾值設(shè)定的缺陷.本文以美國加州Lost Hills油田區(qū)為例，在論證數(shù)據(jù)處理框架的可行性基礎(chǔ)之上，分析了因孔隙流體萃取和孔隙壓力降低引起的地表變形.

合成孔徑雷達(dá)干涉測量； 時(shí)序； 分布式目標(biāo)； 快速同質(zhì)點(diǎn)選取

1　引言

作為雷達(dá)遙感領(lǐng)域的重要分支，合成孔徑雷達(dá)干涉測量(InSAR)已廣泛應(yīng)用于地球科學(xué)領(lǐng)域，成為地表參數(shù)反演不可或缺的對地觀測手段之一(張紅等, 2009; 郭華東等, 2000).其中，InSAR的衍生技術(shù)，差分干涉測量(D-InSAR)，可用于監(jiān)測地表微小形變，理論上可以達(dá)到毫米級的監(jiān)測精度(Bamler and Hartl, 1998).雖然D-InSAR技術(shù)在理論上趨于成熟，也不乏成功的科學(xué)應(yīng)用案例，但在解決實(shí)際工程應(yīng)用問題中，要達(dá)到毫米級精度的形變檢測水平，該技術(shù)仍然需要解決大氣效應(yīng)和時(shí)空失相關(guān)等誤差源的影響(Li et al., 2006; 田馨和廖明生, 2013; 占文俊等, 2015; Zebker and Villasenor, 1992).

為了克服上述困難，在過去的數(shù)年中對于時(shí)間序列InSAR數(shù)據(jù)的點(diǎn)目標(biāo)之后向散射特性的研究成為實(shí)現(xiàn)地表微小形變檢測的突破口(Ferretti et al., 2001; 張紅等, 2009).在這些研究成果之中，永久散射體干涉測量技術(shù)(PS-InSAR)是行之有效的方法之一.它的基本原理是利用多景同一地區(qū)雷達(dá)影像，通過相位和幅度信息，查找不受時(shí)間、空間去相關(guān)影響的點(diǎn)目標(biāo)(例如，人工建筑物、裸露的巖石等).在統(tǒng)計(jì)分析基礎(chǔ)之上，PS-InSAR能夠減輕點(diǎn)目標(biāo)的大氣等誤差的影響，從而估算出每個(gè)PS點(diǎn)上的移動(dòng)速率或時(shí)序位移.然而，這一技術(shù)也存在一些限制.其中最主要的問題是在非人工地表區(qū)域，可用的永久散射體(Persistent Scatterer, PS)數(shù)目十分有限，導(dǎo)致觀測結(jié)果不能客觀反映研究區(qū)域的整體變化.從技術(shù)角度而言，低密度的空間點(diǎn)目標(biāo)在構(gòu)網(wǎng)之后可能給相位解纏和形變解算帶來困難.為提高點(diǎn)密度，改善時(shí)序形變估計(jì)的質(zhì)量，以小基線技術(shù)(SBAS-InSAR)和PS變種技術(shù)(如SqueeSAR，StaMPS)為代表的時(shí)序方法備受矚目(Berardino et al., 2002; Ferretti et al., 2011; Goel and Adam, 2012; Hooper et al., 2007; Wang et al., 2011).除了PS點(diǎn)之外，這些方法還探究分布式目標(biāo)點(diǎn)(distributed scatterer, DS)的相位信息.DS點(diǎn)目標(biāo)機(jī)制的特點(diǎn)涉及分辨率單元內(nèi)所有較小散射體的相干累加，這些散射體中沒有一個(gè)的散射特性是占據(jù)統(tǒng)治地位的(Lee and Pottier, 2009).雖然上述技術(shù)都是通過探索DS點(diǎn)(或PS與DS結(jié)合)為目標(biāo)，但是方法原理卻不同.例如，SBAS技術(shù)通過選取短時(shí)空基線像對減小時(shí)空去相關(guān)的影響.另外，這一方法還需要通過常規(guī)多視或多重空間濾波處理處理增加散射體的信噪比.因此，SBAS應(yīng)該是一種通過犧牲影像分辨率換取估計(jì)精度的方法.當(dāng)影像細(xì)節(jié)變得豐富或空間形變信號復(fù)雜時(shí)，該方法容易損失許多細(xì)節(jié)，并伴隨形變非形變區(qū)域的誤讀.相反地，以SqueeSAR為代表的技術(shù)仍然采用了單一主影像的策略(主要優(yōu)點(diǎn)是減輕主影像大氣誤差和估計(jì)純隨機(jī)噪聲)，通過分析時(shí)間棧樣本的相似程度，以假設(shè)檢驗(yàn)為手段(例如基于兩個(gè)樣本的Kolmogorov-Smirnov檢驗(yàn)(Ferretti et al., 2011)，自適應(yīng)檢驗(yàn)(Jiang et al., 2013)，Anderson-Darling檢驗(yàn)(Goel and Adam, 2012; Parizzi and Brcic, 2011))，選取與目標(biāo)像素具有相同統(tǒng)計(jì)屬性的點(diǎn)并進(jìn)行自適應(yīng)多視.這一方法的理論假設(shè)前提是同質(zhì)點(diǎn)都具有相同的散射特性和相位歷史.之后，再采用常規(guī)PS數(shù)據(jù)處理鏈，聯(lián)合PS和DS點(diǎn)進(jìn)行形變反演.SqueeSAR技術(shù)的優(yōu)點(diǎn)是保持空間分布特征，因而很好的彌補(bǔ)了SBAS技術(shù)的缺點(diǎn).但是，假設(shè)檢驗(yàn)步驟卻十分耗時(shí)，甚至無法應(yīng)用于大數(shù)據(jù)集.當(dāng)影像數(shù)較少時(shí)(一般要大于25景)，假設(shè)檢驗(yàn)功效也很低，這導(dǎo)致偏向零假設(shè)的情況出現(xiàn)，即錯(cuò)誤的接受兩個(gè)樣本是相似的(Jiang et al., 2013)結(jié)果，許多異質(zhì)點(diǎn)都包含在同質(zhì)點(diǎn)之中，造成平均后的影像分辨率的損失(Jiang et al., 2013).另外，目前的DS點(diǎn)選取方案多采用對相干性進(jìn)行經(jīng)驗(yàn)閾值設(shè)定.考慮到樣本相干性自身的隨機(jī)特征，使用經(jīng)驗(yàn)閾值易選取低質(zhì)量的DS點(diǎn)或刪除高質(zhì)量的點(diǎn)目標(biāo).

基于上述問題，本文將提出一種分布式目標(biāo)雷達(dá)干涉數(shù)據(jù)處理框架(命名為FaSHPS-InSAR).FaSHPS-InSAR技術(shù)首先采用快速同質(zhì)選點(diǎn)策略對每個(gè)DS目標(biāo)點(diǎn)進(jìn)行自適應(yīng)多視(Jiang et al., 2015)；然后，在空間相干性精確估計(jì)的條件下，采用無偏相干值和統(tǒng)計(jì)分析對DS點(diǎn)進(jìn)行進(jìn)一步篩選；最后，以穩(wěn)健參數(shù)估計(jì)為手段，恢復(fù)研究區(qū)形變歷史.為減小解纏誤差對時(shí)序變形的影響，將使用重新加權(quán)的最小二乘策略.根據(jù)FaSHPS-InSAR技術(shù)，本文以美國加州Lost Hills油田區(qū)為例，采用27景RadarSAT-1雷達(dá)數(shù)據(jù)對實(shí)驗(yàn)區(qū)進(jìn)行了地表變形反演，并與STAMPS PS-InSAR處理結(jié)果對比分析.研究結(jié)果揭示了該區(qū)域受油田開采影響的時(shí)序形變歷史，為反演因孔隙流體萃取和孔隙壓力降低引起的儲(chǔ)層變形提供數(shù)據(jù)支撐.

2　分布式目標(biāo)干涉數(shù)據(jù)處理核心算法

2.1FaSHPS算法

快速同質(zhì)點(diǎn)選取的核心是在高斯假設(shè)的背景下，將假設(shè)檢驗(yàn)問題轉(zhuǎn)換為置信區(qū)間估計(jì)的一種方法.由于假設(shè)檢驗(yàn)需要逐個(gè)比較參考點(diǎn)與候選點(diǎn)的時(shí)間樣本，計(jì)算P值與比較分布特征需要很大開銷(Jiang et al., 2015).而置信區(qū)間一旦確立，所有候選點(diǎn)只需通過簡單的邏輯運(yùn)算就可以判定與參考點(diǎn)的相似性，因而大大提高運(yùn)算效率.此外，由于模型確定，即便影像數(shù)較小，同質(zhì)點(diǎn)選擇產(chǎn)生的不確定性也小于傳統(tǒng)非參數(shù)假設(shè)檢驗(yàn)方法.這兩點(diǎn)成為這一技術(shù)的主要優(yōu)勢.

(1)

其中P{·}表示概率，z1-α/2為標(biāo)準(zhǔn)正態(tài)分布的1-α/2分位點(diǎn)，Var(A(p))表示p點(diǎn)的真實(shí)方差.根據(jù)分布式目標(biāo)SAR圖像統(tǒng)計(jì)理論，在勻質(zhì)區(qū)域，單視強(qiáng)度影像服從瑞利分布(Lee and Pottier, 2009).其數(shù)字特征，即變異系數(shù)CV是恒量：

(2)

這里E(·)表示期望，在本文表示為μ(p).當(dāng)后向散射特性在時(shí)間上變化不是十分劇烈的時(shí)候，可以假設(shè)時(shí)間上的樣本是勻質(zhì)的.因此，式(1)可以轉(zhuǎn)換為一個(gè)只包含了μ(p)的區(qū)間：

(3)

當(dāng)N較小時(shí)正態(tài)性很難維持，這直接影響著區(qū)間估計(jì)(3)式的有效性.根據(jù)經(jīng)典統(tǒng)計(jì)文獻(xiàn)(Papoulis, 1991)，當(dāng)概率分布平滑時(shí)(例如服從Rayleigh分布的單視強(qiáng)度序列)，N=5就足夠靠近高斯假設(shè).

2.2DS候選點(diǎn)選取

相干性是InSAR數(shù)據(jù)處理中最重要的參數(shù)之一，它貫穿了整個(gè)數(shù)據(jù)處理流程(Jiang et al., 2013; Jiang et al., 2014a).在點(diǎn)目標(biāo)選取中，常規(guī)方法是采用經(jīng)驗(yàn)閾值設(shè)定，即選取高于相干性閾值的那些點(diǎn)作為DS候選點(diǎn).由于干涉圖噪聲分布特征不同，每個(gè)帶估點(diǎn)真實(shí)相干性不同，則估計(jì)的樣本相干性偏差不同(Abdelfattah and Nicolas, 2006; Seymour and Cumming, 1994).因此，在低相干區(qū)域，設(shè)置較小的相干性閾值有可能小于高估系統(tǒng)偏差，使得那些原本干涉信號弱的點(diǎn)目標(biāo)被選取為DS點(diǎn).在高相干區(qū)域，高相干閾值又容易損失DS點(diǎn).根據(jù)每個(gè)像素的參數(shù)空間(視數(shù)，樣本相干性等)自適應(yīng)的選擇閾值判斷DS點(diǎn)便成為本文的主要目標(biāo).這主要由精確相干性估計(jì)和閾值估計(jì)兩個(gè)部分構(gòu)成.

首先相干性觀測可以由以下估計(jì)量獲得(Jiang et al., 2013; Touzi et al., 1999)：

(4)

(5)

一般地，考慮到計(jì)算效率和估計(jì)精度，R=200是一個(gè)折衷的選擇.

(6)

其中k是常數(shù)，代表倍數(shù)關(guān)系，在本文的研究中k=1；σγ表示無偏相干性的Cramer-Rao下界標(biāo)準(zhǔn)偏差(Touzi et al., 1999)：

(7)

3　FaSHPS-InSAR技術(shù)框架構(gòu)建

流程圖1為FaSHPS-InSAR的整體技術(shù)框架.其基本數(shù)據(jù)處理步驟如下：

假設(shè)所有時(shí)序單視復(fù)數(shù)影像(SLC)已被采樣至公共主影像，第一步，將所有影像取模獲得強(qiáng)度序列，采用差分振幅離差指數(shù)初步選擇候選點(diǎn)，再采用2.1節(jié)提出的FaSHPS算法為每個(gè)待估點(diǎn)選取同質(zhì)像元；

第二步，根據(jù)每個(gè)同質(zhì)點(diǎn)集合Ω進(jìn)行自適應(yīng)多視得到差分干涉圖序列IFGs；

第三步，根據(jù)2.2節(jié)的方法對初選點(diǎn)進(jìn)行提煉獲取DS候選點(diǎn)；

第五步，采用3D相位解纏算法解算濾波后的DS點(diǎn)；

圖1　FaSHPS-InSAR數(shù)據(jù)處理鏈流程圖(SLCs表示單視復(fù)數(shù)影像集合)Fig.1　The framework of FaSHPS-InSAR method

4　研究區(qū)概況與實(shí)驗(yàn)數(shù)據(jù)

Lost Hills油田位于美國加州圣華金河谷Kern郡.Kern郡以石油業(yè)和農(nóng)業(yè)為主要經(jīng)濟(jì)體，其石油儲(chǔ)備在加州位居第6，等價(jià)于11億桶原油產(chǎn)量(California Department of Conservation, 1991—2004).在圣華金河谷的局部區(qū)域，與油氣藏開采有關(guān)的儲(chǔ)層壓實(shí)與沉降作用已經(jīng)持續(xù)近60年，而LostHillls地區(qū)自20世紀(jì)50年代早期便觀測到地表沉降.隨著更多生產(chǎn)井的投入使用，該地區(qū)觀測到的沉降在80年代末期呈加速趨勢(Bruno and Bovberg, 1992).根據(jù)美國石油、天然氣和地?zé)豳Y源加州分部數(shù)據(jù)統(tǒng)計(jì)，該區(qū)域數(shù)十年來原油年產(chǎn)量穩(wěn)固提升，從1981年的近6×106桶(約950000 m3)增至2007年的1.2×107桶(約1940000 m3).

油氣藏壓實(shí)作用與沉降源于巖石孔隙壓力和有效應(yīng)力的孔隙塑性力學(xué)特性.當(dāng)石油等流體運(yùn)移出儲(chǔ)層時(shí)，上覆巖層的重量未減小，而孔隙壓力減小使得作用于固體基質(zhì)上的垂向有效應(yīng)力增加，導(dǎo)致壓實(shí).儲(chǔ)層壓實(shí)作用傳播到地面引起地表沉降.而壓實(shí)程度又取決于地層巖石屬性.以Lost Hills油田區(qū)為例，1995年冬季累計(jì)沉降達(dá)到15 cm(Xu, 2002)，沉降量遠(yuǎn)大于地理位置相鄰的且油產(chǎn)量更大的Belridge油區(qū).產(chǎn)生這一現(xiàn)象的主要原因歸咎于Lost Hills地區(qū)絕大部分屬硅藻巖油藏.相比Belridge更硬的砂巖，其儲(chǔ)層更淺(91 m)、生產(chǎn)層更厚(300 m)、滲透度更高(0.1～2 μm2)、黏度更低(0.096 Pa·s)(California Department of Conservation, 1991—2004).被壓實(shí)地層造成的井筒和設(shè)施損壞不僅無益于生產(chǎn)，而且對沿海地面的影響廣泛.

本文采用RadarSAT-1衛(wèi)星2002年2月至2004年2月獲取的27景SAR影像作為數(shù)據(jù)源.波長約為5.6 cm，中心入射角38.4°.圖2為實(shí)驗(yàn)區(qū)概況和成像覆蓋范圍.油藏區(qū)地處中央谷地平原的南部，地勢平坦，屬地中海型氣候.南部降雨少，較干燥，這里本是半沙漠地帶.加州水道系統(tǒng)的興建，此處成為農(nóng)業(yè)的重要產(chǎn)地.在西北至東南沿線，主要地表類型為人工建筑(其中以生產(chǎn)井和平臺(tái)為主)和裸土.沿線兩側(cè)分布著野生植被和農(nóng)作物，F(xiàn)aSHPS-InSAR和STAMPS PS-InSAR技術(shù)將同時(shí)用于獲取該地區(qū)的形變歷史.

5　結(jié)果分析

5.1干涉組合與數(shù)據(jù)處理

兩種數(shù)據(jù)處理方法的干涉組合策略不同，PS-InSAR方法是以一景影像為主影像，其他影像為輔影像形成的星狀圖，如圖3b所示.其中每個(gè)點(diǎn)代表獲取的單視復(fù)數(shù)影像(SLC)，線段表示干涉圖.FaSHPS-InSAR則采用多基線策略，盡可能地形成相干性高的干涉組合以便增加觀測的質(zhì)量和數(shù)量,如圖3a所示.在連接SLC過程中，應(yīng)確保沒有子集出現(xiàn)以便增加解的可靠性.在本文的研究中，兩種方法均配準(zhǔn)至2003年4月4日的影像參考幾何.

在數(shù)據(jù)配準(zhǔn)之后根據(jù)干涉組合和外部SRTM DEM在單視下形成差分干涉圖.PS-InSAR以傳統(tǒng)數(shù)據(jù)處理鏈進(jìn)行處理(Ferretti et al., 2001; Hooper et al., 2007)，此處不再贅述.其中振幅離差指數(shù)初始值設(shè)置為0.3.其他可用參數(shù)均與FaSHPS-InSAR保持一致(詳見表1).在FaSHPS-InSAR處理中，首先采用FaSHPS算法對像元逐個(gè)選取同質(zhì)點(diǎn)并確保點(diǎn)的連通性.本例中原始窗口大小為19×13.圖4展示了FaSHPS方法和傳統(tǒng)規(guī)則窗選點(diǎn)的比較結(jié)果.對于質(zhì)地顯著的區(qū)域，F(xiàn)aSHPS能夠避免與中心像素反射強(qiáng)度不同的散射體.

接著，根據(jù)選取的同質(zhì)點(diǎn)進(jìn)行相干性估計(jì).在相干性偏差糾正計(jì)算中，Bootstrap復(fù)制設(shè)置為R=200.圖5給出了其中一景干涉對的相干性估計(jì)結(jié)果，并與傳統(tǒng)7×7規(guī)則窗的估計(jì)結(jié)果進(jìn)行對照.整體上由于偏差存在，傳統(tǒng)估計(jì)的對比度差在快速去相干的農(nóng)業(yè)用地區(qū)域出現(xiàn)明顯的高估.即便采用的是7×7窗口，空間分辨率仍然較低.從圖5橢圓標(biāo)識(shí)放大的子區(qū)域，可以看出兩種算法的差異.本文方法估計(jì)的相干性在田埂處仍然保證了分辨率，與周邊農(nóng)田的低相干性形成對比.

圖3　不同時(shí)序分析技術(shù)干涉組合策略圖(a) FaSHPS-InSAR干涉組合，色帶表示干涉圖DS候選點(diǎn)平均空間相干性； (b) PS-InSAR技術(shù)干涉組合圖.Fig.3　The configuration of different time series InSAR techniques(a) FaSHPS-InSAR, the colorbar denotes the mean coherence of DS candidates in a interferogram; (b) PS.

圖4　FaSHPS算法與傳統(tǒng)規(guī)則窗口同質(zhì)點(diǎn)選擇方法(紅色點(diǎn)表示參考中心像元，綠色表示同質(zhì)點(diǎn))Fig.4　Homogeneous pixel selection using FaSHPS and conventional boxcar; the reference pixel is shown in red, and its brother neighbors are shown in green

圖5　FaSHPS算法與傳統(tǒng)規(guī)則窗口相干性估計(jì)結(jié)果對比圖(紅色點(diǎn)表示圖4中參考像元位置)Fig.5　The coherence estimation using FaSHPS and conventional boxcar respectively(the red point corresponds to the reference pixel in Fig.4)

圖6展示了基于常規(guī)最小二乘(L2)和ILRS(L1)的勻速率估計(jì)圖.選擇區(qū)域?yàn)閳D7中左下角的農(nóng)田區(qū)域.對于有解纏誤差存在的干涉序列(解纏圖(a)和(b)存在明顯的2pi跳變)，常規(guī)估計(jì)方法(圖6d)中出現(xiàn)了空間非相關(guān)性誤差并伴隨與解纏誤差相關(guān)的勻速率誤差，這在(圖6c)中并不顯著.

5.2Lost Hills油田區(qū)地表沉降監(jiān)測結(jié)果

5.2.1形變監(jiān)測交互驗(yàn)證

在FaSHPS-InSAR和PS-InSAR分別對DS和PS目標(biāo)點(diǎn)進(jìn)行解算之后，可以獲得各自的年平均形變速率圖.從圖7的結(jié)果來看，兩種方法都獲得了相似的沉降輪廓與形變分布.然而，F(xiàn)aSHPS-InSAR估計(jì)的勻速率在空間密度上有明顯優(yōu)勢，例如位于影像中心區(qū)域的沉降盆以及左邊區(qū)域的散射目標(biāo).在自然地表區(qū)域，缺少人工目標(biāo)或強(qiáng)反射散射體成為PS-InSAR低空間密度的主要原因.在增加點(diǎn)密度的情況下最小化分辨率的損失成為FaSHS-InSAR的主要優(yōu)點(diǎn).本例中，在120 km2范圍內(nèi)得到的DS點(diǎn)為148508個(gè)，PS點(diǎn)為10911個(gè).

表1　用于Lost Hills數(shù)據(jù)集的FaSHPS-InSAR數(shù)據(jù)處理參數(shù)Table 1　The parameters used in FaSHPS-InSAR algorithm over Lost Hills area

圖8展示了FaSHPS-InSAR和PS-InSAR處理同名點(diǎn)的時(shí)間序列位移情況.從結(jié)果來看，在相干性差的區(qū)域(如圖8(a，b))，PS點(diǎn)的強(qiáng)度保持穩(wěn)定，但是相位觀測仍具有不確定性，即隨機(jī)相位噪聲.考慮到PS點(diǎn)周圍的點(diǎn)密度更為稀疏，在STAMPS濾波時(shí)很難找到鄰域樣本進(jìn)行平滑.在執(zhí)行3D相位解纏時(shí)，存在的噪聲與稀疏的點(diǎn)密度使得傳統(tǒng)方法無法正確恢復(fù)每個(gè)干涉圖相位的整周數(shù).反演的時(shí)序離散度大.特別是點(diǎn)“2”，該目標(biāo)已經(jīng)遠(yuǎn)離生產(chǎn)井(見圖9)，位于農(nóng)田區(qū)域，持續(xù)的變形和振幅說明誤差的存在.相比之下，在FaSHPS-InSAR處理中，由于引入了更多的DS點(diǎn)，空間平均和增強(qiáng)的點(diǎn)密度可以改善時(shí)空平滑度，從而使得解纏錯(cuò)誤的概率減小.對于影像中心區(qū)域的點(diǎn)目標(biāo)“3”和“4”(圖8(c，d))，由于生產(chǎn)井增加了PS點(diǎn)的密度和點(diǎn)的質(zhì)量，兩種方法因而獲得相近的結(jié)果，這可以從兩個(gè)序列的均方根(RMS)得以驗(yàn)證.

圖6　解纏誤差對勻速率估計(jì)的影響(a)和(b)表示包含了解纏誤差的相位圖； (c) 基于ILRS估計(jì)的勻速率圖； (d) 基于常規(guī)最小二乘估計(jì)的勻速率圖.Fig.6　The impact of unwrapping errors on the deformation rate estimation(a) and (b) show the interferograms including unwrapping errors; Deformation rate estimated by ILRS (a) and conventional least square method (d).

圖7　DS目標(biāo)點(diǎn)和PS目標(biāo)點(diǎn)上獲取的年平均沉降量Fig.7　The deformation rate obtained from (a) DS and (b) PS targets

圖8　FaSHPS-InSAR和PS-InSAR同名點(diǎn)時(shí)間序列圖1,2,3,4點(diǎn)位置對應(yīng)于圖9；4個(gè)目標(biāo)點(diǎn)形變RMS值分別為10.56、20.93、1.95、3.12 mmFig.8　The pattern of displacement time series from FaSHPS-InSAR and PS-InSARThe positions of the points 1,2,3,4 are shown in Fig.9; the RMS value for four selected points is 10.56 mm,20.93 mm,1.95 mm,3.12 mm respectively.

最后，圖10展示了由FaSHPS-InSAR技術(shù)估計(jì)的2002—2004年油田區(qū)累計(jì)沉降二維時(shí)間序列，其中最大沉降量>16 cm.

圖9　FaSHPS-InSAR的DS目標(biāo)點(diǎn)殘差標(biāo)準(zhǔn)偏差圖Fig.9　The STD of the residuals of DS points

5.2.2沉降結(jié)果解譯

以上交互驗(yàn)證能夠說明FaSHPS-InSAR在增加點(diǎn)密度和保持處理精度上的能力.由于缺少外部實(shí)測GPS和水準(zhǔn)測量數(shù)據(jù)，驗(yàn)證工作主要依靠現(xiàn)有的地質(zhì)資料和已取得的研究成果.圖11給出了美國Chevron公司公布的Lost Hills地區(qū)油井分布圖(Xu, 2002).從中可以看出井的空間分布與圖7獲取的沉降位置十分吻合，即形變從西北向東南方向延伸，兩個(gè)沉降中心均與生產(chǎn)井密度最高的位置相對應(yīng).

由孔隙流體萃取和孔隙壓力降低引起的地表變形是一個(gè)非常復(fù)雜的孔隙彈塑性過程.若要深入理解油藏區(qū)沉降還需要了解壓實(shí)地層和巖石屬性，以及油氣田生產(chǎn)資料.一般地，沉降盆比壓實(shí)區(qū)域范圍廣，傳播幅度取決于上覆巖層的物理屬性和壓實(shí)地層的深度.在沉降盆中心，很少涉及水平運(yùn)動(dòng)，主要是垂向運(yùn)動(dòng).盆中心的視線向形變被投影至垂直向.

圖10　基于FaSHPS-InSAR的DS目標(biāo)點(diǎn)二維時(shí)間序列為更清晰地顯示沉降盆構(gòu)架，本文將27景序列進(jìn)行時(shí)間間隔采樣至14景沉降圖.Fig.10　The displacement time series estimated from FaSHPS-InSARFor exhibition purpose, the series are under-sampled from 27 to 14.

文獻(xiàn)(Xu,2002)中獲得的Lost Hills區(qū)域1995年冬季的最大沉降量為15 cm/105 days，等價(jià)于51.14 cm·a-1.相比于本文最大勻速率21.23 cm·a-1，形變速率的減小應(yīng)該與生產(chǎn)歷史相關(guān).為此我們統(tǒng)計(jì)了油年產(chǎn)量、年注水量、生產(chǎn)井?dāng)?shù)量和沉降率的關(guān)系.從圖12中可以看出，油產(chǎn)量在1995年前后有較大波動(dòng)，1994年的產(chǎn)量高達(dá)近1.4×107桶，而注水量僅為1.6×107桶，逆差2×106桶.在1992年達(dá)到注采平衡之前，注水量小于開采量.這些成為高沉降率的主要原因.1995年之后注水量仍在提升.特別是1996—1997年之間，注水量增加了近一倍，而油產(chǎn)量降低了5×106桶.持續(xù)不斷的注水緩解了沉降率，我們推斷這是獲得較低觀測結(jié)果的主要原因.值得注意的是，注水作用不會(huì)立刻發(fā)生.位于北海的Ekofisk油田就是典型的案例之一.根據(jù)測量歷史和巖石屬性，沉降速率的減小在注采平衡后4—5年才發(fā)生(Nagel, 2001; Thomas et al., 1987).這也間接說明了在1997年大幅度注水對后續(xù)緩解形變的影響是顯著的.

圖11　Lost Hills 油田區(qū)生產(chǎn)井分布圖(Xu, 2002)(其中黑色點(diǎn)狀目標(biāo)代表生產(chǎn)井)Fig.11　The distribution of producing well shown in black over Lost Hills area (Xu, 2002).

6　結(jié)論

本文提出了一種基于快速分布式目標(biāo)探測的時(shí)序雷達(dá)干涉數(shù)據(jù)處理方法.該方法由兩個(gè)核心部分構(gòu)成： (1) 采用快速選點(diǎn)方案對初始點(diǎn)進(jìn)行同質(zhì)點(diǎn)選?。?(2) 以真實(shí)相干性為閾值，根據(jù)統(tǒng)計(jì)推斷結(jié)果自適應(yīng)獲得樣本相干性閾值；另外，數(shù)據(jù)處理流程還融合了最新的濾波方法和穩(wěn)健估計(jì)手段.在構(gòu)建的FaSHPS-InSAR技術(shù)框架之下，本文采用27景美國加州Lost Hills油田區(qū)數(shù)據(jù)反演了當(dāng)?shù)氐牡乇沓两登闆r，并與PS-InSAR進(jìn)行了對比分析.實(shí)例結(jié)果表明，F(xiàn)aSHPS-InSAR在抗差、保留影像分辨率和增加空間觀測密度上具有明顯優(yōu)勢，這將有利于非人工地表的對地監(jiān)測.根據(jù)研究結(jié)果和當(dāng)?shù)厣a(chǎn)歷史，本文還揭示了油藏生產(chǎn)井密度與地表形變的關(guān)系.相比1995年的高沉降率，2002—2004年緩解的沉降量應(yīng)歸因于持續(xù)的回灌水量增加.然而沉降還在繼續(xù)，表明壓實(shí)作用依然存在.這將對地表基礎(chǔ)設(shè)施如井故障或套管損壞產(chǎn)生深遠(yuǎn)的影響.

除了上述特點(diǎn)之外，運(yùn)算效率是需要重點(diǎn)關(guān)注的問題.就同質(zhì)點(diǎn)選擇而言，對于80萬個(gè)初始點(diǎn)目標(biāo)，F(xiàn)aSHPS-InSAR方法在Linux MATLAB 2011b

圖12　Lost Hills 油田的生產(chǎn)歷史左坐標(biāo)軸表示油產(chǎn)量、注水量以及生產(chǎn)井?dāng)?shù)量；右坐標(biāo)軸表示沉降速率；柱狀圖表示油產(chǎn)量和井?dāng)?shù)量，藍(lán)色曲線表示注水情況，墨綠色的點(diǎn)和虛線表示沉降率變化.Fig.12　The history of oil production over Lost HillsThe axis in left denotes the oil production, injection (waterflood) and the number of producing well; the bars show the oil and producing well; the curve in blue denotes the injection; the points and the dot line in green denote the change of the subsidence rate.

Intel i7 2.80 GHz CPU 平臺(tái)上的運(yùn)算效率僅為約5 min.對于相同數(shù)量的點(diǎn)，SqueeSAR耗時(shí)將超過7h.這表明FaSHPS-InSAR更適用于處理大數(shù)據(jù)集和寬幅影像序列.因此，高效運(yùn)算成為FaSHPS-InSAR的主要技術(shù)優(yōu)勢.

Abdelfattah R, Nicolas J M. 2006. Interferometric SAR coherence magnitude estimation using second kind statistics.IEEETransactionsonGeoscienceandRemoteSensing, 44(7): 1942-1953. Bamler R, Hartl P. 1998. Synthetic aperture radar interferometry. Inverse Problems, 14(4): R1.

Berardino P, Fornaro G, Lanari R, et al. 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms.IEEETransactionsonGeoscienceandRemoteSensing, 40(11): 2375-2383.

Bruno M, Bovberg C A. 1992. Reservoir compaction and surface subsidence above the Lost Hills Field, California. ∥The 33th US Symposium on Rock Mechanics (USRMS). Santa Fe, New Mexico, USA: American Rock Mechanics Association.

California Department of Conservation. 1991—2004. Oil and gas statistics, annual report. California: California Department of oil, Gas and Geothermal Resources.

Efron B, Tibshirani R J. 1994. An Introduction to the Bootstrap. New York: CRC Press.

Ferretti A, Prati C, Rocca F. 2001. Permanent scatterers in SAR interferometry.IEEETransactionsonGeoscienceandRemoteSensing, 39(1): 8-20.

Ferretti A, Fumagalli A, Novali F, et al. 2011. A new algorithm for processing interferometric data-stacks: SqueeSAR.IEEETransactionsonGeoscienceandRemoteSensing, 49(9): 3460-3470. Goel K, Adam N. 2012. An advanced algorithm for deformation estimation in non-urban areas.ISPRSJournalofPhotogrammetryandRemoteSensing, 73: 100-110. Guo H D, Shao Y, Wang C L. 2000. Theory and Application of Radar Earth Observation (in Chinese). Beijing: Science Press.

Hooper A, Segall P, Zebker H. 2007. Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos.JournalofGeophysicalResearch:SolidEarth(1978—2012), 112(B7), doi: 10.1029/2006JB004763.

Jiang M, Ding X L, Li Z W. 2014. Hybrid approach for unbiased coherence estimation for multitemporal InSAR.IEEETransactionsonGeoscienceandRemoteSensing, 52(5): 2459-2473. Jiang M, Ding X L, Li Z W, et al. 2014a. The improvement for baran phase filter derived from unbiased InSAR coherence.IEEEJournalofSelectedTopicsinAppliedEarthObservationsandRemoteSensing, 7(7): 3002-3010. Jiang M, Ding X L, Tian X, et al. 2014b. A hybrid method for optimization of the adaptive Goldstein filter.ISPRSJournalofPhotogrammetryandRemoteSensing, 98: 29-43.

Jiang M, Ding X L, Hanssen R F, et al. 2015. Fast statistically homogeneous pixel selection for covariance matrix estimation for multitemporal InSAR.IEEETransactionsonGeoscienceandRemoteSensing, 53(3): 1213-1224.

Lee J S, Pottier E. 2009. Polarimetric Radar Imaging: From Basics to Applications. Boca Raton: CRC Press.

Li Z H, Fielding E J, Cross P, et al. 2006. Interferometric synthetic aperture radar atmospheric correction: GPS topography-dependent turbulence model.JournalofGeophysicalResearch:SolidEarth(1978—2012), 111(B2), doi: 10.1029/2005JB003711. Nagel N B. 2001. Compaction and subsidence issues within the petroleum industry: from wilmington to ekofisk and beyond.PhysicsandChemistryoftheEarth,PartA:SolidEarthandGeodesy, 26(1-2): 3-14.

Papoulis A. 1991. Probability, Random Variables, and Stochastic Processes. 3rd ed. New York: McGraw-Hill Book Company.

Parizzi A, Brcic R. 2011. Adaptive InSAR stack multilooking exploiting amplitude statistics: a comparison between different techniques and practical results.IEEEGeoscienceandRemoteSensingLetters, 8(3): 441-445.

Seymour M S, Cumming I G. 1994. Maximum likelihood estimation for SAR interferometry. ∥Proceedings of the Surface and Atmospheric Remote Sensing: Technologies, Data Analysis and Interpretation, International Geoscience and Remote Sensing Symposium. Pasadena: IEEE, 2272-2275.

Thomas L K, Dixon T N, Evans C E, et al. 1987. Ekofisk waterflood pilot.JournalofPetroleumTechnology, 39(2): 221-232.

Tian X, Liao M S. 2013. The analysis of conditions for InSAR in the field of deformation monitoring.ChineseJ.Geophys. (in Chinese), 56(3): 812-823, doi: 10.6038/cjg20130310.

Touzi R, Lopes A, Bruniquel J, et al. 1999. Coherence estimation for SAR imagery.IEEETransactionsonGeoscienceandRemoteSensing, 37(1): 135-149.

Wang T, Perissin D, Rocca F, et al. 2011. Three Gorges Dam stability monitoring with time-series InSAR image analysis.ScienceChinaEarthSciences, 54(5): 720-732.

Xu H B. 2002. Production induced reservoir compaction and surface subsidence, with applications to four-dimensional seismic[Ph. D. thesis]. Stanford: Stanford University.

Zebker H A, Villasenor J. 1992. Decorrelation in interferometric radar echoes.IEEETransactionsonGeoscienceandRemoteSensing, 30(5): 950-959.

Zebker H A, Chen K. 2005. Accurate estimation of correlation in InSAR observations.IEEEGeoscienceandRemoteSensingLetters, 2(2): 124-127.

Zhan W J, Li Z W, Wei J C, et al. 2015. A strategy for modeling and estimating atmospheric phase of SAR interferogram.ChineseJournalofGeophysics(in Chinese), 58(7): 2320-2329, doi: 10.6038/cjg20150710.

Zhang H, Wang C, Wu T, et al. 2009. The Study of DINSAR Technique Based on the Coherent Targets (in Chinese). Beijing: Science Press.

附中文參考文獻(xiàn)

郭華東, 邵蕓, 王長林. 2000. 雷達(dá)對地觀測理論與應(yīng)用. 北京: 科學(xué)出版社.

田馨, 廖明生. 2013. InSAR技術(shù)在監(jiān)測形變中的干涉條件分析. 地球物理學(xué)報(bào), 56(3): 812-823, doi: 10.6038/cjg20130310.

占文俊, 李志偉, 韋建超等. 2015. 一種InSAR大氣相位建模與估計(jì)方法. 地球物理學(xué)報(bào), 58(7): 2320-2329, doi: 10.6038/cjg20150710.

張紅, 王超, 吳濤等. 2009. 基于相干目標(biāo)的DInSAR方法研究. 北京: 科學(xué)出版社.

(本文編輯汪海英)

FaSHPS-InSAR technique for distributed scatterers: A case study over the lost hills oil field, California

JIANG Mi1, DING Xiao-Li2, HE Xiu-Feng1, LI Zhi-Wei3, SHI Guo-Qiang1

1SchoolofEarthScienceandEngineering,HohaiUniversity,Nanjing210098,China2Dept.ofLandSurveying&Geo-Informatics,TheHongKongPolytechnicUniversity,HongKong,China3SchoolofGeoscienceandInfo-PhysicsandGeomaticsEngineering,CentralSouthUniversity,Changsha410083,China

Multitemporal Synthetic Aperture Radar Interferometry (MT-InSAR) technique for the distributed scatterers (DSs) has been widely developed since 2010. The main limitations of this technique in the current state are high computing burden and the difficulty to define DSs with a suitable threshold. In this paper, we present a novel framework of data processing of radar interferometry for DSs. Compared with the state-of-the-art, the significant advantage of the new method is its low computing burden with high self-adaptability, and improved DSs density without loss of resolution. Besides, the DS candidates are selected according to the statistics of spatial sample coherence magnitude rather than the empirical threshold. The real SAR dataset over Lost Hillsoil field, California, is used to demonstrate the value of our method.

Interferometric synthetic aperture radar (InSAR); Time series; Distributed scatterer; Fast homogeneous statistically pixels selection

10.6038/cjg20161007.

國家自然科學(xué)基金(41404009)，江蘇省交通運(yùn)輸科技與成果轉(zhuǎn)化項(xiàng)目(2016)(16Y08)和中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(2013/B15020007)資助.

蔣彌，男，1982年生，博士，主要從事InSAR研究.E-mail：mijiang@hhu.edu.cn

10.6038/cjg20161007

P223,P631

2015-12-05，2016-08-23收修定稿

蔣彌， 丁曉利， 何秀鳳等. 2016. 基于快速分布式目標(biāo)探測的時(shí)序雷達(dá)干涉測量方法:以Lost Hills油藏區(qū)為例. 地球物理學(xué)報(bào)，59(10):3592-3603,

Jiang M, Ding X L, He X F, et al. 2016. FaSHPS-InSAR technique for distributed scatterers: A case study over the lost hills oil field, California.ChineseJ.Geophys. (in Chinese),59(10):3592-3603,doi:10.6038/cjg20161007.