金重陽(yáng) 劉 毓 鄧云開(kāi) 田衛(wèi)明，3 胡政權(quán)

（1.重慶三峽學(xué)院電子與信息工程學(xué)院，重慶 404130；2.北京理工大學(xué)信息與電子學(xué)院，北京 100081；3.北京理工大學(xué)重慶創(chuàng)新中心，重慶 401120）

1 引言

GB-InSAR 差分干涉處理分析形變時(shí)，需要以高質(zhì)量像素點(diǎn)作為研究對(duì)象，即PS 點(diǎn)。PS 點(diǎn)干涉相位主要包含形變相位、大氣相位及噪聲相位［1］。雷達(dá)發(fā)射的電磁波在大氣傳輸過(guò)程中，均會(huì)產(chǎn)生一定程度的折射，不同時(shí)刻所對(duì)應(yīng)的大氣條件發(fā)生變化，改變?cè)械膫鞑シ较蚺c路徑，產(chǎn)生大氣相位延遲［2］。理論數(shù)據(jù)表明在溫度為25 ℃，斜距1 km 處，僅1%的濕度變化可引起約1.3 mm形變測(cè)量誤差。大氣相位為GB-InSAR形變測(cè)量的主要誤差源［3-4］。

針對(duì)大氣相位的補(bǔ)償，目前主要有GCP（Ground Control Points，地面控制點(diǎn)）補(bǔ)償［5］、氣象數(shù)據(jù)校正［6］，以及參數(shù)模型法［7］三種方法。第1 種方法主要通過(guò)人工在場(chǎng)景內(nèi)布設(shè)角反射器作為GCP，通過(guò)估計(jì)GCP的大氣相位實(shí)現(xiàn)監(jiān)測(cè)目標(biāo)點(diǎn)的相位補(bǔ)償。在大氣勻質(zhì)性假設(shè)條件下，一般采用距離線性模型擬合各PS點(diǎn)的大氣相位。對(duì)于較為復(fù)雜的監(jiān)測(cè)場(chǎng)景，大氣不滿足勻質(zhì)性假設(shè)時(shí)，采用高階斜距模型或空間維插值實(shí)現(xiàn)區(qū)域內(nèi)PS點(diǎn)的大氣相位估計(jì)。2015年，汪學(xué)琴等學(xué)者［8］以不同分布的GCP剔除了目標(biāo)點(diǎn)的大部分大氣相位，與正垂線所得數(shù)據(jù)進(jìn)行了比對(duì)分析，驗(yàn)證了該方法的有效性。2017 年，黃其歡等學(xué)者［9］利用GCP 推算出的大氣擾動(dòng)因子，根據(jù)距離加權(quán)獲取所有待測(cè)點(diǎn)大氣擾動(dòng)相位，實(shí)現(xiàn)了大氣相位補(bǔ)償?？紤]到傳統(tǒng)的GCP補(bǔ)償需要人工布設(shè)角反射點(diǎn)，依賴于布設(shè)點(diǎn)的數(shù)量與質(zhì)量，常常面臨場(chǎng)景中存在危險(xiǎn)區(qū)域難以布設(shè)反射點(diǎn)的情況。第2種方法需要在觀測(cè)場(chǎng)景內(nèi)建立氣象站，獲得雷達(dá)數(shù)據(jù)對(duì)應(yīng)時(shí)刻的氣象參數(shù)（溫度、濕度、大氣壓），基于大氣折射率模型，定量求解大氣相位，實(shí)現(xiàn)相位補(bǔ)償。2008 年，Pipia 等［10］分析了大氣對(duì)電磁波傳輸?shù)挠绊?，證實(shí)了氣象數(shù)據(jù)的變化與干涉相位有直接關(guān)系。2014年，董杰等學(xué)者［6］使用角反射器作為監(jiān)測(cè)目標(biāo)，利用觀測(cè)區(qū)域內(nèi)的氣象數(shù)據(jù)構(gòu)建了大氣折射模型，定量求解了大氣擾動(dòng)相位，修正了監(jiān)測(cè)目標(biāo)的干擾相位。氣象數(shù)據(jù)校正依賴于監(jiān)測(cè)區(qū)域中氣象站的分布，數(shù)量較少的氣象站難以準(zhǔn)確獲得各區(qū)域的氣象數(shù)據(jù)，補(bǔ)償效果較差。第3 種方法則是對(duì)PS 點(diǎn)的大氣干涉相位數(shù)學(xué)建模，利用最小二乘回歸估計(jì)出模型參數(shù)，實(shí)現(xiàn)相位誤差的補(bǔ)償。2005年，Noferini等［7］首次提出了基于PS技術(shù)，估計(jì)大氣相位模型參數(shù)的方法。2014年，Iglesias等［11］分析了高程對(duì)于氣象條件的影響，將高程引入大氣相位參數(shù)模型，提高了測(cè)量精度。2016 年，徐亞明等［12］將該方法有效地應(yīng)用于施工期高危邊坡的監(jiān)測(cè)中。當(dāng)雷達(dá)監(jiān)測(cè)區(qū)域過(guò)大或天氣惡劣時(shí)，通常表現(xiàn)為較強(qiáng)的時(shí)空變性，部分干涉圖無(wú)法建立多參數(shù)模型來(lái)有效估計(jì)大氣相位，參數(shù)模型法誤差較大，不再適用［13］。

