張金翼，王 輝，吳思利，鄭世超，顧約翰

（1.上海衛(wèi)星工程研究所，上海 201109；2.上海市毫米波空天信息獲取及應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室，上海 201109）

0 引言

近年來，地下資源的過度開采、水土流失等問題愈發(fā)嚴(yán)重，地質(zhì)災(zāi)害（如地震、泥石流、礦區(qū)塌陷、城市地面沉降）頻發(fā)，對(duì)我國(guó)人民群眾的生命及財(cái)產(chǎn)安全造成威脅。在該類災(zāi)害中典型的表征為地表形變，因此監(jiān)測(cè)地表形變是防災(zāi)減災(zāi)的重要手段。傳統(tǒng)的測(cè)量地表形變技術(shù)手段為水準(zhǔn)測(cè)量、GPS 測(cè)量等［1］，但上述方法存在費(fèi)時(shí)、費(fèi)力、數(shù)據(jù)點(diǎn)孤立的缺陷，無法及時(shí)、高效地預(yù)測(cè)和評(píng)估地質(zhì)災(zāi)害。差分干涉合成孔徑雷達(dá)（Differential Interferometric Synthetic Aperture Radar，DInSAR）是天基遙感中較為常用的測(cè)量地表形變的技術(shù)手段，因?yàn)槿鞎r(shí)、全天候工作，速度快、范圍廣、測(cè)量密度高等特點(diǎn)逐漸受到人們的重視［2-3］，差分干涉合成孔徑雷達(dá)形變測(cè)量主要利用相位的變化，直接測(cè)得的干涉相位被包裹在［-Π，Π）)[-π，π 內(nèi)，要獲得可用的相位信號(hào)須對(duì)其進(jìn)行解纏繞，其效果對(duì)測(cè)量的精度有較大的影響［4-5］。

常用的解纏方法按照原理分為以下3 類：1）基于路徑跟蹤解纏，選擇不同的路徑對(duì)相鄰像元進(jìn)行差值積分恢復(fù)真實(shí)相位，較有代表性的是GOLSTEIN［6］于1988 年提出的枝切法，可在計(jì)算出殘差點(diǎn)后，快速有效地尋找連接相鄰殘差點(diǎn)的最優(yōu)路徑；2）基于最小范數(shù)思想，通過擬合函數(shù)求得纏繞相位和解纏相位的離散偏微分差，進(jìn)而求解相位解纏的整體最佳估計(jì)值，其中較為常用的是最小二乘法、Jacobi 迭代法、Gauss-Seidel 迭代法和逐次超松弛（Successive Over Relaxation，SOR）迭代法，這些算法的缺陷在于數(shù)據(jù)量較大時(shí)收斂速度慢，不適用于大規(guī)模數(shù)據(jù)［7-8］；3）基于網(wǎng)絡(luò)規(guī)劃的算法，由COSTANTINI［9］于1996 年提出，網(wǎng)絡(luò)規(guī)劃算法 在相位解纏中引入最小代價(jià)流的概念，通過搜索全區(qū)域最短路徑生成最優(yōu)枝切線進(jìn)行解纏，可以將誤差限制在相干質(zhì)量差的區(qū)域內(nèi)，保證其余區(qū)域結(jié)果的準(zhǔn)確性，精度較高，因此常被用作其他算法的精度參照，但由于最小費(fèi)用網(wǎng)絡(luò)流（Minimum Cost Network Flow，MCF）的算法特性，當(dāng)相位差圖像相干性差、殘差點(diǎn)數(shù)量多時(shí)，算法的復(fù)雜度較高，效率較低［10-12］。CHEN［13］提出改進(jìn)的MCF 算法，分割整個(gè)圖像，與傳統(tǒng)的MCF 算法相比，提高了算法效率，但分塊的規(guī)模為人工設(shè)定，沒有具體的標(biāo)準(zhǔn)，不同的分塊大小會(huì)對(duì)該算法的精度產(chǎn)生較大影響。

本文提出一種基于自適應(yīng)分塊的改進(jìn)MCF 解纏算法，該方法根據(jù)質(zhì)量圖所表現(xiàn)的殘差信息，將分塊大小進(jìn)行自適應(yīng)尋優(yōu)，更多地將高相干區(qū)域劃分在一個(gè)子塊。實(shí)驗(yàn)結(jié)果表明，該算法在提升算法效率的同時(shí)保證了準(zhǔn)確度，有效地拓展了MCF 算法的應(yīng)用場(chǎng)景。

1 MCF 算法解纏原理

相位解纏的基本思路是將已獲取的干涉圖中的纏繞相位，通過疊加纏繞數(shù)恢復(fù)相位的真實(shí)值，進(jìn)而反映真實(shí)的形變情況［14-15］。

定義纏繞相位為ψi，j，真實(shí)相位為?i，j，對(duì)于一個(gè)M×N的矩陣，可得方程式如下：

式中：k為整數(shù)，0 ≤i≤M-1，0 ≤J≤N-1。

定義纏繞算子W，對(duì)相鄰像素間的纏繞相位進(jìn)行±2kπ 的操作，使其變?yōu)檎鎸?shí)相位。定義x方向和y方向上的絕對(duì)相位差分別為，其方程式如下：

