顧及協(xié)方差函數的自適應四叉樹InSAR數據壓縮算法

張 靜,張 勤,2,趙超英,2,張菊清

1.長安大學地質工程與測繪學院,陜西西安 710054;2.西部礦產資源與地質工程教育部重點實驗室,陜西西安 710054

顧及協(xié)方差函數的自適應四叉樹InSAR數據壓縮算法

張 靜1,張 勤1,2,趙超英1,2,張菊清1

1.長安大學地質工程與測繪學院,陜西西安 710054;2.西部礦產資源與地質工程教育部重點實驗室,陜西西安 710054

利用InSAR變形監(jiān)測結果進行形變機理反演時,由于InSAR獲取的數據點眾多,且往往含有較多的誤差乃至粗差點,嚴重制約了反演計算的效率和可靠性。為此,本文提出顧及InSAR變形監(jiān)測數據的物理空間相關性來設立協(xié)方差函數,并依據協(xié)方差函數確定四叉樹象限分解閾值和最大象限大小的自適應四叉樹分解InSAR數據壓縮算法。本算法能夠在盡可能保留形變信號特征細節(jié)信息的同時,極大地降低InSAR數據量。論文以西安地區(qū)地面沉降InSAR形變監(jiān)測結果為例進行了試驗分析,驗證了該算法的有效性。結果表明,該方法能夠在不損失形變信號特征的同時,有效地實現(xiàn)InSAR數據壓縮和噪聲消除的目的。

四叉樹;數據壓縮;空間相關;協(xié)方差函數

1 引 言

合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)是近幾十年迅速發(fā)展起來的一種能夠精確測定地表三維信息和高程變化信息的微波遙感新技術。它具有高分辨率、高覆蓋和不受區(qū)域環(huán)境與氣候約束的特征,能較好地解決傳統(tǒng)地表形變觀測中三維空間上由于人工布設觀測點獲取的形變信息量不足和惡劣區(qū)域信息采集困難等問題[1-2],故可獲取更多的形變信號和揭示更多的地球物理現(xiàn)象,為地球物理反演研究提供一種全新的途徑,因此InSAR被眾多學者認為是一種極具潛力的空間對地觀測新技術[3]。然而,由于InSAR獲取的是基于像素的連續(xù)面狀形變數據,其包含的數據量較大,不但包含了大量冗余數據,而且還包含不少強噪聲和偽信號等,這些冗余數據和粗差噪聲不僅影響數據后處理的計算效率,而且還會因誤差的存在對InSAR結果解譯和形變機理反演等的可靠性帶來極為不利的影響[4]。因此,在保持形變信息基本特征不變的前提下,在有效地壓縮InSAR數據量的同時能消除一定的噪聲影響,是InSAR數據后處理的重要任務之一[5]。

無論何種影像的數據壓縮都應遵循盡可能保留原數據的基本特征信號的原則。常用的影像壓縮算法是通過將影像變換到頻率域或小波域[6-9],通過僅存儲小部分的頻域或小波系數來達到影像壓縮的目的,但是這類方法是對數據存儲大小的壓縮而非數據量大小本身的壓縮。而基于InSAR形變結果的反演研究則需要對數據量本身進行壓縮,因此,這類方法并不適用。針對此,文獻[10]提出了InSAR數據點的規(guī)則采樣法,該方法計算簡單,但需要在數據壓縮率和形變細節(jié)間進行取舍。文獻[4,11]分別利用象限分解的方法實現(xiàn)InSAR數據的重采樣以達到數據壓縮目的。文獻[12]則提出了基于像素分辨率的InSAR數據減采樣方法。這些方法均能在一定程度上達到InSAR數據量的有效壓縮,但其更多的是針對地震、火山等大尺度形變,需設定某些形變震源參數等,不具有普適性?；诖?本文將四叉樹分解法引入InSAR監(jiān)測數據的壓縮中,四叉樹分解是一種基于均勻性檢測的圖像分割算法[13-14];同時,考慮到InSAR數據自身誤差的特點,本文提出了顧及InSAR監(jiān)測數據的物理空間相關性設定協(xié)方差函數對四叉樹分解中的參數進行設定的方法,提高了算法的可靠性。這種顧及協(xié)方差函數的自適應四叉樹壓縮算法在較好保持形變信息的同時,可對數據量進行有效的壓縮,且能在一定程度上達到噪聲消除的目的,具有較好的普適性。

