基于二分法的山區PS-InSAR形變參數估計

孫路遙，李向新,2*

（1. 昆明理工大學 國土資源工程學院；2. 高原山區測繪技術應用工程研究中心）

0 引言

限制山區 PS-InSAR應用的一個原因是，其在第一次模型反演時，需要估計該地區的形變速率及殘余高度。以往多憑經驗來大致估計這些參數的范圍，因而往往存在較大的誤差，并嚴重影響最終的精度。由于SBAS-InSAR的處理結果包含形變速率和殘余高度，那么，將其作為 PS-InSAR待估參數的初值自然是很好的選擇。但SBAS-InSAR的處理結果較真實值也存在一定地偏差，尤其當研究區相干條件較差導致結果不可靠時，這種偏差未必可以忽略，這意味著需要對初值進一步調整。

不同于選取 GCP時有特定的指標（殘差絕對值）可以評估其內部的離差程度[1]，形變速率和殘余高度的估計卻很難通過某個特定的量化指標予以評價，僅能從中間結果（形變速率圖、高程殘差圖）中觀察到些許的變化，而這種變化是更好或是更差，尚未有統一的標準。這一點十分重要，因為只有建立了標準，才會有評價的依據，參數的調整才可能是一個收斂的過程。

但大部分情況下，估計初值與最優取值相差并不算懸殊，這意味著每次調整參數的幅度十分有限，當參數變化很小時，中間結果的變化往往也是細微的，那么在這種情況下，也就很難確定一個標準，作為評判結果圖像質量變好（差）的依據。

最后，即使可以通過一定的標準判斷圖像質量的變化，還需要考慮調參的策略，即在最優值未知的情況下，如何通過最少的實驗次數，向最優值逼近。

1 研究區概況

新建村位于云南省麗江市玉龍縣鳴音鄉的下胖落村（右岸）與寧蒗縣翠依鄉庫枝村（左岸）交界處的金沙江中游河段上。該地雖降水充沛，但由于蒸發量大加之貧瘠的土壤條件，使得該地鮮有大型喬木及連片的灌木生長，是典型的熱干河谷，一定程度上減少了因植被茂密造成的影像間的失相干現象。裁剪出的研究區是一個3′×3′的矩形區域，其經度范圍為 100°25′30″E～100°28′30″E，緯度范圍為27°35′30″N～27°38′30″N。圖 1 是研究區的影像圖，圖2是2019年4月時的地表植被情況。

圖1 新建地區影像圖Fig.1 Image map of xinjian

圖2 新建地區地表植被情況Fig.2 Surface vegetation in xinjian

2 數據源及實驗平臺

此次試驗選用的是Sentinel-1A自2018年6月17日至2019年4月13日拍攝的25景升軌SAR影像。文獻[2]的研究表明，當衛星入射角大于 50度時，相干系數和后向散射系數都會大幅降低，造成干涉效果變差，區域數據不可信，由于衛星成像幾何及坡度和坡向的原因，選擇升軌時衛星視線方向的入射角較降軌更小，有更好的可視性，因此這里選擇的是升軌數據。此外，還選用了ALOS World 3D的30 m分辨率的DEM以及精密定軌星歷數據。此外，本文使用的是ENVI5.5環境下的SARScape5.5.1模塊。

3 實驗

3.1 問題分析

嚴格來講，引言中提出了三個問題：（1）調參幅度有限導致結果變化不明顯；（2）即使變化明顯，能歸納出什么樣的標準；（3）即使標準明確，怎樣使參數快速收斂。

前文提到，當參數變化很小時，中間結果的變化往往也是細微的。但反過來想，當參數變化很大時，是否就意味著中間結果的變化會很大，進而可以確定一個標準，作為評判結果圖像質量變好（差）的依據。因而下面的實驗展示的是一個如何從零向最優取值逼近的過程，而非從初值（即SBAS-InSAR的結果）開始估計，就是為了可以觀察到更明顯的中間結果變化過程，以更好地理解質量判定標準的由來。但使用這些準則時，就需要從中間結果的細微變化中仔細區分了。