因此，本文提出了基于PS 技術(shù)的改進(jìn)方法，首先分析各鄰近PS 對(duì)的差分相位序列標(biāo)準(zhǔn)差，實(shí)現(xiàn)PS 集合中噪聲點(diǎn)以及形變點(diǎn)濾除，對(duì)剩余PS 采用MIC 結(jié)合空間維相干系數(shù)的雙閾值實(shí)現(xiàn)穩(wěn)定PS 的選取。利用穩(wěn)定PS空間非均勻分布的特點(diǎn)，通過(guò)空間維插值擬合大氣相位，從而實(shí)現(xiàn)大氣相位的有效補(bǔ)償。采用直線式掃描邊坡監(jiān)測(cè)雷達(dá)作為實(shí)驗(yàn)平臺(tái)，對(duì)30 幅雷達(dá)圖像進(jìn)行分析，對(duì)比驗(yàn)證了本文方法在復(fù)雜大氣條件下相位補(bǔ)償?shù)挠行浴?/p>

2 補(bǔ)償方法

2.1 參數(shù)模型補(bǔ)償法

該方法基于PS 技術(shù)，利用雷達(dá)圖像中PS 點(diǎn)估計(jì)模型參數(shù)。目前對(duì)于PS 點(diǎn)的選取廣泛采用幅度離差法，通過(guò)評(píng)估圖像中像素點(diǎn)的幅度穩(wěn)定性來(lái)代替其相位穩(wěn)定性，實(shí)現(xiàn)PS點(diǎn)的選擇。

任意PS點(diǎn)的大氣干涉相位Δφatm建模為［14］：

其中，λ為發(fā)射電磁波的波長(zhǎng)，Δn(r，t)為與時(shí)間t和空間r相關(guān)的大氣折射指數(shù)，L為電磁波傳輸路徑。在大氣空間同質(zhì)性良好的前提下，Δn在空間r不發(fā)生變化，因此公式（1）可改寫(xiě)為：

式中，β0、β1分別表示線性方程的常量與線性相關(guān)系數(shù)，R表示發(fā)射天線與PS點(diǎn)的斜距。

基于所有PS點(diǎn)建立線性方程組

ΔΦ為n個(gè)PS 點(diǎn)解纏相位構(gòu)成的n維向量，β為待估計(jì)線性系數(shù)所構(gòu)成的2 維向量，ε為隨機(jī)誤差向量，表示一次模型的誤差相位。利用最小二乘回歸對(duì)β0、β1進(jìn)行估計(jì)，得到

其中，T表示矩陣的轉(zhuǎn)置。因此，線性模型估計(jì)的大氣相位為：

