一種搜索魯棒的中國(guó)余數(shù)定理的多基線相位解纏繞技術(shù)*

蔣志彪，王 建，宋 千，周智敏

(國(guó)防科技大學(xué) 電子科學(xué)學(xué)院， 湖南 長(zhǎng)沙 410073)

合成孔徑雷達(dá)干涉測(cè)量(Interferometric Synthetic Aperture Radar, InSAR)是目前極具潛力的測(cè)繪技術(shù)之一，能夠獲得高分辨率和高精度的地形高程。同時(shí)多基線和多波段的InSAR技術(shù)進(jìn)一步提高了InSAR 技術(shù)獲取數(shù)字高程模型(Digital Elevation Model, DEM)的測(cè)量精度[1-2]。

目前，在多基線干涉合成孔徑雷達(dá)信號(hào)處理步驟中，多基線相位解纏繞仍是干涉處理的關(guān)鍵和核心問(wèn)題所在。在過(guò)去的幾十年里，許多的多基線相位解纏繞方法被先后提出來(lái)[3]。文獻(xiàn)[3]首次提出了一種基于聚類分析(Cluster Analysis, CA)的多基線相位解纏繞算法。文獻(xiàn)[4]提出了以最大似然和最小二乘思想為基礎(chǔ)的兩種相位解纏繞方法。文獻(xiàn)[5]提出了一種基線迭代的相位解纏繞算法，該方法通過(guò)不同基線之間的迭代，實(shí)現(xiàn)最終的相位解纏繞。1999年，F(xiàn)erretti等將小波技術(shù)運(yùn)用到相位解纏繞技術(shù)中[6]，提出基于小波思路的相位解纏繞算法。隨后，一些基于參數(shù)估計(jì)的多相位解纏繞方法也被相繼提出來(lái)，這些方法以統(tǒng)計(jì)學(xué)為架構(gòu)，如最大似然架構(gòu)[7-8]和最大后驗(yàn)架構(gòu)[9]，但是這類以統(tǒng)計(jì)學(xué)為架構(gòu)的多基線相位解纏繞方法都要求其干涉圖數(shù)據(jù)是獨(dú)立的，而且運(yùn)算耗時(shí)比較大。在文獻(xiàn)[9]的基礎(chǔ)上，F(xiàn)erraioli等為了進(jìn)一步降低運(yùn)算時(shí)間，提出了基于全變分(Total Variation, TV)模型和最優(yōu)圖割思路的優(yōu)化算法[10]。

文獻(xiàn)[11]首次將經(jīng)典中國(guó)余數(shù)定理(Classical Chinese Remainder Theorem, C-CRT)運(yùn)用到多基線InSAR相位解纏繞問(wèn)題中，后來(lái)文獻(xiàn)[12]對(duì)基于C-CRT的多基線相位解纏繞方法進(jìn)行了應(yīng)用推廣，另外文獻(xiàn)[13]將這一方法運(yùn)用到了機(jī)載毫米波InSAR多基線相位解纏繞中，分析了基線構(gòu)型誤差和基線構(gòu)型優(yōu)化設(shè)計(jì)等問(wèn)題。無(wú)疑，基于C-CRT的相位解纏繞方法較好解決了干涉相位欠采樣處的解纏難題，然而該方法對(duì)噪聲特別敏感，而實(shí)際中獲取的干涉數(shù)據(jù)又不可避免地存在噪聲，此問(wèn)題的存在極大地限制了該方法在實(shí)際中的廣泛運(yùn)用。

值得注意的是，Xia等在中國(guó)余數(shù)定理方面進(jìn)行了深入的研究，提出了一種魯棒的中國(guó)余數(shù)定理[14]，該方法通過(guò)搜索的方式來(lái)確定模糊數(shù)，故又稱之為搜索魯棒的中國(guó)余數(shù)定理(Searched-form Robust Chinese Remainder Theorem, SR-CRT)。該方法利用余數(shù)的冗余來(lái)魯棒地恢復(fù)原來(lái)的數(shù)，因而具有較好的魯棒性，故該方法廣泛地運(yùn)用于動(dòng)目標(biāo)檢測(cè)、測(cè)距和頻率估計(jì)等問(wèn)題中[15-18]。

本文通過(guò)創(chuàng)新性地引入公因子，構(gòu)建了新的同余方程組，將SR-CRT首次運(yùn)用到多基線相位解纏繞中，并重點(diǎn)分析了其抗噪聲性能。

1 中國(guó)余數(shù)定理原理

