InSAR形變場最佳插值算法對比研究

李 巖，韓 蕾

(1.中國地震局第二監(jiān)測中心，陜西西安 710054；2.西安歐亞學(xué)院人居環(huán)境學(xué)院，陜西西安 710065)

1 引 言

合成孔徑雷達(dá)干涉測量(InSAR)因具備全天時(shí)、全天候、大范圍監(jiān)測，有一定穿透能力等優(yōu)勢，已廣泛應(yīng)用于地震同震形變等地殼運(yùn)動領(lǐng)域[1-2]。通過對形變前后SAR影像進(jìn)行干涉處理，提取地形和地表形變信息，可以實(shí)現(xiàn)地表運(yùn)動的有效監(jiān)測[3-4]。

在應(yīng)用過程中，由于衛(wèi)星原始影像畸變、處理參數(shù)設(shè)置不當(dāng)、地形地貌復(fù)雜等因素造成干涉圖像相干性差異很大，導(dǎo)致形變信息不全，影響后續(xù)分析及反演結(jié)果的精度。

為了解決上述問題，多種圖像插值算法被引進(jìn)到InSAR處理過程以提高形變場質(zhì)量，但不同的插值算法會導(dǎo)致獲取的最終形變信息不一致。針對這一問題，本文選取了4種常用的圖像插值算法，并通過真實(shí)震例對比分析了各算法優(yōu)缺點(diǎn)，為準(zhǔn)確獲取插值形變場提供技術(shù)參考。

2 常用插值算法

InSAR形變場的精度和保真度是進(jìn)行地震同震破裂和應(yīng)力變化反演的重要保障。由于InSAR形變場像元的物理屬性，以及圖像插值算法理論模型的不同，使得最終插值結(jié)果與實(shí)際形變存在差異。當(dāng)前常用的插值算法包括最鄰近插值、雙線性插值、立方卷積插值和正弦函數(shù)插值[5-7]。

2.1 最鄰近插值法

最鄰近插值法又稱零階插值，該方法通過反向變換，選取計(jì)算像元周圍4個(gè)原始像元中距離最近的一個(gè)原始像元的灰度值作為新值點(diǎn)。此算法簡單且直觀，但是當(dāng)圖像中包含像元之間灰度級的細(xì)微結(jié)構(gòu)變化時(shí)，應(yīng)用最鄰近插值法難以保證結(jié)果的精度。

2.2 雙線性插值法

雙線性插值法又稱一階插值，其核心思想是在兩個(gè)方向分別進(jìn)行一次線性插值，計(jì)算像元的灰度值，根據(jù)兩個(gè)方向上各兩個(gè)最鄰近像元的灰度值加權(quán)平均計(jì)算得到。因此該算法比最鄰近插值法計(jì)算量大，精度明顯提高，不會出現(xiàn)像元值不連續(xù)的情況。但是雙線性插值法具有低通濾波器的性質(zhì)，會使得高頻分量受損，導(dǎo)致對比度明顯的分界線在一定程度上變得模糊。

2.3 立方卷積插值法

立方卷積插值法能夠克服上述兩種算法的不足，其主要思想是取計(jì)算像元點(diǎn)周圍的16個(gè)原始像元點(diǎn)作為參考值，在水平方向、垂直方向上分別進(jìn)行三階插值，最后計(jì)算得到計(jì)算像元的灰度值。由于考慮了計(jì)算像元點(diǎn)周圍更多原始像元點(diǎn)的相關(guān)性，因此計(jì)算得到的像元灰度值也會更接近真實(shí)值，可以比雙線性插值法更加清楚地表現(xiàn)細(xì)節(jié)。但由于該算法考慮了16個(gè)像元點(diǎn)的值，且插值函數(shù)是三階函數(shù)，因此計(jì)算量也明顯增大。

2.4 正弦函數(shù)插值法

正弦函數(shù)插值法是對原始像元進(jìn)行函數(shù)運(yùn)算后用曲線將各個(gè)像元點(diǎn)連接起來，通過這個(gè)重建的函數(shù)求出任意位置處像元的灰度值。較其他插值方法，正弦函數(shù)法能有效減少插值所造成的信號損失，較好地保留形變場的有效信息。但是為了精確計(jì)算某一像元值，正弦函數(shù)插值需要覆蓋無限多個(gè)點(diǎn)，這在實(shí)際應(yīng)用中是無法實(shí)現(xiàn)的，并且大量的數(shù)據(jù)點(diǎn)導(dǎo)致插值過程十分耗時(shí)，而精度提高效果并不顯著。