由于部分PS 點(diǎn)中解纏相位存在噪聲或變形的影響，導(dǎo)致模型與觀測(cè)相位之間存在偏差。采用觀測(cè)相位和估計(jì)相位差的無(wú)偏估計(jì)來(lái)去除有偏點(diǎn)，即不可靠PS點(diǎn)［15］。

通過(guò)式（6）準(zhǔn)則剔除不可靠PS 點(diǎn)，再次迭代估計(jì)，直至不可靠PS點(diǎn)全部剔除，得到β0、β1準(zhǔn)確估計(jì)。

ΔTatm取值一般在0.1～0.2 rad內(nèi)。

除上述線性模型外，常用模型還包括方位角-斜距模型、二階斜距模型等，如式（7）、（8）。

其中，β0、β1、β2為待估計(jì)的模型參數(shù)，R、θ分別為PS點(diǎn)的斜距與方位角。當(dāng)大氣變化較為復(fù)雜時(shí)，上述模型均無(wú)法做到有效補(bǔ)償。

2.2 氣象數(shù)據(jù)校正法

由式（1）可知，大氣折射指數(shù)n僅與時(shí)間t和空間r相關(guān)。監(jiān)測(cè)區(qū)域較小或天氣良好時(shí)，大氣隨空間均勻變化，氣象參數(shù)可視為保持一致。t時(shí)刻回波中的大氣相位可以寫(xiě)為：

式中，φatm大氣擾動(dòng)相位只與距離、波長(zhǎng)、大氣折射指數(shù)有關(guān)。

根據(jù)大氣折射經(jīng)驗(yàn)?zāi)Ｐ?，可得t時(shí)刻大氣折射指數(shù)n(t)表達(dá)式為［16］：

式中，T為溫度，單位為K。P為大氣壓強(qiáng)，單位為mbar。h為相對(duì)濕度，單位為%。

2.3 改進(jìn)方法

2.3.1 噪點(diǎn)剔除

通過(guò)求解空間中相鄰PS 點(diǎn)差分干涉相位序列標(biāo)準(zhǔn)差，設(shè)定合理閾值，剔除受噪聲影響較大的PS點(diǎn)。

任意PS點(diǎn)差分干涉相位Δφ建模為：

其中，φdefo為形變相位分量，φatm為大氣相位分量，φnoi為噪聲相位分量?？紤]到三個(gè)分量彼此獨(dú)立，PS點(diǎn)的干涉相位序列標(biāo)準(zhǔn)差σInP可以表示為

其中，σdefo、σatm和σnoi分別對(duì)應(yīng)各相位分量標(biāo)準(zhǔn)差。假定空間中相鄰PS 點(diǎn)p1、p2，其差分相位Δφ(p1，p2)可表示為

在天氣良好條件下，大氣相位呈現(xiàn)大面積區(qū)域的高度空間相干性，即使復(fù)雜大氣條件下存在空變性區(qū)域，局部空間中仍然呈現(xiàn)高度相干性，即趨近于0。而噪聲相位在空間中則呈現(xiàn)無(wú)規(guī)則性，則臨近PS 點(diǎn)對(duì)p1，p2 差分相位序列標(biāo)準(zhǔn)差可以表示為［17］

基于Delaunay 算法對(duì)PS 點(diǎn)構(gòu)建三角形網(wǎng)絡(luò)。任意一個(gè)PS 點(diǎn)均能由一個(gè)或多個(gè)三角形找到臨近PS 點(diǎn)，定義為臨近PS 點(diǎn)對(duì)。去除三角形中較長(zhǎng)的邊，且保證任意PS 點(diǎn)至少存在一個(gè)臨近PS 點(diǎn)，通過(guò)計(jì)算臨近PS點(diǎn)對(duì)差分相位序列標(biāo)準(zhǔn)差，設(shè)定合理閾值，去除與其相鄰PS 點(diǎn)相位差異較大的點(diǎn)，即受噪聲影響較大或部分發(fā)生形變的PS點(diǎn)。

2.3.2 基于雙閾值的穩(wěn)定點(diǎn)選擇