本文提出的InSAR數據壓縮算法,主要針對InSAR相位監(jiān)測數據進行分析,該數據既可是解纏相位,也可是經過地理編碼等轉換后的形變數據,算法具有通用性。本文主要以InSAR解纏相位數據為例進行算法的解釋和分析。

2 四叉樹InSAR數據壓縮

在形變模型反演中,為獲取最優(yōu)的模型參數,常常需要采用非線性迭代法進行大量的模型預測估計。而基于像素的連續(xù)面狀覆蓋的InSAR變形監(jiān)測數據,則含有成千上萬數以百萬計的數據點,其中包含了大量的冗余數據和強噪聲等引起的誤差乃至粗差等,利用這些形變數據點反演需要進行成千上萬次的模型迭代估計,這無疑給模型的解算帶來了不可估量的計算困難和負擔,而且異常噪聲的存在,將導致反演結果的嚴重失真[12]。因此有必要在盡可能地保持形變重要細節(jié)信息不變的條件下,尋求合適的InSAR變形監(jiān)測數據壓縮算法,這是InSAR數據后處理和形變模型反演研究的一個重要內容。本文研究將自適應四叉樹分解方法引入InSAR數據壓縮,該算法可在保持形變基本特征信息的同時,對數據量進行有效的壓縮。

2.1 四叉樹分解原理

四叉樹分解是一種基于均勻性檢測的圖像分割算法。自適應四叉樹分解,則是為了以最少的數據量,得到各個像塊的均值,從而獲得對原始圖像的最佳逼近[13]。對影像進行四叉樹分解的基本流程為:首先通過添加零值將影像的維數擴充為2的冪次,四叉樹分解要求影像的大小為N＝2n,其中,n為正整數,N為影像的維數。然后將影像分成大小相等的4個象限,分別計算各個象限內解纏值的方差,并將其與預先設定的方差閾值進行比較,若某象限方差超過閾值,則將該象限進一步分解為4個新的象限;重復上述過程直到收斂,即各個象限方差均小于給定的閾值;最后將分解后的各個象限的均值或中值賦給象限的中心坐標,僅保留該象限中心坐標點即可[15]。

2.2 算法的優(yōu)缺點分析

四叉樹壓縮算法的主要優(yōu)點是可在形變梯度變化較大處保留較多數據點,而在數據平坦即連續(xù)處保留較少點,因此相較于常規(guī)的規(guī)則降采樣等方法,它能更好地在保留形變細節(jié)信息不損失的同時,還可實現(xiàn)有效的數據壓縮;同時由于采用均值算法為象限賦值,其本質上也是一種低通濾波,可有效地消除影像中噪聲占主導地位的高頻部分,實現(xiàn)噪聲消除。但是,該算法也有不足之處:若某個分解象限內的形變值整體具有較好的一致性而僅在某個較小的區(qū)域內存在形變特征,盡管該小形變區(qū)域內的方差可能會明顯高于給定的方差閾值,但由于該象限整體的方差值比分解閾值低很多,在最后運算時僅將一個值賦給該象限的中心坐標,從而忽略該小形變區(qū)域,易造成形變細節(jié)信息的損失。為解決該問題,本文在進行四叉樹分解時,通過設定最大象限大小參數對其進行改進。

結合以上可知,在進行四叉樹分解時,需確定兩個重要參數:象限分解閾值和最大象限大小。這兩個參數的確定,需要對數據噪聲的統(tǒng)計特性進行分析,為此,本文引入噪聲協(xié)方差函數進行計算。

3 InSAR噪聲協(xié)方差函數估計