3 算例及結(jié)果分析

3.1 數(shù)據(jù)選擇和InSAR處理

為了對上述4種插值算法進(jìn)行統(tǒng)一的分析、比較及驗(yàn)證，本文選用2幅覆蓋2003年12月26日伊朗巴姆地震的C波段 ENVISAT ASAR降軌數(shù)據(jù)和2幅覆蓋2010年4月14日青海玉樹地震的L波段ALOS PALSAR升軌數(shù)據(jù)(見表1)，采用高級雷達(dá)圖像處理軟件ENVI/SARscape(試用版)[8]進(jìn)行兩軌法差分干涉測量[9]，濾波方法選用較為健壯的Goldstein法[10]，相位解纏選用最小費(fèi)用流法[11]。兩例試驗(yàn)數(shù)據(jù)覆蓋性好，相干性高，得到的干涉圖連續(xù)完整(圖2a、圖3a)，試驗(yàn)流程如圖1所示。通過設(shè)定相干閾值為0.5，即只對0.5以上的差分干涉圖的相位值進(jìn)行解纏，得到包含空值的試驗(yàn)形變場樣本(圖2b、圖3b)，分別進(jìn)行上述4種插值算法的內(nèi)插與地理編碼。

圖1 SAR 影像處理流程圖Fig.1 Flow chart of SAR image processing

圖2 2003年12月26巴姆地震InSAR結(jié)果及不同插值結(jié)果Fig.2 InSAR results and different interpolation results of the Bam earthquake on 26 December 2003

圖3 2010年4月14日玉樹地震InSAR結(jié)果及不同插值結(jié)果Fig.3 InSAR results and different interpolation results of the Yushu earthquake on 14 April 2010

3.2 插值結(jié)果分析對比

偏斜度和陡峭度分布是對統(tǒng)計(jì)數(shù)據(jù)分布形態(tài)對稱性和陡緩程度的度量。從試驗(yàn)結(jié)果可以看出，地震造成的地物地貌破壞，在InSAR干涉圖表現(xiàn)為低相干性，導(dǎo)致初始形變樣本和實(shí)驗(yàn)樣本之間數(shù)據(jù)分布形態(tài)明顯差異，且形變場中空值越大，差異越大。表1為試驗(yàn)用SAR數(shù)據(jù)參數(shù)，表2、表3為不同插值算法的形變結(jié)果對比。

表1 試驗(yàn)用SAR數(shù)據(jù)參數(shù)Tab.1 SARdataparametersusedintests地震衛(wèi)星波長/cm可探測最小形變量/m影像時(shí)間空間基線/m時(shí)間基線/d2003-12-26巴姆ENVISATASAR5.6(C波段)0.028VV2003-12-32004-02-11-3.069702010-04-14玉樹ALOSPALSAR23.6(L波段)0.118HH2010-01-152010-04-17699.57492

表2 巴姆地震不同插值算法的形變結(jié)果對比Tab.2 ComparisonofdeformationresultsoftheBamearthquakewithdifferentinterpolationalgorithms數(shù)據(jù)量/個(gè)最大值/m最小值/m平均值/m均方根標(biāo)準(zhǔn)差/10-5偏斜度陡峭度初始形變樣本8888090.3085-0.16430.0036360.06744297.143280.974765.95775實(shí)驗(yàn)樣本(相干性≥0.5)8510550.3085-0.16420.0050930.06772647.320621.001715.96939最鄰近插值法8869450.3086-0.16420.0035050.06754637.162550.977975.94411雙線性插值法8867150.3084-0.16410.0035160.06754667.163460.978025.94410三階立方卷積插值法8867150.3085-0.16420.0035160.06754797.163590.978005.94406四階立方卷積插值法8867150.3085-0.16420.0035160.06754807.163600.977995.94404正弦函數(shù)插值法8868260.3088-0.17050.0035050.06742957.150620.977935.94493