如果最優取值的中間結果與偏差值存在明顯差異，且會隨偏差的增大放大這種差異，此時已隱晦地指出了參數收斂的必要前提是：在最優取值的兩側，成果圖的質量是遞減的（不一定是線性）。在此背景下，使參數快速收斂的問題可近似地理解為：在一正態分布曲線中，x范圍有界，且對于任意給出的12xx、可比較函數值的大小，如何最快地向期望值μ逼近。

注意到正態分布曲線中，期望值μ的左（右）側是有序遞增（減）的，基于此特點，自然地想到了計算科學中常見的二分法。二分查找是一種非常高效的搜索方法，主要原理是每次搜索可以拋棄一半的值來縮小范圍，其時間復雜度是（）2O log n，一般用于對普通搜索方法的優化。盡管二分查找與正態分布曲線的適用條件尚存一定差距（前者多適用于有序數組，而后者是連續曲線；前者適用的數組是單向有序的，而后者則以期望值μ為界，雙向有序），但稍加調整，還是可以克服的。

事實上，本文是先做了大量的實驗，進而歸納出了一些標準，但為了把重點落在二分法的闡述上，下面會先給出這些標準，以便于對實驗思路的理解。

3.2 中間結果質量判定標準

PS-InSAR第一次模型反演時，其得到的中間結果包含了形變速率圖及高程改正圖，前者是對形變速率的初步估計，后者是對輸入 DEM 殘余高度的修正。由于SAR是相干系統，斑點噪聲是其固有特性。即使區域均勻，圖像也會表現出明顯的亮度隨機變化，因而就很難通過某塊區域的灰度直方圖判斷其噪聲類型[3]。幸運的是，有些噪聲可以通過對圖像的觀察直接或間接地分辨出來，如椒鹽噪聲和均勻分布噪聲。從實驗結果來看，在第一次模型反演時處理得越好（這里指參數越靠近最優取值），這些結果圖會顯現出一些優良的特性：

（1）空間連續性顯著。文獻[1]首次提出了空間連續性的概念，作者認為，無論是滑坡的整體失穩亦或是局部失穩，其在空間維度上都應具有一定的整體性，即達到一定的規模，在平面內連通一定的面積，而非零散的幾個像元。那么類似地，對于形變速率和高程改正值而言，其也應有一定的空間連續性，即，若某片區域失穩，其整體的形變速率以及因形變累計造成的高程改正累積在一定空間范圍內也應觀察到連片現象。從另一個角度看，這種空間連續性不僅適用于失穩區域，對于穩定區域也應是適用的。

（2）分塊性不明顯。當研究區面積大于 52km時，一般的處理策略是對其進行分塊處理。假設某次模型反演沒有任何瑕疵，得到的形變速率圖跟高程改正圖也比較接近真實的位移，那么塊與塊之間，也應是自然流暢地過渡，不應存在明顯的像素值的跳躍。如同光學影像間的拼接一般，要盡量做到沒有痕跡。

為激發進口潛力，加快跨境電商等新業態新模式發展，2016年5月以來，我國對跨境電商零售進口實行“暫按個人物品監管”的過渡期安排，有效促進了行業穩定發展，但也存在各方權責不清、政策預期不穩定等問題。6部門聯合發布的通知中，一大特點是取消“暫按”，明確按個人自用進境物品監管，不執行首次進口許可批件、注冊或備案要求。

（3）椒鹽噪聲很小。胡椒噪聲和鹽噪聲分別對應低灰度噪聲和高灰度噪聲。在數字圖像處理過程中，相較于其他噪聲類型，椒鹽噪聲是最容易直觀辨別的[4]。對應于實際的物理概念，其反應的是在連續的區域內交替出現極大值與極小值，顯然，這是不實際的。