由上述方法剔除受噪聲影響較大的PS點(diǎn)后，剩余PS點(diǎn)干涉相位可表示為：

其中，φnoi為噪聲相位分量，數(shù)值一般在0.1 rad 左右。若PS 點(diǎn)不發(fā)生形變，即為穩(wěn)定PS 點(diǎn)，則φdefo趨近于0，干涉相位Δφ與大氣相位φatm呈高度相關(guān)。通過(guò)計(jì)算PS 點(diǎn)大氣相位序列與差分干涉相位序列的MIC，設(shè)定合理閾值，實(shí)現(xiàn)部分穩(wěn)定PS點(diǎn)選擇。

MIC，即最大互信息系數(shù)，用于衡量?jī)蓚€(gè)變量X和Y之間的相關(guān)程度，包括線性、非線性的強(qiáng)度。相較于Pearson 相關(guān)系數(shù)只能衡量變量線性相關(guān)程度的局限，MIC能夠捕獲到變量之間各種線性、非線性關(guān)系，對(duì)噪聲不敏感，且具有更高的準(zhǔn)確度，符合復(fù)雜大氣條件下，各區(qū)域大氣折射率呈現(xiàn)各種線性或非線性關(guān)系的實(shí)際情況，能實(shí)現(xiàn)空間大部分區(qū)域穩(wěn)定PS點(diǎn)選擇。其表達(dá)式為［18］：

其中，X、Y為離散隨機(jī)變量，I[X；Y]為X、Y的互信息，表示為

式中，ρ(X)、ρ(Y)分別為X和Y的邊緣概率，ρ(X，Y)是X和Y的聯(lián)合概率。其計(jì)算方式通過(guò)將X、Y以散點(diǎn)圖的形式離散在二維空間中，將當(dāng)前空間沿X、Y方向上劃分為|X|，|Y|個(gè)格子，以網(wǎng)格中當(dāng)前散點(diǎn)個(gè)數(shù)的頻率作為邊緣概率與聯(lián)合概率。B是變量，通常為n0.6，其中n為數(shù)據(jù)量的個(gè)數(shù)。

式中，X、Y分別為各PS 點(diǎn)的大氣相位序列與差分干涉相位序列，通過(guò)計(jì)算MIC 并設(shè)定合理的閾值，實(shí)現(xiàn)部分穩(wěn)定PS點(diǎn)的選擇。由于大氣常常表現(xiàn)為復(fù)雜時(shí)空變性，部分區(qū)域部分時(shí)間段內(nèi)出現(xiàn)空變性，導(dǎo)致單一的MIC 閾值過(guò)低，出現(xiàn)部分區(qū)域缺少穩(wěn)定PS點(diǎn)，對(duì)空間大氣相位插值擬合過(guò)程中將會(huì)出現(xiàn)較大誤差。

為此，文章引入相干系數(shù)，評(píng)估PS 點(diǎn)的穩(wěn)定程度。PS 點(diǎn)主要有兩種失相干因素，即時(shí)間失相干（大氣改變）、空間失相干（位置改變）。部分區(qū)域在監(jiān)測(cè)過(guò)程中的部分時(shí)間段呈現(xiàn)空變性，為避免時(shí)間失相干，通過(guò)求解多組連續(xù)獲取的兩張雷達(dá)圖像的相干系數(shù)，設(shè)定較高閾值，實(shí)現(xiàn)部分空變區(qū)域穩(wěn)定PS點(diǎn)選擇。

相干系數(shù)法主要思路為在干涉相位對(duì)移動(dòng)窗口并計(jì)算窗口內(nèi)全體像元的相干信息作為當(dāng)前窗口中心像元的相干系數(shù)值。其數(shù)學(xué)表達(dá)式為［19］：

其中S1，S2分別為干涉相位對(duì)的兩個(gè)窗口，*指的是共軛相乘，m，n分別對(duì)應(yīng)窗口的行、列。