表3 玉樹地震不同插值算法的形變結(jié)果對比Tab.3 ComparisonofdeformationresultsoftheYushuearthquakewithdifferentinterpolationalgorithms數(shù)據(jù)量/個(gè)最大值/m最小值/m平均值/m均方根標(biāo)準(zhǔn)差/10-4偏斜度陡峭度初始形變樣本22392000.3315-0.4290.0621240.11427201.11013-2.007417.39076實(shí)驗(yàn)樣本(相干性≥0.5)20050960.3245-0.49410.0460170.12460421.34053-1.608514.85072最鄰近插值法22262580.3244-0.48080.0560890.12345741.27301-2.147097.32938雙線性插值法22250180.3208-0.47640.0560060.12348551.27386-2.145287.31159三階立方卷積插值法22250370.3257-0.47900.0560220.12351331.27413-2.144657.31145四階立方卷積插值法22250360.3260-0.47960.0560220.12351611.27416-2.144577.31139正弦函數(shù)插值法22254340.3090-0.47920.0559680.12322561.27071-2.148297.33467

從內(nèi)插結(jié)果來看，4種插值方法都使得干涉形變場連續(xù)清晰。但不同的插值方法對形變邊緣梯度處理存在差異。由表2和表3所列的不同插值算法的形變結(jié)果對比可以看出，最鄰近插值法的數(shù)據(jù)量最多，由于該方法只考慮最近的像元，因此最大值和最小值接近試驗(yàn)樣本，平均值遠(yuǎn)小于試驗(yàn)樣本，也小于初始樣本，偏斜度與陡峭度也較其他方法與試驗(yàn)樣本相差較大，在插值形變場邊緣出現(xiàn)了明顯的鋸齒現(xiàn)象。雙線性插值法是最鄰近插值法的一種改進(jìn)，基本克服了插值結(jié)果跳躍的缺點(diǎn)，其插值結(jié)果與離散程度較穩(wěn)定，但是數(shù)據(jù)量最少，邊緣和細(xì)節(jié)插值被抑制，出現(xiàn)輪廓模糊現(xiàn)象。三階立方卷積插值法與四階立方卷積插值法的結(jié)果較一致，最大值、最小值和平均值最接近初始樣本，由于四階收斂的核函數(shù)比三階收斂的核函數(shù)更加接近于階躍函數(shù)，因此四階立方卷積插值法的結(jié)果較平滑，精度更加可靠。正弦函數(shù)插值法的數(shù)據(jù)量僅次于最鄰近插值法，均方根和數(shù)據(jù)分布形態(tài)都最接近初始形變樣本，但是最大值、最小值和平均值都與初始樣本和實(shí)驗(yàn)樣本差異較大。

4 結(jié)束語

InSAR干涉圖的形變信息不全會對后續(xù)的分析反演結(jié)果產(chǎn)生不良影響，隨著形變空值越多，反演破裂結(jié)果和真實(shí)地震破裂的差異越明顯。為了提高形變場質(zhì)量，將多種圖像插值算法引進(jìn)到InSAR處理過程中。本文利用4種常用的圖像插值算法，通過真實(shí)震例進(jìn)行試驗(yàn)對比分析，得到以下結(jié)論：

1)最鄰近插值法只考慮最近的像元值，計(jì)算簡單，但插值后的形變場邊緣容易產(chǎn)生鋸齒效應(yīng)，造成形變值的不平滑?；谄渌惴ㄔ恚钹徑逯捣ㄌ貏e適用于數(shù)據(jù)分類的抽樣。

2)雙線性插值法、三階立方卷積插值法和四階立方卷積插值法分別考慮周圍4、8和16個(gè)像元。一般來說，高階插值可以提供更好的結(jié)果。本文試驗(yàn)結(jié)果表明，雙線性插值法插值的形變場邊緣受到平滑作用，出現(xiàn)輪廓模糊現(xiàn)象，四階立方卷積插值法與三階立方卷積插值法的結(jié)果幾乎無差異。

3)正弦函數(shù)插值法考慮其周圍的256個(gè)像元，因此運(yùn)算量會大大增加，數(shù)據(jù)數(shù)量、精度和分布形態(tài)都最接近初始形變樣本，但是數(shù)據(jù)值與樣本值差異較大。