（4）均勻分布噪聲不明顯。相較于椒鹽噪聲，均勻分布噪聲并不易被直接觀察到。但其仍然有著一些明顯的特點，即各個亮度區間的像元值大體上是均勻分布的。從視覺上來看，圖像內邊界不夠銳利，且亮度值從高到低都有存在，不會顯得整體過亮或過暗。反應到形變速率圖和高程改正圖上，就是形變零散區與背景區差異不明顯。

3.3 實驗方法

這里采用了逐次逼近的方法，即根據中間結果圖像質量的變化情況，不斷地調整參數，由兩側向最優值逼近。該方法的假設前提是：在最優取值的兩側，成果圖的質量是遞減的（不一定是線性）。具體可分為四步：第一步，固定其他參數，對待定參數取過大值a和過小值b以及二者的某一中間值c分別模型反演；第二步，依據分塊性、時空連續性、椒鹽噪聲和均勻分布噪聲等質量標準評判三組成果圖的優劣次序；第三步,根據質量遞減原則，判斷最優值在前半段或后半段，之后將a、c或c、b作為新的過大、過小值，重復步驟一至三，直至結果圖質量穩定；第四步，固定該參數，并按步驟一至三依次找出其他參數的最優取值。

需要說明的是，實驗并沒有照此標準把所有參數依次調試一遍，而是根據實際情況，適當地減少了一些可預估結果的參數，以減少工作量。

3.4 實驗過程

對于城市地區，當不確定有明顯沉降時，形變速率和殘余高度一般設置為±25 mm/year和±70 m，當有已知沉降或高層建筑物較多時，可以適當增加取值至±50 mm/year和±100 m。而新建地區地處峽谷深處，鮮有人煙，只有散落的低矮民房，因而將殘余高度暫設為±20 m。由于兩景影像可探測到的最大形變為波長的一半（28.15 mm），換算成年為867.6 mm/year，故將其設為過大值。這里將形變速率參數設為±400 mm/year、±500 mm/year、±867 mm/year，得到的形變速率圖和高程殘差圖如下。

圖3 a 不同參數的形變速率圖(±400mm/year、 ±500 mm/year、 ±867 mm/year )Fig.3a Deformation rate graph under different parameters

很明顯，三者的高程殘差圖都有明顯的分塊現象，即都很差，這意味著應適當修正殘余高度的值。此外，圖3.1a和圖3.2a中都觀察到了明顯的空間連續性現象，即存在連續的區域像元值相等，而圖3.3a中僅有模糊的輪廓，因而可以說明，圖3.3a質量較差。由于圖3.1a和圖3.2a比較相近，因而暫時擱置，先調整殘余高度的取值，這里將形變速率參數統一設為±500 mm/year，殘余高度設為±10 m、±25 m、±50 m、±60 m、±70 m，得到的形變速率圖和高程殘差圖如下。

圖4 a 不同參數的形變速率圖(a1～a5：±500 mm/year )Fig.4a Deformation rate graph under different parameters

圖4b 不同參數的高程殘差圖(±10 m、 ±25 m 、 ±50 m 、 ±60 m 、 ±70 m )Fig.4b Elevation residual plots with different parameters

圖4 .1b較圖4.3b，分塊性更為明顯，塊內區分性進一步降低，可見將殘余高度調小是一種錯誤的方向。圖4.3b、圖4.4b、圖4.5b中，高程殘差圖的分塊性依次遞減，細節更加突出，可見將該值調高是正確的方向。較之圖4.4b和圖4.5b，圖4.3b連續區域內雜點（椒鹽噪聲）略多一點。此外，圖 4.4b背景（非連續區域）內的差異較其他兩幅圖小一些，但連續區與背景區的差異也表現出同樣的趨勢，因此可以認為，圖4.4b的均勻分布噪聲會更大一點。為了進一步區分三者的差異，移動窗口至邊緣區域，結果如下。