通過(guò)設(shè)定MIC 與相干系數(shù)的合理閾值，求取并集，實(shí)現(xiàn)空間各區(qū)域范圍內(nèi)穩(wěn)定PS 點(diǎn)選擇，避免存在部分空變區(qū)域缺失穩(wěn)定PS 點(diǎn)而引起較大的插值擬合誤差。

2.3.3 Kriging插值擬合

選出的穩(wěn)定PS 點(diǎn)在圖像各區(qū)域中呈現(xiàn)非均勻分布，通過(guò)空間維插值擬合，實(shí)現(xiàn)所有PS 點(diǎn)大氣相位估計(jì)。本文采用Kriging插值擬合。

Kriging 插值方法基于變異函數(shù)和結(jié)構(gòu)分析理論。本文根據(jù)待求PS 點(diǎn)一定鄰域范圍內(nèi)的穩(wěn)定PS點(diǎn)大氣相位，在考慮了其自身大小、形狀以及與待求樣點(diǎn)之間的空間位置信息，通過(guò)變異函數(shù)提供的結(jié)構(gòu)信息對(duì)待求PS 點(diǎn)進(jìn)行的一種最優(yōu)的線性無(wú)偏估計(jì)。其數(shù)學(xué)表達(dá)式為［20］：

其中，φatm(si)為第i個(gè)穩(wěn)定PS點(diǎn)處的大氣相位，N表示穩(wěn)定PS點(diǎn)個(gè)數(shù)。λi為對(duì)應(yīng)樣本點(diǎn)的未知權(quán)重，其滿足估計(jì)值與真實(shí)值φatm(s0)的差值最小，即

其中，γ(h)為半變異函數(shù)，表示待求PS 點(diǎn)與穩(wěn)定PS點(diǎn)距離的不同，影響程度不同。一般來(lái)說(shuō)，在對(duì)氣象要素場(chǎng)插值時(shí)半變異函數(shù)選用球狀模型，其兼顧了儲(chǔ)層參數(shù)的隨機(jī)性與相關(guān)性。數(shù)學(xué)表達(dá)式為：

其中，C為拱高，C0為塊金值，a為變程，表示空間中具有相關(guān)性的范圍。

通過(guò)穩(wěn)定PS 點(diǎn)空間位置以及差分干涉相位信息進(jìn)行Kriging 插值擬合，實(shí)現(xiàn)所有PS 點(diǎn)大氣相位的線性無(wú)偏最優(yōu)估計(jì)。

3 實(shí)驗(yàn)信息

本次實(shí)驗(yàn)區(qū)域?yàn)橹貞c市萬(wàn)州區(qū)九道拐，所監(jiān)測(cè)的山體邊坡縱向高程約170 m、橫向?qū)挾燃s300 m，且地處長(zhǎng)江北岸，時(shí)常受江風(fēng)襲擾，該區(qū)域大氣條件較為復(fù)雜［21］。

如圖1（a）所示，采用實(shí)驗(yàn)平臺(tái)為北京理工雷科公司研發(fā)的直線式掃描邊坡監(jiān)測(cè)雷達(dá)。該雷達(dá)工作波段為Ku，測(cè)量周期3～10 min，合成孔徑約1.8 m，1 km 處空間分辨率為0.3 m×4 m，測(cè)量精度較高，達(dá)±0.1 mm。圖1（b）、（c）所示分別為雷達(dá)監(jiān)測(cè)區(qū)域與氣象站。

圖1 實(shí)驗(yàn)設(shè)備與場(chǎng)地Fig.1 Experimental equipments and venue

本次采用數(shù)據(jù)為邊坡雷達(dá)連續(xù)獲取的30 張雷達(dá)圖像，時(shí)間從2020年12月18日10點(diǎn)40分～2020年12月18日15點(diǎn)13分，每幅圖像測(cè)量周期約為9分鐘。通過(guò)氣象站獲得該時(shí)間段內(nèi)的氣象數(shù)據(jù)（溫度、濕度、大氣壓），其變化曲線如圖2（a）、（b）、（c）所示。

圖2 氣象數(shù)據(jù)變化曲線圖Fig.2 Meteorological data change curve