1.1 經(jīng)典中國(guó)余數(shù)定理

設(shè)兩兩互質(zhì)的正整數(shù)Γ1，Γ2，…，ΓL，且都大于1；a1，a2，…，aL為任意整數(shù)，n1,n2,…,nL為模糊數(shù)，一次同余方程組：

(1)

在0≤x<Γ?Γ1Γ2…ΓL時(shí)，同余方程組有唯一解x，解的形式可以寫(xiě)成如下形式：

(2)

(3)

1.2 搜索魯棒的中國(guó)余數(shù)定理

設(shè)正整數(shù)M1，M2，…，ML存在最大公約數(shù)C，且都大于1；r1，r2，…，rL為任意整數(shù)，一次同余方程組：

(4)

式中，正整數(shù)Mi滿足

Mi=CΓii=1,2,…,L

(5)

考慮噪聲時(shí)，即余數(shù)中存在誤差，如

(6)

其中，Δri是余數(shù)誤差。

(7)

式中，[·]表示四舍五入運(yùn)算。模糊數(shù)ni是通過(guò)式(8)中的二維搜索方式求解，故稱為搜索魯棒的中國(guó)余數(shù)定理(SR-CRT)。

(8)

其中，

2≤i≤L

(9)

通過(guò)上述兩種定理原理的對(duì)比可知， C-CRT要求同余方程組的模兩兩互質(zhì)，而SR-CRT要求模之間具有最大公約數(shù)，且除去該最大公約數(shù)后所得的數(shù)都兩兩互質(zhì)。

顯然，兩種方法對(duì)運(yùn)用的同余方程組的模具有不同的要求，對(duì)于滿足C-CRT進(jìn)行相位解纏繞要求的基線數(shù)據(jù)，是否也可以通過(guò)SR-CRT進(jìn)行相位解纏繞呢？

2 基于SR-CRT的相位解纏繞技術(shù)

對(duì)于上述問(wèn)題，可通過(guò)人為引入一個(gè)公因子，構(gòu)建新的同余方程組，使之滿足SR-CRT的條件，下面以兩幅干涉圖的情形詳細(xì)介紹其過(guò)程。

2.1 相位解纏繞原理

根據(jù)多基線干涉幾何關(guān)系，不同基線長(zhǎng)度下的兩幅干涉圖中，每一個(gè)像素都近似滿足[3]：

(10)

式中：φ1，φ2為不同基線長(zhǎng)度下的纏繞相位(去平地相位)；B1，B2為不同的垂直基線長(zhǎng)度；n1，n2分別代表兩幅干涉相位圖的模糊數(shù)。

利用C-CRT進(jìn)行相位解纏的條件是兩幅干涉圖中的基線長(zhǎng)度比值必須互質(zhì)，即滿足

(11)

就可以構(gòu)建出同余方程組：

(12)

根據(jù)定義可得余數(shù)ai和剩余數(shù)ei分別為：

(13)

其中，?·」為取整運(yùn)算。

對(duì)于滿足C-CRT進(jìn)行相位解纏的基線條件，如式(11)所示，將分子和分母同時(shí)乘以一個(gè)公因子C，可得：

(14)

聯(lián)立式(10)和式(14)，得：

M1(φ1+2πn1)=M2(φ2+2πn2)

(15)

將式(15)中各項(xiàng)提出2π，變形可得：

(16)

分別定義新的余數(shù)ri和剩余數(shù)fi為：

(17)

將式(17)的余數(shù)和剩余數(shù)代入式(16)可得：

2π(f1+r1+M1n1)=2π(f2+r2+M2n2)

(18)

將式(18)中各項(xiàng)對(duì)2π進(jìn)行求余，可得剩余數(shù)f1和f2是相等的，最后可以構(gòu)建出新的同余方程組：

(19)

只要公因子C取任意正整數(shù)，就可以利用SR-CRT求解模糊數(shù)n1，n2，得出方程組的解。

對(duì)i=1,2，考慮噪聲時(shí)，即余數(shù)中存在誤差，如下式

(20)

式中，Δri為余數(shù)誤差。

當(dāng)余數(shù)誤差滿足|Δri|≤τ

φi=φi+2niπ

(21)

以上就是基于SR-CRT的多基線相位解纏繞的基本原理。

2.2 抗噪聲性能分析

通過(guò)相位解纏原理分析可知，引入公因子后，新的同余方程組能夠運(yùn)用SR-CRT，那么其抗噪聲性能如何？

當(dāng)纏繞相位中存在噪聲的情況下，可以得出：