圖5a 邊緣區域的形變速率圖(a1～a5：±500 mm/year )Fig.5a Deformation rate graph under different parameters in Marginal area

圖5 b 邊緣區域不同參數的高程殘差圖(±50 m、 ± 6 0 m、 ±70 m )Fig.5b Elevation residual plots with different parameters in Marginal area

首先看到圖5.3b中再次觀察到較為明顯的分塊現象，說明高程殘差設到±70 m已經開始遠離最優取值了。同時還注意到，圖5.2a和圖5.3a觀察到了新的細節，即藍色的條帶，由于其存在些許蜿蜒且不位于分塊處理的邊界，且故該條帶應為河流或存在一定高差的連續邊坡。因此，殘余高度設為±60 m時結果會更好一些。為進一步確定殘余高度的最優取值并分析其特點，下面給出殘余高度設為±52 m、±55 m、±58 m、±62 m、±65 m得到的形變速率圖和高程殘差圖。

圖6 a 不同參數的形變速率圖(a1～a5：±500 mm/year )Fig.6a Deformation rate graph under different parameters

首先看到圖6.1a依然沒有觀察到藍色的條帶，以 及圖6.5a、圖6.5b中不太明顯的分塊邊界，這說明±52 m以及±65 m仍不太合適，同時注意到，±55 m～±62 m之間的差異并不明顯，說明最佳參數不見得是一個值，而是一個范圍。其次，觀察到±50 m與±52 m，±62 m與±65 m之間，取值很接近但結果的細節差異顯著，說明在最優取值的邊緣，模型反演質量會極速下滑。此外，后十次測試設置的形變參數速率相同，僅不斷調試殘余高度的取值，但形變速率圖也發生了一些變化，說明二者是互相影響的。

4 總結

本文首先提出用SBAS-InSAR的處理結果作為PS-InSAR的形變參數估計初值是可行的，之后從初值精度不高出發，引出了基于中間結果圖像質量改進參數精度所面臨的三個問題：結果變化不明顯、缺乏評價標準以及如何快速收斂。給出的示例，直接或間接地回答了以上問題：增加調參幅度以放大結果的差異、從結果的差異中總結一些標準、據此標準利用改進的二分法達到快速收斂。

這里對二分法的改進之處是：通過中間結果圖像質量與兩端取值的靠近程度，來判斷最優取值在前半段還是后半段，之后對該段再次取中間值以同樣的方法處理，直至圖像質量穩定。

5 結論與展望

本實驗通過增加調參幅度放大 PS-InSAR模型反演結果的差異，進而總結了一些質量評價標準，用以在調整估計參數時，判斷調整結果的優劣。這些標準取材于參數的實際物理意義及常見的圖像質量指標，并無過于模糊之處，但不可否認，該標準仍有很大改進及完善的余地。

同時還注意到，本方法效率不高，工作量大，每一幅結果都要做一次處理，費時費力。若純粹為了PS-InSAR輸入的估計參數精確一些，顯然沒有必要。由于SqueeSAR本質上仍是一種PS-InSAR，也就意味著，其同樣面臨著參數估計的問題。不同于滑坡探測，滑坡監測以及滑坡預測對測量精度都有著極高的要求，而 SqueeSAR在滑坡監測以及滑坡預測方面展現出的巨大應用潛力[5][6]以及當前并不足夠成熟、開放的技術細節，都要求我們盡可能地提高每一步流程的精度，因而本文的工作對后人開展 SqueeSAR在滑坡預測方面的應用研究亦可提供一定的參考價值。

最后需要指出的是，本文關注的重點，是因應估計參數設置有偏，之后根據怎樣的標準如何一步步修正的問題，與文獻[7]參數設置無誤但區域參考點落到失相干區域時造成局部圖像中間結果質量變差是有本質區別的。