由于大氣和電離層等影響通常在數十到數百千米的尺度范圍內是相關的,且和高程相關的大氣水汽變化也會引起和地形相關的衛(wèi)星視向信息延遲等[16],因此InSAR獲取的干涉信號中存在空間相關信息[17],在測量中,這些信號被認為是觀測噪聲,而且這些噪聲的空間相關度同樣會影響InSAR監(jiān)測數據的質量。為有效地利用InSAR觀測結果進行形變機理的反演和比較等,必須對這些噪聲信號進行估計和去除。特別在InSAR監(jiān)測數據中的這種相關性,更主要地表現(xiàn)為一種物理空間相關,不能直接用函數擬合等簡單的數學模型構建,對于此,本文給出了利用協(xié)方差函數對這些噪聲信號進行估計和分析的方法。

在進行噪聲協(xié)方差函數Cn估計時,假設噪聲滿足二階平穩(wěn)且各向同性,即兩點間的協(xié)方差僅與其距離相關而和點的位置及方向無關。設r為兩點間距離,則其協(xié)方差可表示為

協(xié)方差函數的計算通常有兩種方法[12],一種是參照變異函數的計算,利用一個簡單的經驗協(xié)方差函數式求解

式中,P(r)為距離r的點對數目;f(xi)和f(xi＋ri)為影像在xi和xi＋ri處的取值;為影像的局部均值。

另一種則是根據二維自相關計算。在利用二維自相關計算協(xié)方差時,首先需對數據中的無效值進行0值填充,然后利用功率譜的逆傅里葉變換計算二維自相關,最后利用二維自相關的旋轉平均值計算一維協(xié)方差。在頻域中會在影像的邊界處引入常數變化,但是由于影像的均值為0,其基本上是正負相抵的,故可忽略。唯一需要的是利用數據中的非0點數對功率譜進行加權以使自相關正則化。影像gij的功率譜可通過離散二維傅里葉變換Gkl計算

式中,Gkl為gij的二維離散傅里葉變換;G為Gkl的共軛;Nnz為gij中的有效值個數。

求得自相關后,即可利用二維自相關函數的旋轉平均值估計一維協(xié)方差函數。如圖1所示即為利用二維自相關計算協(xié)方差函數的流程圖。

圖1 利用二維自相關計算協(xié)方差函數流程Fig.1 Covariance function calculation based on 2D autocorrelation

協(xié)方差函數給出的是一個兩點間相似性的測度,當兩點間的距離大于某個特定值時,協(xié)方差值趨于0,表明當兩點間的距離大于該特定值時,數據點間不相關。

4 算法的實現(xiàn)

在求取噪聲協(xié)方差函數時,僅關心噪聲的統(tǒng)計特性,故需要分離進行形變信息和噪聲的估計,可通過分析干涉條紋確定可能的形變區(qū)域,再利用掩膜處理將該區(qū)域形變從干涉圖中去除;另外,干涉圖中通常還含有因干涉基線估計不正確所引起的軌道殘余誤差[18-19],由于軌道誤差通常呈現(xiàn)為明顯的線性趨勢分布,故可利用線性擬合估計的方法提取一個最佳擬合平面將其剔除[20]。由此,可認為剩余的信息即為局部均值接近于0的噪聲數據,則可利用文中上節(jié)提到的方法進行噪聲協(xié)方差函數的估計。

由估計的噪聲協(xié)方差函數即可確定四叉樹分解的兩個重要參數:象限分解閾值和最大象限大小。象限分解閾值原則上應遠大于噪聲方差,但足夠小以保留形變細節(jié)特征不損失,通?？扇?倍噪聲方差;最大象限大小可設為影像中使數據點間不相關的距離。這兩個參數都可以根據估計的噪聲協(xié)方差函數確定。利用確定的參數,采用四叉樹對InSAR數據進行分解:即當某個分解象限的形變方差小于給定閾值時,將該象限內影像的均值賦給該象限中心坐標,在數據輸出時僅保留該象限中心值。圖2即為本文算法的完整流程圖。

圖2 顧及協(xié)方差函數的InSAR數據四叉樹壓縮算法流程圖Fig.2 Flow chart of the quadtree InSAR data reduction method based on covariance function

5 應用實例