圖3（a）所示為雷達(dá)圖像，采用最大的幅值作為參考值進(jìn)行了dB 處理。采用幅度離差法對(duì)30 幅圖像進(jìn)行PS 點(diǎn)選擇，圖3（b）為PS 點(diǎn)分布情況。通過(guò)設(shè)置幅度離差門(mén)限0.2，幅度門(mén)限?35 dB，共篩選出14842個(gè)PS點(diǎn)。

圖3 雷達(dá)圖像與PS點(diǎn)選擇結(jié)果Fig.3 Radar image and PSs selection result

分析雷達(dá)圖像時(shí)，采用連續(xù)獲取的兩張雷達(dá)圖像作為主、副影像進(jìn)行差分干涉，本次采用30 幅雷達(dá)圖像前后差分干涉，得到29張干涉相位圖。

短時(shí)間內(nèi)，PS點(diǎn)一般不發(fā)生形變或形變量較小，即φdefo在0 rad 左右。如圖4（a）、（b）所示，分別對(duì)應(yīng)連續(xù)兩幅雷達(dá)圖像的干涉相位圖與干涉相位散點(diǎn)圖。

圖4 相鄰時(shí)間A組干涉相位圖與散點(diǎn)圖Fig.4 Interferometric phase diagram and scatter plot of group A at adjacent times

可以看出，PS點(diǎn)的干涉相位主要分布在-0.5～0.5 rad 范圍內(nèi)，隨斜距發(fā)生線性變化。這說(shuō)明此段時(shí)間范圍內(nèi)，大氣變化較小，氣象條件較為穩(wěn)定，可以將大氣延遲帶來(lái)的相位誤差建模為隨斜距線性變化的分量，紅色實(shí)線代表的線性大氣相位估計(jì)可以起到很好的相位補(bǔ)償作用。

圖5（a）、（c）所示分別對(duì)應(yīng)B、C 兩組相鄰時(shí)間干涉相位圖?？梢钥闯?，干涉相位圖中呈明顯空變性且兩組干涉相位散點(diǎn)圖（b）、（d）中PS點(diǎn)的干涉相位較為分散，均隨斜距發(fā)生了非線性改變，且變化趨勢(shì)各不相同。說(shuō)明在兩組圖像的獲取時(shí)間范圍內(nèi)，大氣具有較強(qiáng)的時(shí)變與空變性，無(wú)法通過(guò)多項(xiàng)式模型估計(jì)大氣相位，常規(guī)方法不再適用。

圖5 相鄰時(shí)間B、C組干涉相位圖與散點(diǎn)圖Fig.5 Interference phase diagram and scatter diagram of two groups B and C at adjacent times

4 實(shí)驗(yàn)結(jié)果

4.1 噪點(diǎn)剔除

本次通過(guò)Delaunay三角網(wǎng)建立臨近PS點(diǎn)對(duì)，以所有三角形邊長(zhǎng)的平均值長(zhǎng)度作為參考，剔除超過(guò)該參考值（4.35 m）的長(zhǎng)邊，且保證任意PS點(diǎn)至少存在一個(gè)臨近PS點(diǎn)。

根據(jù)式（14）計(jì)算出的所有PS點(diǎn)的臨近標(biāo)準(zhǔn)差，若一個(gè)PS 點(diǎn)存在多個(gè)臨近PS 點(diǎn)，則取其平均。如圖6（a）所示，為所有PS 點(diǎn)臨近標(biāo)準(zhǔn)差圖。圖6（b）所示為臨近標(biāo)準(zhǔn)差散點(diǎn)圖，由于隨著斜距不斷增加，PS 點(diǎn)信噪比逐漸降低，臨近標(biāo)準(zhǔn)差也會(huì)相應(yīng)增大。因此臨近標(biāo)準(zhǔn)差采用隨距離呈線性相關(guān)的閾值門(mén)限。本次在最短觀測(cè)距離400 m采用0.35 rad，最遠(yuǎn)距離1100 m 處采用0.4 rad，如圖中紅色虛線所示。共篩選出597 個(gè)點(diǎn)，圖6（c）中黑色像素點(diǎn)即為PS點(diǎn)中篩選出的噪點(diǎn)或部分發(fā)生形變的點(diǎn)，無(wú)規(guī)則分布在各個(gè)區(qū)域中，符合噪聲特性。