(22)

式中，Δφi是噪聲相位，Δφi∈[0,2π)。

根據(jù)式(17)和式(22)，可以得到有噪聲的余數(shù)表達(dá)式為：

(23)

又0

(24)

根據(jù)SR-CRT可知，當(dāng)余數(shù)誤差Δri滿足|Δri|≤τ

(25)

式(25)就是基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能表達(dá)式，其抗噪聲性能受參數(shù)Γi、剩余數(shù)fi、公因子C的影響。

需要注意的是，由于公因子C可以取任意的正整數(shù)，又0

(26)

此時(shí)噪聲相位Δφi只需要滿足：

0≤Δφi<π/Γi

(27)

這時(shí)，基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能僅受參數(shù)Γi的影響，滿足式(27)條件的噪聲相位的區(qū)間分布示意如圖1所示。

通過(guò)噪聲性能原理分析可知，基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能受參數(shù)Γi、剩余數(shù)fi、公因子C的影響，而通過(guò)選擇合適公因子C可以消除剩余數(shù)fi、公因子C對(duì)噪聲相位的限制，從而提高該方法的抗噪聲性能。

圖1 噪聲相位分布Fig.1 Distribution of phase noise

2.3 多維情況推廣

上述相位解纏原理是以兩幅干涉圖為例分析的，本節(jié)將給出其多維情況的推廣。

根據(jù)多基線干涉幾何關(guān)系，對(duì)于L幅不同基線的干涉相位圖，同理，每一個(gè)像素也都滿足[3]：

(28)

式中：φ1，φ2，…，φL為不同基線長(zhǎng)度下的纏繞相位(去平地相位)；B1，B2，…，BL為不同垂直基線長(zhǎng)度；n1，n2，…，nL分別代表在不同干涉相位圖中的模糊數(shù)。

L幅干涉圖中的任意兩副基線的長(zhǎng)度比值必須互質(zhì)，即滿足：

i=1,2,…,L，j=1,2,…,L,i≠j

(29)

同樣，將式(29)中的分子和分母同時(shí)乘以一個(gè)公因子C，構(gòu)建出新的同余方程組：

(30)

定義新的余數(shù)ri和剩余數(shù)fi為：

(31)

(32)

只要公因子C為正整數(shù)，就可以利用SR-CRT求解模糊數(shù)，最后求出絕對(duì)相位。

通過(guò)上述基于SR-CRT的多基線相位解纏繞的基本原理分析可知，無(wú)論是二維還是多維情況，該方法都需要引入公因子C，而且公因子C可取任意正整數(shù)。下面通過(guò)試驗(yàn)分析公因子C的取值對(duì)基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能的影響。

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

在上一節(jié)中，詳細(xì)分析了基于SR-CRT方法的相位解纏原理和抗噪聲性能，為了進(jìn)一步驗(yàn)證該方法的有效性和抗噪聲性能，本節(jié)進(jìn)行性能驗(yàn)證實(shí)驗(yàn)和算法執(zhí)行時(shí)間分析，并與基于C-CRT方法的相位解纏結(jié)果和文獻(xiàn)[3]中基于CA方法的相位解纏繞結(jié)果進(jìn)行對(duì)比分析，實(shí)驗(yàn)中驗(yàn)證算法的軟件運(yùn)行環(huán)境為MATLAB2012，計(jì)算機(jī)的CPU主頻為3.2 GHz，內(nèi)存為8 GB。

本實(shí)驗(yàn)采用由美國(guó)Long′s Peak國(guó)家公園真實(shí)數(shù)字地形高程圖(圖像大小為152×458)生成的干涉相位圖來(lái)驗(yàn)證算法，實(shí)驗(yàn)以三幅干涉圖為例來(lái)進(jìn)行分析，主要系統(tǒng)參數(shù)如下：雷達(dá)系統(tǒng)天線頻段為C波段，S1、S2、S3和S4分別表示四副自發(fā)自收天線的中心位置，以天線S1為參考，即S1為主天線，雷達(dá)主天線到場(chǎng)景中心的斜距為570 km，主天線的入射角為10.5°。三組基線的構(gòu)建形式為：基線B1代表天線S1和S2對(duì)應(yīng)基線的基線長(zhǎng)度，基線B1長(zhǎng)度為120 m；基線B2代表天線S1和S3對(duì)應(yīng)基線的基線長(zhǎng)度，基線B2長(zhǎng)度為150 m；基線B3代表天線S1和S4對(duì)應(yīng)基線的基線長(zhǎng)度，基線B3的長(zhǎng)度為200 m。