式中：ψi+1，j為M×N矩陣中(i+1，j)處點(diǎn)的纏繞相位，rad；ψi，j+1為M×N矩陣中(i，j+1)處點(diǎn)的纏繞相位，rad。

由式（2）變換可得方程式如下：

對(duì)于干涉圖中大部分點(diǎn)，沿著點(diǎn)(i，j)、(i+1，j)、(i+1，j+1)、(i，j+1)的纏繞相位差環(huán)路積分為0，而實(shí)際應(yīng)用中，由于噪聲等因素的存在導(dǎo)致出現(xiàn)殘差點(diǎn)，殘差點(diǎn)處的纏繞相位差環(huán)路積分為±2π，即定義殘差為

式 中：ei，j為殘差；為(i+1，j)點(diǎn)在y方向上的絕對(duì)相位差，單位為rad；為(i，j+1)點(diǎn) 在x方向上的絕對(duì)相位差，單位為rad。

由于存在殘差點(diǎn)，無法確定干涉圖全圖的積分，其變成與積分路徑有關(guān)的變量，因此需考慮殘差點(diǎn)的影響，使其降至最小。在最小費(fèi)用流解纏算法中，將正/負(fù)殘差點(diǎn)定義為供應(yīng)/需求節(jié)點(diǎn)，殘差為供應(yīng)/需求量，將其他點(diǎn)定義為轉(zhuǎn)運(yùn)節(jié)點(diǎn)。相鄰的節(jié)點(diǎn)之間通過帶有流的弧連接，以調(diào)節(jié)供需平衡，最小費(fèi)用流解纏的目的是在所有節(jié)點(diǎn)供需平衡的基礎(chǔ)上，確?？傎M(fèi)用取最小。達(dá)到供需平衡即［16-18］：

式 中：為 點(diǎn)(i，j)到點(diǎn)(i+1，j)的流量；為點(diǎn)(i，j)到點(diǎn)(i，j+1)的流量，1 ≤i≤M-1，1 ≤j≤N-1。

總費(fèi)用流為

式中：為點(diǎn)(i，j)到點(diǎn)(i+1，j)對(duì)應(yīng)的費(fèi)用；為點(diǎn)(i，j)到點(diǎn)(i，j+1)的費(fèi)用。

式中：ψ0，0為（0，0）位置處點(diǎn)的纏繞相位，rad。

綜上所述，應(yīng)用MCF 算法的解纏步驟如下。

1）對(duì)相干圖做預(yù)處理，選定閾值，提取相干系數(shù)高于閾值的相位。