本文以西安地區(qū)的InSAR解纏相位數據為例進行了試驗分析,圖3即為研究區(qū)解纏干涉圖,其中形變發(fā)生的主要區(qū)域為西安的高新開發(fā)區(qū),其沉降量級約為10 cm左右。圖4為利用二維自相關估計的噪聲協(xié)方差函數,從該協(xié)方差函數可看到,最大噪聲方差約為0.18 cm2,取四叉樹分解閾值為4倍噪聲方差,即0.72 cm2。進一步從圖上可看到,當點間距離大于0.8 km時,協(xié)方差函數接近于0,認為數據點間不相關,故可取象限的最大距離為0.8 km,即最大象限的面積為(0.64 km2)?？紤]到影像的像素大小為20 m,相當于40個像素大小,故取最大象限為64(需為2N)。圖5為利用四叉樹壓縮算法提取的數據散點圖,其數據量從原始的560 000個點降為44 414個點,壓縮后的數據量僅約為原數據的7.93%。從圖上可看到在形變發(fā)生區(qū)數據點采樣密集,充分保留了形變的細節(jié)特征信息量,而在數據變化平緩或無形變處僅保留了少量的數據點,從而大大地降低了數據的量級;且利用四叉樹壓縮數據能夠很好地對原數據進行重構,如圖6所示,利用四叉樹壓縮數據得到的解纏圖很好地保留了原解纏圖的形變細節(jié)信息,且相較于原圖其連續(xù)性更好,在部分噪聲明顯區(qū)域進行了有效的噪聲消除。為進一步比較四叉樹壓縮算法的效果,對圖3和圖6中小框內的區(qū)域分別進行放大顯示,結果如圖7和圖8所示。

從圖7和圖8的研究區(qū)域放大視圖可看出,其很好地保留了形變信號的細節(jié)信息,幾乎無細節(jié)損失,但數據點數卻降低了2～3個數量級,可有效地提高后續(xù)模型優(yōu)化反演等操作的運算效率。另外,通過對圖7和圖8的觀察分析可看到,利用顧及協(xié)方差的四叉樹壓縮算法得到的解纏圖噪聲明顯變小,相位變化連續(xù)性更好。為對比說明本文方法的效果,進一步提取了圖3和圖6中的兩條剖線A和B進行分析,結果如圖9所示,從剖線圖上可更直觀地看到利用該算法在一定程度上對數據中的部分噪聲進行了有效消弱,數據整體更光滑,達到了噪聲消除的目的。

圖3 西安解纏影像Fig.3 Unwrapping image of Xi’an

圖4 協(xié)方差函數Fig.4 Covariance function

圖5 四叉樹分解數據點圖Fig.5 Scatter points of Quadtree decomposition

圖6 四叉樹壓縮數據重構影像Fig.6 The image by quadtree reduction data

圖7 裁剪區(qū)原影像Fig.7 Sub-image of the full image

圖8 裁剪區(qū)四叉樹分解影像Fig.8 Quadtree decomposition image of the sub-image

圖9 研究區(qū)剖線分析(左為剖線A,右為剖線B)Fig.9 Profile analysis of the area

6 結 論

本文提出了顧及InSAR數據物理空間相關特性設立協(xié)方差函數的自適應四叉樹分解InSAR數據壓縮算法,并以西安地區(qū)的InSAR解纏相位數據為例進行了試驗分析,充分地證明了該算法的有效性。從結果可看到,該算法能夠在形變變化明顯處進行密集采樣,而在形變變化平坦處進行稀疏采樣,從而能夠在較好地保留In-SAR數據的細節(jié)信息不損失的條件下,達到有效的InSAR數據量壓縮目的,使InSAR數據點的數量降低2～3個量級,能有效地提高后續(xù)模型優(yōu)化反演等操作的運算效率,大大節(jié)省計算時間和硬件損耗;同時,由于采用了均值運算作為四叉樹分解的算子,在保證了該算法為線性的條件下,可有效地消除數據中的部分高頻噪聲,減小相位的突變,使數據整體連續(xù)性增強,數據更加平滑。