圖2(a)～(c)分別為基線B1、B2、B3對(duì)應(yīng)的絕對(duì)相位參考圖。圖2(d)～(f)分別為基線B1、B2、B3對(duì)應(yīng)的纏繞相位圖，從圖中可以看出基線越長(zhǎng)，干涉條紋越密集，其中，纏繞相位圖中的相位噪聲為改進(jìn)的高斯噪聲，考慮了地形和視數(shù)因素。

為了便于分析基于SR-CRT算法的性能，給出了基于C-CRT方法的相位解纏繞結(jié)果和基于CA方法的相位解纏繞結(jié)果做對(duì)比，并給出了公因子分別取值為6和200時(shí)，基于SR-CRT方法的相位解纏繞結(jié)果，下面以B3基線的相關(guān)解纏結(jié)果為例進(jìn)行分析。

圖3是基于C-CRT的多基線相位解纏繞技術(shù)的求解結(jié)果。從相位誤差圖中，可以明顯看出解纏繞出現(xiàn)錯(cuò)誤區(qū)域比較密集，并且成條紋形狀，通過(guò)統(tǒng)計(jì)直方圖可知，絕對(duì)相位誤差均值為1.45，標(biāo)準(zhǔn)差為30.60，絕對(duì)相位重建誤差大于2π的像素占總像素的28.22%。

圖4是基于SR-CRT方法、公因子C取值為6時(shí)的求解結(jié)果。通過(guò)絕對(duì)相位圖和相位誤差圖可知，圖像兩側(cè)區(qū)域求解絕對(duì)相位與參考相位相比偏差較大，主要原因是這些區(qū)域地形比較陡峭，造成圖像的相干性很低，即在圖像表現(xiàn)出比較大的噪聲，使式(25)中條件不滿足，從而無(wú)法通過(guò)SR-CRT搜索到正確的模糊數(shù)，造成解纏繞結(jié)果出現(xiàn)錯(cuò)誤。通過(guò)誤差統(tǒng)計(jì)直方圖可知，絕對(duì)相位重建誤差均值為0.19，標(biāo)準(zhǔn)差為10.89，其中重建誤差大于2π的像素占總像素的5.48%。

圖5是基于SR-CRT方法、公因子C取值200時(shí)的求解結(jié)果，通過(guò)絕對(duì)相位圖和相位誤差圖可知，除了圖像兩側(cè)區(qū)域孤立的噪聲點(diǎn)外，絕對(duì)相位圖與參考相位圖基本相似。通過(guò)誤差統(tǒng)計(jì)直方圖可知，誤差均值為0.15，標(biāo)準(zhǔn)差為7.92，其中重建誤差大于2π的像素占總像素的1.88%。

圖6是基于CA方法的求解結(jié)果，通過(guò)絕對(duì)相位圖和相位誤差圖可知，除了圖像兩側(cè)區(qū)域孤立的噪聲點(diǎn)外，絕對(duì)相位圖與參考相位圖基本相似。通過(guò)誤差統(tǒng)計(jì)直方圖可知，誤差均值為0.06，標(biāo)準(zhǔn)差為8.85，其中重建誤差大于2π的像素占總像素的2.32%。

通過(guò)圖3、圖4、圖5的實(shí)驗(yàn)結(jié)果對(duì)比可知，基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能明顯優(yōu)于基于C-CRT的相位解纏繞技術(shù)；通過(guò)圖4與圖5的實(shí)驗(yàn)結(jié)果對(duì)比可知，公因子C取值200時(shí)的求解結(jié)果要優(yōu)于公因子C取值6時(shí)的求解結(jié)果，因而選擇較大的公因子可進(jìn)一步提高抗噪聲性能。

圖6 基于CA方法的解纏結(jié)果Fig.6 Results of the CA method

通過(guò)圖4、圖5和圖6的實(shí)驗(yàn)結(jié)果對(duì)比分析可知，基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能比較接近于基于CA相位解纏繞技術(shù)的，而且當(dāng)公因子C取值合理時(shí)，基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能甚至要優(yōu)于基于CA相位解纏繞技術(shù)的。