圖6 PS點(diǎn)臨近標(biāo)準(zhǔn)差圖與噪點(diǎn)選擇結(jié)果Fig.6 PSs near standard deviation map and noise selection results

4.2 基于雙閾值的穩(wěn)定PS點(diǎn)選擇

經(jīng)過(guò)噪點(diǎn)剔除后，將剩余PS點(diǎn)差分干涉相位序列分別與對(duì)應(yīng)斜距下的大氣干涉相位序列進(jìn)行MIC計(jì)算。如圖7（a）所示為大氣折射指數(shù)差分曲線，由式（10）計(jì)算得到，根據(jù)式（9）得到對(duì)應(yīng)PS 點(diǎn)大氣干涉相位序列。圖7（b）所示為剩余PS 點(diǎn)MIC 圖，設(shè)定閾值為0.92，得到與大氣相位變化高度相關(guān)的穩(wěn)定PS 點(diǎn)8436 個(gè)，其分布圖如圖7（c）所示。圖中可以看出，穩(wěn)定PS 點(diǎn)非均勻分布在圖像中，但圖中紅色橢圓部分受大氣空變性影響無(wú)法選出穩(wěn)定PS點(diǎn)，在對(duì)其空間維插值補(bǔ)償時(shí)，將會(huì)產(chǎn)生較大誤差。

圖7 大氣折射指數(shù)差分曲線與高相關(guān)穩(wěn)定PS點(diǎn)選擇結(jié)果Fig.7 Atmospheric refractive index difference curve and high correlation stable PSs selection results

本文通過(guò)對(duì)29 張連續(xù)時(shí)間的干涉相位對(duì)的相干系數(shù)取平均，實(shí)現(xiàn)部分空變區(qū)域高相干穩(wěn)定PS點(diǎn)的選擇。圖8（a）所示，為圖像中所有像素點(diǎn)的相干系數(shù)圖。設(shè)定閾值為0.99，在經(jīng)噪點(diǎn)剔除后的剩余PS 點(diǎn)中共選取到高相干穩(wěn)定PS 點(diǎn)4288 個(gè)，如圖8（b）所示，圖中可以看出在紅色橢圓區(qū)域均選到穩(wěn)定PS 點(diǎn)，彌補(bǔ)MIC 單閾值選點(diǎn)的不足。穩(wěn)定PS點(diǎn)總數(shù)由雙閾值取并集得到，共計(jì)10185 個(gè)，占比68.62%。圖8（c）所示，即為本次選取的穩(wěn)定PS 點(diǎn)分布圖，穩(wěn)定點(diǎn)非均勻分布在圖像各個(gè)區(qū)域。

圖8 相干系數(shù)圖與穩(wěn)定PS點(diǎn)選擇結(jié)果Fig.8 Coherence coefficient map and stable PSs selection results

4.3 插值擬合結(jié)果

以本文所提方法篩選出穩(wěn)定PS點(diǎn)后，利用其位置信息（X、Y坐標(biāo)值）以及29 張差分干涉相位序列作為已知樣本點(diǎn)，半變異函數(shù)選用球狀模型，參考式（22），進(jìn)行Kriging 插值擬合實(shí)現(xiàn)所有PS 點(diǎn)大氣相位的線性無(wú)偏最優(yōu)估計(jì)，從而進(jìn)行相位補(bǔ)償。