2）在相位的集合中建立Delaunay 三角網(wǎng)，進(jìn)而構(gòu)建對(duì)偶圖，將殘差映射到對(duì)偶網(wǎng)絡(luò)中。

3）在對(duì)偶網(wǎng)絡(luò)中，應(yīng)用最小費(fèi)用流法連接正負(fù)殘差點(diǎn)對(duì)，計(jì)算最小費(fèi)用流集合。

4）根據(jù)流的大小和方向?qū)ο辔痪仃嚪e分，得到解纏結(jié)果，再從高質(zhì)量區(qū)域向低質(zhì)量區(qū)域積分［12］。

2 自適應(yīng)分塊的改進(jìn)MCF 解纏算法

MCF 算法一經(jīng)推出，其準(zhǔn)確度便得到肯定，但隨著殘差點(diǎn)的增多，該算法的效率有所下降，因此提高算法效率是優(yōu)化該算法的一個(gè)方向［16］。CHEN［13］提出改進(jìn)的MCF 算法，將整個(gè)相干圖劃分為不重疊的矩形塊，將矩形塊作為完全獨(dú)立的圖像進(jìn)行相位解纏，降低解纏過程中所需的內(nèi)存資源，提高效率。

原始相干圖使用MCF 算法解纏的算法復(fù)雜度可表示為O(N2)，其中N為圖中殘點(diǎn)的數(shù)量，將相干圖分為i個(gè)矩形塊，則分割后算法復(fù)雜度可表示為O(N12+N22+N32+…+Ni2)，其 中N1+N2+N3+…+Ni=N，后者的算法復(fù)雜度遠(yuǎn)小于前者。

該種方法雖然提升了相位解纏的效率，但由于分割大小的不同會(huì)導(dǎo)致解纏精度的不穩(wěn)定。在本文中引入一種基于質(zhì)量圖的自適應(yīng)分塊改進(jìn)MCF算法，切割時(shí)，根據(jù)質(zhì)量圖自適應(yīng)選取分割大小，保留相干性較強(qiáng)的區(qū)域，在提高效率的同時(shí)保證精度。

本文提出的算法流程如下。

1）對(duì)相干圖做預(yù)處理，選定閾值，提取相干系數(shù)高于閾值的相位。

2）根據(jù)相干系數(shù)的計(jì)算方式計(jì)算質(zhì)量圖。

在干涉SAR 中殘差點(diǎn)的存在是由于SAR 疊掩或去相關(guān)噪聲引起的，因此相關(guān)系數(shù)較低的區(qū)域其殘差點(diǎn)較多，相關(guān)系數(shù)較高的區(qū)域殘差點(diǎn)較少。一般將表征干涉圖像各點(diǎn)相位的特征圖稱為質(zhì)量圖，質(zhì)量圖可劃分為相干系數(shù)圖、偽相干系數(shù)圖、相位導(dǎo)數(shù)變化圖及最大相位梯度圖4 種，本文采用使用范圍較廣的相干系數(shù)圖［20］。

相干系數(shù)圖中各個(gè)區(qū)域之間的相干性通過相干系數(shù)表示，相干系數(shù)值較高表示2 個(gè)區(qū)域之間的相關(guān)性較好，即相位信息準(zhǔn)確；相干系數(shù)值較低說明2 個(gè)區(qū)域之間相關(guān)性較差，即相位信息有誤。A、B2 張復(fù)圖像的相干系數(shù)可表示為

式中：E為求期望；γ為相干系數(shù)；B*為B的共軛。

3）根據(jù)質(zhì)量圖選定應(yīng)進(jìn)行分割的塊的尺寸。

質(zhì)量圖表示相干系數(shù)的大小，在最小費(fèi)用流解纏方法中，影響解纏效率與準(zhǔn)確度的是干涉圖中集中存在的殘差點(diǎn)。因此在分割時(shí)，將相干性較強(qiáng)的區(qū)域集中在一起，進(jìn)而使得殘差點(diǎn)被零散的分割到各個(gè)子塊內(nèi)，提高每個(gè)分割塊的解纏效率。