[1] FERRETTI A,MONTI G A,PRATI C,et al.InSAR Principles:Guidelines for SAR Interferometry Processing and Interpretation[M].AG Noordwijk:ESA Publications,2007.

[2] LUO Haibin,HE Xiufeng,LIU Yanxiong.Estimation of Three-dimensional Surface Motion Velocities Using Integration of DInSAR and GPS[J].Acta Geodaetica et Cartographica Sinica,2008,37(2):168-171.(羅海濱,何秀鳳,劉焱雄.利用DInSAR和GPS綜合方法估計地表3維形變速率[J].測繪學報,2008,37(2):168-171.)

[3] LIU Guoxiang.Surface Deformation Monitoring with Radar Interferometry[M].Beijing:Surveying and Mapping Press,2006.(劉國祥.利用雷達干涉技術監(jiān)測區(qū)域地表形變[M].北京:測繪出版社,2006.)

[4] JONSSON S.Modeling Volcano and Earthquake Deformation from Satellite Radar Interferometric Observations, [D].California:Stanford University,2002.

[5] ZH ANG Jing,ZH AO Chaoying,ZH ANG Qin,et al.Research on Deformation Contour Mapping from Interferometry[J].Journal of Geodesy and Geodynamcs,2009,29 (5):143-146.(張靜,趙超英,張勤,等.InSAR形變等值線圖繪制方法研究[J].大地測量與地球動力學,2009,29 (5):143-146.)

[6] DU Xutao.Discussed on the Image Compression Algorithm Based on Transformation[J].Software Guide,2012,11(1):145-148.(獨旭濤.基于變換的圖像壓縮算法探討[J].軟件導刊,2012,11(1):145-148.)

[7] LI Sha,HUANG Lin,ZHOU Jianqi,et al.SAR Image Compression Algorithm Based on Wavelet Transform[J].Computer Measurement and Control,2012,20(8):2310-2313.(李莎,黃琳,周劍奇,等.基于小波變換的SAR圖像壓縮算法研究[J].計算機測量與控制,2012,20(8): 2310-2313.)

[8] XIE Haihui.Speckle Reduction and Data Compression in SAR Images Using Wavelet-Based Method[D].Chengdu:University of Electronic Science and Technology,2004.(謝?；?干斑抑制及圖像壓縮的小波域方法[D].成都:電子科技大學,2004.)

[9] WELSTEAD S T.Fractal and Wavelet Image Compression Techniques[D].Washington:SPIE Optical Engineering Press, 1999,232.

[10] PRITCHARD M E,SIMONS M,ROSEN P A,et al.Coseismic Slip from the 1995 July 30 Mw＝8.1 Antofagasta,Chile,Earthquake as Constrained by InSAR and GPS Observations[J].Geophysical Journal International,2002,150(2):362-376.

[11] SIMONS M,FIALKO Y,RIVERA L.Coseismic Deformation from the 1999 Mw 7.1 Hector Mine,California,Earthquake as Inferred from InSAR and GPS Observations[J].Bulletin of the Seismological Society of America,2002,92: 1390-1402.

[12] LOH MAN R B,SIMONS M.Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation:Noise Structure and Data Downsampling[J].Geochemistiy,Geophysics,Geosystems,2005,Q01007.DOI: 10.1029/2004GC000841.

[13] NI Lin.Compression of Remote Sensing Images Based on Adaptive Quadtree Partitioning[J].Journal of Remote Sensing,2002,6(5):343-351.(倪林.基于自適應四叉樹分割的遙感圖像壓縮算法[J].遙感學報,2002,6(5): 343-351.)

[14] LIU Yang,GONG Adu,LI Jing.A Model for Massive 3D Terrain Simplification Based on Data Block Partition and Quad-tree[J].Acta Geodaetica et Cartographica Sinica, 2010,39(4):410-415.(劉揚,宮阿都,李京.基于數據分層分塊的海量三維地形四叉樹簡化模型[J].測繪學報, 2010,39(4):410-415.)