如圖9（a）所示為本次對(duì)29 張干涉相位圖補(bǔ)償后的干涉相位累積圖，相較于圖9（b）常規(guī)線性補(bǔ)償與圖9（c）氣象數(shù)據(jù)校正方法下的干涉相位累積圖，所有PS 點(diǎn)補(bǔ)償后的相位均在0 rad 左右，大氣的空變性得到了有效改善。而兩種常規(guī)方法補(bǔ)償效果較差，噪點(diǎn)較多，且圖中左側(cè)PS 點(diǎn)干涉相位在-0.5 rad 左右，右側(cè)點(diǎn)在0.5～1 rad 之間，此為29張干涉相位序列中大氣相位補(bǔ)償精度較低，殘余大氣相位不斷疊加而導(dǎo)致。

圖9 不同方法下的干涉圖補(bǔ)償結(jié)果Fig.9 Interferogram compensation results under different methods

為了更加直觀地看出殘余大氣相位的疊加影響，在不同斜距、方位角處選擇出幅度離差最小的PS點(diǎn)作為參考點(diǎn)，分析隨時(shí)間序列的干涉相位變化圖。本次共選取6 個(gè)參考點(diǎn)，其空間分布圖如圖10（a）所示。圖10（b）為常規(guī)線性補(bǔ)償后參考點(diǎn)的相位變化曲線，可以看出參考點(diǎn)相位明顯偏離0 rad，且變化趨勢(shì)各不相同。這說(shuō)明常規(guī)線性補(bǔ)償誤差較大、精度較低，存在嚴(yán)重的過(guò)、欠補(bǔ)償現(xiàn)象，由大氣帶來(lái)的相位誤差隨時(shí)間和空間劇烈變化，波動(dòng)幅度在-1.5～2 rad 之間，約2.9 mm 形變誤差。圖10（c）為氣象數(shù)據(jù)校正后參考點(diǎn)相位變化曲線，其波動(dòng)較常規(guī)線性補(bǔ)償更為劇烈，補(bǔ)償效果最差，波動(dòng)幅度在-1～3 rad 之間，約4.3 mm 形變誤差。如圖10（d）所示為本次改進(jìn)方法下的參考點(diǎn)相位變化曲線，圖中看出參考點(diǎn)相位波動(dòng)均在0 rad 左右，且相位序列變化隨機(jī)，不受大氣相位疊加影響，補(bǔ)償效果最優(yōu)。其中波動(dòng)原因主要來(lái)自于噪聲，噪聲相位一般在±0.1 rad之間，影響較小。

圖10 不同方法下的參考點(diǎn)相位累積變化圖Fig.10 The cumulative change diagram of the reference point phase under different methods

5 結(jié)論

GB-InSAR 形變監(jiān)測(cè)中，大氣變化所帶來(lái)的相位延遲是主要的誤差源之一。在復(fù)雜大氣條件下，常規(guī)相位補(bǔ)償方法不再適用，文章所提的一種基于PS技術(shù)的改進(jìn)方法，解決了常規(guī)方法無(wú)法對(duì)大氣波動(dòng)劇烈時(shí)間段內(nèi)的干涉相位圖進(jìn)行有效補(bǔ)償?shù)膯?wèn)題。通過(guò)對(duì)九道拐實(shí)測(cè)數(shù)據(jù)分析，比較了改進(jìn)方法與常規(guī)方法下的參考點(diǎn)的相位變化曲線和累積相位圖，驗(yàn)證了本文方法在復(fù)雜大氣條件下相位補(bǔ)償?shù)挠行浴?/p>

本文方法仍存在一些欠缺之處。首先在穩(wěn)定PS點(diǎn)選擇過(guò)程中，需要人為憑經(jīng)驗(yàn)設(shè)定MIC 以及相干系數(shù)閾值；其次當(dāng)雷達(dá)圖像較少時(shí)，MIC無(wú)法有效捕捉到PS 點(diǎn)的大氣相位序列與差分干涉相位序列之間的線性、非線性關(guān)系，導(dǎo)致穩(wěn)定PS點(diǎn)的漏選；最后在插值擬合過(guò)程中，Kriging 插值算法復(fù)雜度較高，插值擬合效率較慢，對(duì)于應(yīng)用于時(shí)序GB-InSAR圖像實(shí)時(shí)處理中還需要進(jìn)一步研究。