在實(shí)際操作過程中，為濾去噪聲求整張質(zhì)量圖的均值，將其作為門限濾波，表達(dá)式為

式中：T為門限值；M、N為質(zhì)量圖規(guī)模；f(i，j)為(i，j)位置點(diǎn)的相干系數(shù)。

在選擇分塊尺寸的過程中，選用密度峰值聚類算法。該算法是目前較為先進(jìn)的聚類算法，需要較少的輸入?yún)?shù)，能檢測(cè)聚類中心，具有簡(jiǎn)單、聚類速度快、穩(wěn)定性好、聚類高效的優(yōu)點(diǎn)，且對(duì)噪聲不敏感，是較為理想的聚類算法［21-22］。

對(duì)于質(zhì)量圖中的數(shù)據(jù)，采用歐氏距離計(jì)算兩點(diǎn)之間的相似性：

式中：(xi，xj)為第i個(gè)點(diǎn)和第j個(gè)點(diǎn)；l為維度；m為維度的數(shù)量。

計(jì)算數(shù)據(jù)集中每個(gè)數(shù)據(jù)點(diǎn)的密度，表達(dá)式如下：

式中：S(m)為判斷函數(shù)；dc為截?cái)嗑嚯x，設(shè)定為5。

當(dāng)m≥0 時(shí)S(m)=1；當(dāng)m≤0 時(shí)S(m)=0。

依照式（11），即可計(jì)算得到以每一點(diǎn)為中心的與其距離小于5 的點(diǎn)的個(gè)數(shù)，再求取每一數(shù)據(jù)點(diǎn)與其他密度比其大的數(shù)據(jù)點(diǎn)間的最小距離，如該點(diǎn)密度最大，則取與其他數(shù)據(jù)點(diǎn)距離的最大值為

式中：di，j為第i個(gè)點(diǎn)與第j個(gè)點(diǎn)的距離。

求得密度及最小距離后，將橫軸表示數(shù)據(jù)點(diǎn)的密度，縱軸表示數(shù)據(jù)點(diǎn)的最小距離，選取密度高且距離大的點(diǎn)作為聚類中心點(diǎn)，再將其他點(diǎn)按密度值降序排列，歸于聚類中心點(diǎn)［23-25］。

得到聚類結(jié)果后，將聚類點(diǎn)所在的矩形區(qū)域作為參考的干涉圖分塊區(qū)域，如聚類后產(chǎn)生多個(gè)子塊，將其大小取均值作為最終的分塊參考大小。

依照此方法得到的分塊可將相干性較強(qiáng)的點(diǎn)集中在同一塊區(qū)域，并依照大小劃分其他區(qū)域，提高算法的效率和準(zhǔn)確度。

4）在相位的集合中建立Delaunay 三角網(wǎng)，進(jìn)而構(gòu)建對(duì)偶圖，將殘差映射到對(duì)偶網(wǎng)絡(luò)中。

5）在對(duì)偶網(wǎng)絡(luò)中應(yīng)用最小費(fèi)用流法連接正負(fù)殘差點(diǎn)對(duì)，計(jì)算最小費(fèi)用流集合。

6）根據(jù)流的大小和方向?qū)ο辔痪仃嚪e分，得到解纏結(jié)果，再從高質(zhì)量區(qū)域向低質(zhì)量區(qū)域積分［12］。

算法流程如圖1 所示。

圖1 自適應(yīng)分塊的改進(jìn)MCF 解纏算法流程Fig.1 Flow chart of the improved MCF unwrapping algorithm with adaptive chunking

3 仿真實(shí)驗(yàn)及結(jié)果分析計(jì)算機(jī)模擬結(jié)果