[15] HE Xingheng,CHEN Hui.Method Based on Quarter-tree Decomposition of Image Compression[J].Computer Engineering and Applications,2008,44(9):181-183.(何興恒,陳慧.一種基于四叉樹分解的圖像壓縮方法[J].計算機工程與應用,2008,44(9):181-183.)

[16] WAN Qing,ZHANG Lu,JIANG Houjun,et al.Estimation and Removal of Atmospheric Effects in InSAR Topographic Mapping[J].Journal of Remote Sensing,2012,16 (5):1074-1088.

[17] KNOSPE S H G,JONSSON S.Covariance Estimation for DInSAR Surface Deformation Measurements in the Pressence of Anisotropic Atmospheric Noise[J].IEEE Transactions on Geoscience and Remote Sensing,2010,48: 2057-2065.

[18] JI Lingyun,WANG Qingliang,CUI Duxin,et al.Influence of SAR Orbital Data to the Accuracy of DEM by Doris[J].Science of Surveying and Mapping,2009,34(5):57-59.(季靈運,王慶良,崔篤信,等.SAR衛(wèi)星軌道數據精度對Doris獲取DEM精度的影響研究[J].測繪科學,2009,34 (5):57-59.)

[19] JIN Guowang,WU Yirong,XIANG Maosheng,et al.Baseline Estimation Algorithm of InSAR with Block Adjustment [J].Acta Geodaetica et Cartographica Sinica,2011,40(5): 616-622.(靳國旺,吳一戎,向茂生,等.基于區(qū)域網平差的InSAR基線估計方法[J].測繪學報,2011,40(5): 616-622.)

[20] HETLAND E A,MUSéP,SIMONS M,et al.Multiscale InSAR Time Series(MIn TS)Analysis of Surface Deformation[J]Journal of Geophysical Research,2012,117, B02404.DOI:10.1029/2011JB008731.

(責任編輯:陳品馨)

A Quadtree InSAR Data Reduction Method Based on Covariance Function

ZHANG Jing1,ZHANG Qin1,2,ZHAO Chaoying1,2,ZHANG Juqing1

1.College of Geology Engineering and Geomatics,Chang’an University,Xi’an 710054,China;2.Key Laboratory of Western China’s Mineral Resources and Geological Engineering,Ministry of Education,Xi’an 710054,China

A major problem in inversion of deformation mechanism using InSAR data is that the InSAR results often contain thousands to millions of data points.Furthermore,there always exist errors and even some blunders,which make the data inversion be lower efficient and lower reliable.Thus,an adaptive quadtree decomposition method for InSAR data reduction is proposed in order to reduce the data numbers without losing the significant information about the deformation.The two important parameters of quadtree decomposition by covariance function are determined,which are eatablished by taking account of the physical spatial correlation of InSAR data.The algorithm can preserve details of deformation as much as possible and achieve efficient data reduction.This method is evaluated with InSAR data over Xi’an land subsidence.The results indicate that the algorithm proposed in this manuscript can not only reduce InSAR data number efficiently under a very good preservation of deformation signal,but can eliminate the noise of deformation results efficiently.

quadtree;data reduction;spatial correlation;covariance function

ZHANG Jing(1984—),female,PhD candidate,majors in InSAR data error analysis.

ZHANG Qin

P237

A

1001-1595(2014)05-0486-07

國家自然科學基金(41072266;41274004);國土資源大調查(1212011220186);中央高校基本科研業(yè)務費專項(2013G5260006);地震行業(yè)科研專項(201208009)

2012-12-21

張靜(1984—),女,博士生,研究方向為合成孔徑雷達干涉測量數據誤差分析。

E-mail:woshijing-001＠163.com

張勤

E-mail:zhangqinle＠263.net.cn

ZHANG Jing,ZHANG Qin,ZHAO Chaoying,et al.A Quadtree InSAR Data Reduction Method Based on Covariance Function [J].Acta Geodaetica et Cartographica Sinica,2014,43(5):486-492.(張靜,張勤,趙超英,等.顧及協(xié)方差函數的自適應四叉樹InSAR數據壓縮算法[J].測繪學報,2014,43(5):486-492.)

10.13485/j.cnki.11-2089.2014.0068

修回日期:2013-12-25