下面對(duì)基于SR-CRT的多基線相位解纏技術(shù)(簡(jiǎn)稱搜索方法)的算法執(zhí)行時(shí)間進(jìn)行仿真分析。由原理分析可知該方法是以單像素為單位進(jìn)行相位解纏繞的，同基于C-CRT的多基線相位解纏技術(shù)(簡(jiǎn)稱經(jīng)典方法)一樣，都屬于單像素相位解纏繞。而文獻(xiàn)[3]中基于CA的多基線相位解纏技術(shù)是以像素類為單位進(jìn)行相位解纏(簡(jiǎn)稱聚類分析方法)。圖7給出了這三種方法的執(zhí)行時(shí)間跟像素?cái)?shù)目的變化關(guān)系。

圖7 三種算法耗時(shí)Fig.7 Time consuming of three algorithms

從圖7中可知，三種方法的執(zhí)行時(shí)間跟像素?cái)?shù)目是成正比的。通過(guò)三種方法的執(zhí)行時(shí)間結(jié)果對(duì)比可知，聚類分析方法耗時(shí)最優(yōu)，其執(zhí)行時(shí)間隨著圖像尺寸增大而緩慢增加，相對(duì)而言，經(jīng)典方法和搜索方法執(zhí)行時(shí)間隨著圖像尺寸增大而快速增加。造成上述結(jié)果的主要原因是聚類分析方法是以像素類為單位進(jìn)行解纏繞的，受圖像尺寸的影響相對(duì)較小，而經(jīng)典方法和搜索方法都是以單個(gè)像素為單位進(jìn)行相位解纏繞的，因而隨著圖像尺寸增大，算法的執(zhí)行時(shí)間明顯增加。另外，搜索方法的耗時(shí)要比經(jīng)典方法的耗時(shí)長(zhǎng)，其主要原因是搜索方法需要搜索模糊數(shù)，每一個(gè)像素都需要進(jìn)行式(8)這樣的二維搜索，盡管文獻(xiàn)[19-20]已經(jīng)將式(8)這樣的二維搜索過(guò)程簡(jiǎn)化成了一維搜索，但算法的計(jì)算復(fù)雜度在搜索次數(shù)較大的情況下仍然比較高，其耗時(shí)隨著圖像尺寸增大而急劇增大。

4 結(jié)論

將SR-CRT首次運(yùn)用到多基線相位解纏中，重點(diǎn)分析和研究了其相位解纏繞原理和抗噪聲性能，最后通過(guò)實(shí)驗(yàn)驗(yàn)證該方法的有效性，并與基于C-CRT的多基線相位解纏繞技術(shù)進(jìn)行了對(duì)比分析，總結(jié)得出了以下四點(diǎn)重要結(jié)論：

1)與基于C-CRT的相位解纏繞技術(shù)相比，基于SR-CRT的相位解纏繞技術(shù)中，通過(guò)公因子這一可控因素的引入(公因子可控是指公因子可取任意正整數(shù))，可以基本消除剩余數(shù)對(duì)噪聲相位的限制，從而提高算法的抗噪聲性能。

2)基于SR-CRT的相位解纏繞技術(shù)的抗噪聲性能明顯優(yōu)于基于C-CRT的相位解纏繞技術(shù)的抗噪聲性能，而通過(guò)選擇較大的公因子可進(jìn)一步提高基于SR-CRT的相位解纏繞抗噪聲性能。

3)對(duì)于公因子C的取值，結(jié)合式(25)和式(26)可知，單從理論角度來(lái)說(shuō)，公因子C取值越大(即正無(wú)窮大)越好，然而在實(shí)際過(guò)程中，公因子C取值可以根據(jù)式(26)中的第二個(gè)條件來(lái)確定，一般取滿足條件的所有C取值中的較小值就可以使得基于SR-CRT的相位解纏繞技術(shù)獲得比較滿意的抗噪聲性能。

4)基于SR-CRT的相位解纏繞技術(shù)是以單個(gè)像素為單位進(jìn)行相位解纏繞的，對(duì)于每一個(gè)像素，都需要通過(guò)搜索來(lái)求解模糊數(shù)，算法的復(fù)雜度和搜索的次數(shù)有關(guān)，當(dāng)搜索次數(shù)較大的情況下，其算法耗時(shí)隨著像素?cái)?shù)目增大而急劇增加。

上述抗噪聲性能分析和結(jié)論為基于SR-CRT的相位解纏繞技術(shù)進(jìn)一步在實(shí)際中的運(yùn)用奠定了理論基礎(chǔ)。后續(xù)研究將著重探討對(duì)含有噪聲的實(shí)測(cè)多基線干涉數(shù)據(jù)進(jìn)行相位解纏繞以及進(jìn)一步有效降低基于SR-CRT的相位解纏繞技術(shù)耗時(shí)的方法。