仿真實(shí)驗(yàn)在一臺(tái)配備Intel Core i7-10750H 處理器（2.6 GHz，6 核12 線程）、NVIDIA GeForce GTX 1660 Ti GPU（6 GB 顯存）和16 GB DDR4 內(nèi)存的計(jì)算機(jī)上進(jìn)行，運(yùn)行于Windows 10 操作系統(tǒng)，使用Matlab R2019b 仿真，在該平臺(tái)生成干涉條紋圖，并按照本文方法對(duì)其解纏繞，解纏前干涉如圖2所示。

圖2 解纏前干涉Fig.2 Interference image before phase unwrapping

3.1 仿真實(shí)驗(yàn)

上文中所提到的方法對(duì)干涉圖像做處理，經(jīng)自適應(yīng)算法選取可知高相干點(diǎn)位于（80～170，50～180），塊大小為90×130，將上圖保留高相干塊并依照其大小劃分其他區(qū)域，進(jìn)而對(duì)其解纏繞，與從小到大隨機(jī)選取其他分塊大小做解纏對(duì)比，解纏后實(shí)驗(yàn)結(jié)果如圖3 所示。圖3（a）～圖3（f）分別為尋優(yōu)后分塊的解纏結(jié)果、選取分割塊大小為50×50 的解纏結(jié)果、選取分割塊大小為100×100 的解纏結(jié)果、選取分割塊大小為150×150 的解纏結(jié)果、選取分割塊大小為200×200 的解纏結(jié)果和選取分割塊大小為256×256 的解纏結(jié)果。圖3（a）～圖3（f）的解纏全流程處理時(shí)間分別為4.07、2.18、3.32、4.55、4.51 和8.92 s。

圖3 不同分塊方式解纏后圖像Fig.3 Images after unwrapping with different chunking methods

不同分塊解纏結(jié)果誤差如圖4 所示。其中圖4（a）～圖4（f）的解纏后均方根誤差（Root Mean Square Error，RMSE）分別為4.002 8、5.564 1、5.574 6、5.232 5、5.232 5 和3.962 9 rad。

圖4 不同分塊方式解纏圖像后均方根誤差Fig.4 RMSE results after unwrapping the images with different chunking methods

仿真試驗(yàn)結(jié)果統(tǒng)計(jì)見表1。

表1 仿真試驗(yàn)結(jié)果Tab.1 Simulation results

3.2 結(jié)果分析

由上述仿真結(jié)果可知，采用文中算法選取到的分割大小可在4.07 s 內(nèi)完成解纏，誤差為4.002 8 rad，對(duì)比其他隨機(jī)選取的分塊大小，本算法在保證解纏效率的同時(shí)，其準(zhǔn)確度沒有大幅下滑。對(duì)比誤差最小的分塊方式（256×256），解纏誤差惡化了1.01%，但效率提升了54.37%。證明本算法在平衡效率和準(zhǔn)確度方面具有一定的優(yōu)勢(shì)，對(duì)提升改進(jìn)MCF 算法的應(yīng)用場(chǎng)景與計(jì)算效率有一定的貢獻(xiàn)。

4 結(jié)束語

針對(duì)改進(jìn)的基于分割的最小費(fèi)用網(wǎng)絡(luò)流算法準(zhǔn)確率及效率受分割塊大小影響的問題，提出利用質(zhì)量圖表征相干性優(yōu)劣的特性自適應(yīng)確定分割塊大小的方法，進(jìn)而將相干性較好的區(qū)域聚集在一起，將殘差點(diǎn)散落地分布在不同分割塊中。仿真結(jié)果表明，本方法通過自適應(yīng)尋優(yōu)方式選取合適的分割塊，使得準(zhǔn)確度在不受過多影響的情況下，有效地提高算法效率。不足之處是在分割干涉圖進(jìn)行的過程中，存在將其分塊過小的可能，進(jìn)而導(dǎo)致解纏后圖像出現(xiàn)較為嚴(yán)重的馬賽克現(xiàn)象。在后續(xù)研究中，將探究如何對(duì)自適應(yīng)分塊后的圖像做針對(duì)性的解纏、融合處理。