孫路遙,李向新,2*
(1. 昆明理工大學(xué) 國土資源工程學(xué)院;2. 高原山區(qū)測繪技術(shù)應(yīng)用工程研究中心)
限制山區(qū) PS-InSAR應(yīng)用的一個原因是,其在第一次模型反演時,需要估計該地區(qū)的形變速率及殘余高度。以往多憑經(jīng)驗來大致估計這些參數(shù)的范圍,因而往往存在較大的誤差,并嚴(yán)重影響最終的精度。由于SBAS-InSAR的處理結(jié)果包含形變速率和殘余高度,那么,將其作為 PS-InSAR待估參數(shù)的初值自然是很好的選擇。但SBAS-InSAR的處理結(jié)果較真實(shí)值也存在一定地偏差,尤其當(dāng)研究區(qū)相干條件較差導(dǎo)致結(jié)果不可靠時,這種偏差未必可以忽略,這意味著需要對初值進(jìn)一步調(diào)整。
不同于選取 GCP時有特定的指標(biāo)(殘差絕對值)可以評估其內(nèi)部的離差程度[1],形變速率和殘余高度的估計卻很難通過某個特定的量化指標(biāo)予以評價,僅能從中間結(jié)果(形變速率圖、高程殘差圖)中觀察到些許的變化,而這種變化是更好或是更差,尚未有統(tǒng)一的標(biāo)準(zhǔn)。這一點(diǎn)十分重要,因為只有建立了標(biāo)準(zhǔn),才會有評價的依據(jù),參數(shù)的調(diào)整才可能是一個收斂的過程。
但大部分情況下,估計初值與最優(yōu)取值相差并不算懸殊,這意味著每次調(diào)整參數(shù)的幅度十分有限,當(dāng)參數(shù)變化很小時,中間結(jié)果的變化往往也是細(xì)微的,那么在這種情況下,也就很難確定一個標(biāo)準(zhǔn),作為評判結(jié)果圖像質(zhì)量變好(差)的依據(jù)。
最后,即使可以通過一定的標(biāo)準(zhǔn)判斷圖像質(zhì)量的變化,還需要考慮調(diào)參的策略,即在最優(yōu)值未知的情況下,如何通過最少的實(shí)驗次數(shù),向最優(yōu)值逼近。
新建村位于云南省麗江市玉龍縣鳴音鄉(xiāng)的下胖落村(右岸)與寧蒗縣翠依鄉(xiāng)庫枝村(左岸)交界處的金沙江中游河段上。該地雖降水充沛,但由于蒸發(fā)量大加之貧瘠的土壤條件,使得該地鮮有大型喬木及連片的灌木生長,是典型的熱干河谷,一定程度上減少了因植被茂密造成的影像間的失相干現(xiàn)象。裁剪出的研究區(qū)是一個3′×3′的矩形區(qū)域,其經(jīng)度范圍為 100°25′30″E~100°28′30″E,緯度范圍為27°35′30″N~27°38′30″N。圖 1 是研究區(qū)的影像圖,圖2是2019年4月時的地表植被情況。
圖1 新建地區(qū)影像圖Fig.1 Image map of xinjian
圖2 新建地區(qū)地表植被情況Fig.2 Surface vegetation in xinjian
此次試驗選用的是Sentinel-1A自2018年6月17日至2019年4月13日拍攝的25景升軌SAR影像。文獻(xiàn)[2]的研究表明,當(dāng)衛(wèi)星入射角大于 50度時,相干系數(shù)和后向散射系數(shù)都會大幅降低,造成干涉效果變差,區(qū)域數(shù)據(jù)不可信,由于衛(wèi)星成像幾何及坡度和坡向的原因,選擇升軌時衛(wèi)星視線方向的入射角較降軌更小,有更好的可視性,因此這里選擇的是升軌數(shù)據(jù)。此外,還選用了ALOS World 3D的30 m分辨率的DEM以及精密定軌星歷數(shù)據(jù)。此外,本文使用的是ENVI5.5環(huán)境下的SARScape5.5.1模塊。
嚴(yán)格來講,引言中提出了三個問題:(1)調(diào)參幅度有限導(dǎo)致結(jié)果變化不明顯;(2)即使變化明顯,能歸納出什么樣的標(biāo)準(zhǔn);(3)即使標(biāo)準(zhǔn)明確,怎樣使參數(shù)快速收斂。
前文提到,當(dāng)參數(shù)變化很小時,中間結(jié)果的變化往往也是細(xì)微的。但反過來想,當(dāng)參數(shù)變化很大時,是否就意味著中間結(jié)果的變化會很大,進(jìn)而可以確定一個標(biāo)準(zhǔn),作為評判結(jié)果圖像質(zhì)量變好(差)的依據(jù)。因而下面的實(shí)驗展示的是一個如何從零向最優(yōu)取值逼近的過程,而非從初值(即SBAS-InSAR的結(jié)果)開始估計,就是為了可以觀察到更明顯的中間結(jié)果變化過程,以更好地理解質(zhì)量判定標(biāo)準(zhǔn)的由來。但使用這些準(zhǔn)則時,就需要從中間結(jié)果的細(xì)微變化中仔細(xì)區(qū)分了。
如果最優(yōu)取值的中間結(jié)果與偏差值存在明顯差異,且會隨偏差的增大放大這種差異,此時已隱晦地指出了參數(shù)收斂的必要前提是:在最優(yōu)取值的兩側(cè),成果圖的質(zhì)量是遞減的(不一定是線性)。在此背景下,使參數(shù)快速收斂的問題可近似地理解為:在一正態(tài)分布曲線中,x范圍有界,且對于任意給出的12xx、可比較函數(shù)值的大小,如何最快地向期望值μ逼近。
注意到正態(tài)分布曲線中,期望值μ的左(右)側(cè)是有序遞增(減)的,基于此特點(diǎn),自然地想到了計算科學(xué)中常見的二分法。二分查找是一種非常高效的搜索方法,主要原理是每次搜索可以拋棄一半的值來縮小范圍,其時間復(fù)雜度是()2O log n,一般用于對普通搜索方法的優(yōu)化。盡管二分查找與正態(tài)分布曲線的適用條件尚存一定差距(前者多適用于有序數(shù)組,而后者是連續(xù)曲線;前者適用的數(shù)組是單向有序的,而后者則以期望值μ為界,雙向有序),但稍加調(diào)整,還是可以克服的。
事實(shí)上,本文是先做了大量的實(shí)驗,進(jìn)而歸納出了一些標(biāo)準(zhǔn),但為了把重點(diǎn)落在二分法的闡述上,下面會先給出這些標(biāo)準(zhǔn),以便于對實(shí)驗思路的理解。
PS-InSAR第一次模型反演時,其得到的中間結(jié)果包含了形變速率圖及高程改正圖,前者是對形變速率的初步估計,后者是對輸入 DEM 殘余高度的修正。由于SAR是相干系統(tǒng),斑點(diǎn)噪聲是其固有特性。即使區(qū)域均勻,圖像也會表現(xiàn)出明顯的亮度隨機(jī)變化,因而就很難通過某塊區(qū)域的灰度直方圖判斷其噪聲類型[3]。幸運(yùn)的是,有些噪聲可以通過對圖像的觀察直接或間接地分辨出來,如椒鹽噪聲和均勻分布噪聲。從實(shí)驗結(jié)果來看,在第一次模型反演時處理得越好(這里指參數(shù)越靠近最優(yōu)取值),這些結(jié)果圖會顯現(xiàn)出一些優(yōu)良的特性:
(1)空間連續(xù)性顯著。文獻(xiàn)[1]首次提出了空間連續(xù)性的概念,作者認(rèn)為,無論是滑坡的整體失穩(wěn)亦或是局部失穩(wěn),其在空間維度上都應(yīng)具有一定的整體性,即達(dá)到一定的規(guī)模,在平面內(nèi)連通一定的面積,而非零散的幾個像元。那么類似地,對于形變速率和高程改正值而言,其也應(yīng)有一定的空間連續(xù)性,即,若某片區(qū)域失穩(wěn),其整體的形變速率以及因形變累計造成的高程改正累積在一定空間范圍內(nèi)也應(yīng)觀察到連片現(xiàn)象。從另一個角度看,這種空間連續(xù)性不僅適用于失穩(wěn)區(qū)域,對于穩(wěn)定區(qū)域也應(yīng)是適用的。
(2)分塊性不明顯。當(dāng)研究區(qū)面積大于 52km時,一般的處理策略是對其進(jìn)行分塊處理。假設(shè)某次模型反演沒有任何瑕疵,得到的形變速率圖跟高程改正圖也比較接近真實(shí)的位移,那么塊與塊之間,也應(yīng)是自然流暢地過渡,不應(yīng)存在明顯的像素值的跳躍。如同光學(xué)影像間的拼接一般,要盡量做到?jīng)]有痕跡。
為激發(fā)進(jìn)口潛力,加快跨境電商等新業(yè)態(tài)新模式發(fā)展,2016年5月以來,我國對跨境電商零售進(jìn)口實(shí)行“暫按個人物品監(jiān)管”的過渡期安排,有效促進(jìn)了行業(yè)穩(wěn)定發(fā)展,但也存在各方權(quán)責(zé)不清、政策預(yù)期不穩(wěn)定等問題。6部門聯(lián)合發(fā)布的通知中,一大特點(diǎn)是取消“暫按”,明確按個人自用進(jìn)境物品監(jiān)管,不執(zhí)行首次進(jìn)口許可批件、注冊或備案要求。
(3)椒鹽噪聲很小。胡椒噪聲和鹽噪聲分別對應(yīng)低灰度噪聲和高灰度噪聲。在數(shù)字圖像處理過程中,相較于其他噪聲類型,椒鹽噪聲是最容易直觀辨別的[4]。對應(yīng)于實(shí)際的物理概念,其反應(yīng)的是在連續(xù)的區(qū)域內(nèi)交替出現(xiàn)極大值與極小值,顯然,這是不實(shí)際的。
(4)均勻分布噪聲不明顯。相較于椒鹽噪聲,均勻分布噪聲并不易被直接觀察到。但其仍然有著一些明顯的特點(diǎn),即各個亮度區(qū)間的像元值大體上是均勻分布的。從視覺上來看,圖像內(nèi)邊界不夠銳利,且亮度值從高到低都有存在,不會顯得整體過亮或過暗。反應(yīng)到形變速率圖和高程改正圖上,就是形變零散區(qū)與背景區(qū)差異不明顯。
這里采用了逐次逼近的方法,即根據(jù)中間結(jié)果圖像質(zhì)量的變化情況,不斷地調(diào)整參數(shù),由兩側(cè)向最優(yōu)值逼近。該方法的假設(shè)前提是:在最優(yōu)取值的兩側(cè),成果圖的質(zhì)量是遞減的(不一定是線性)。具體可分為四步:第一步,固定其他參數(shù),對待定參數(shù)取過大值a和過小值b以及二者的某一中間值c分別模型反演;第二步,依據(jù)分塊性、時空連續(xù)性、椒鹽噪聲和均勻分布噪聲等質(zhì)量標(biāo)準(zhǔn)評判三組成果圖的優(yōu)劣次序;第三步,根據(jù)質(zhì)量遞減原則,判斷最優(yōu)值在前半段或后半段,之后將a、c或c、b作為新的過大、過小值,重復(fù)步驟一至三,直至結(jié)果圖質(zhì)量穩(wěn)定;第四步,固定該參數(shù),并按步驟一至三依次找出其他參數(shù)的最優(yōu)取值。
需要說明的是,實(shí)驗并沒有照此標(biāo)準(zhǔn)把所有參數(shù)依次調(diào)試一遍,而是根據(jù)實(shí)際情況,適當(dāng)?shù)販p少了一些可預(yù)估結(jié)果的參數(shù),以減少工作量。
對于城市地區(qū),當(dāng)不確定有明顯沉降時,形變速率和殘余高度一般設(shè)置為±25 mm/year和±70 m,當(dāng)有已知沉降或高層建筑物較多時,可以適當(dāng)增加取值至±50 mm/year和±100 m。而新建地區(qū)地處峽谷深處,鮮有人煙,只有散落的低矮民房,因而將殘余高度暫設(shè)為±20 m。由于兩景影像可探測到的最大形變?yōu)椴ㄩL的一半(28.15 mm),換算成年為867.6 mm/year,故將其設(shè)為過大值。這里將形變速率參數(shù)設(shè)為±400 mm/year、±500 mm/year、±867 mm/year,得到的形變速率圖和高程殘差圖如下。
圖3 a 不同參數(shù)的形變速率圖(±400mm/year、 ±500 mm/year、 ±867 mm/year )Fig.3a Deformation rate graph under different parameters
很明顯,三者的高程殘差圖都有明顯的分塊現(xiàn)象,即都很差,這意味著應(yīng)適當(dāng)修正殘余高度的值。此外,圖3.1a和圖3.2a中都觀察到了明顯的空間連續(xù)性現(xiàn)象,即存在連續(xù)的區(qū)域像元值相等,而圖3.3a中僅有模糊的輪廓,因而可以說明,圖3.3a質(zhì)量較差。由于圖3.1a和圖3.2a比較相近,因而暫時擱置,先調(diào)整殘余高度的取值,這里將形變速率參數(shù)統(tǒng)一設(shè)為±500 mm/year,殘余高度設(shè)為±10 m、±25 m、±50 m、±60 m、±70 m,得到的形變速率圖和高程殘差圖如下。
圖4 a 不同參數(shù)的形變速率圖(a1~a5:±500 mm/year )Fig.4a Deformation rate graph under different parameters
圖4b 不同參數(shù)的高程殘差圖(±10 m、 ±25 m 、 ±50 m 、 ±60 m 、 ±70 m )Fig.4b Elevation residual plots with different parameters
圖4 .1b較圖4.3b,分塊性更為明顯,塊內(nèi)區(qū)分性進(jìn)一步降低,可見將殘余高度調(diào)小是一種錯誤的方向。圖4.3b、圖4.4b、圖4.5b中,高程殘差圖的分塊性依次遞減,細(xì)節(jié)更加突出,可見將該值調(diào)高是正確的方向。較之圖4.4b和圖4.5b,圖4.3b連續(xù)區(qū)域內(nèi)雜點(diǎn)(椒鹽噪聲)略多一點(diǎn)。此外,圖 4.4b背景(非連續(xù)區(qū)域)內(nèi)的差異較其他兩幅圖小一些,但連續(xù)區(qū)與背景區(qū)的差異也表現(xiàn)出同樣的趨勢,因此可以認(rèn)為,圖4.4b的均勻分布噪聲會更大一點(diǎn)。為了進(jìn)一步區(qū)分三者的差異,移動窗口至邊緣區(qū)域,結(jié)果如下。
圖5a 邊緣區(qū)域的形變速率圖(a1~a5:±500 mm/year )Fig.5a Deformation rate graph under different parameters in Marginal area
圖5 b 邊緣區(qū)域不同參數(shù)的高程殘差圖(±50 m、 ± 6 0 m、 ±70 m )Fig.5b Elevation residual plots with different parameters in Marginal area
首先看到圖5.3b中再次觀察到較為明顯的分塊現(xiàn)象,說明高程殘差設(shè)到±70 m已經(jīng)開始遠(yuǎn)離最優(yōu)取值了。同時還注意到,圖5.2a和圖5.3a觀察到了新的細(xì)節(jié),即藍(lán)色的條帶,由于其存在些許蜿蜒且不位于分塊處理的邊界,且故該條帶應(yīng)為河流或存在一定高差的連續(xù)邊坡。因此,殘余高度設(shè)為±60 m時結(jié)果會更好一些。為進(jìn)一步確定殘余高度的最優(yōu)取值并分析其特點(diǎn),下面給出殘余高度設(shè)為±52 m、±55 m、±58 m、±62 m、±65 m得到的形變速率圖和高程殘差圖。
圖6 a 不同參數(shù)的形變速率圖(a1~a5:±500 mm/year )Fig.6a Deformation rate graph under different parameters
首先看到圖6.1a依然沒有觀察到藍(lán)色的條帶,以 及圖6.5a、圖6.5b中不太明顯的分塊邊界,這說明±52 m以及±65 m仍不太合適,同時注意到,±55 m~±62 m之間的差異并不明顯,說明最佳參數(shù)不見得是一個值,而是一個范圍。其次,觀察到±50 m與±52 m,±62 m與±65 m之間,取值很接近但結(jié)果的細(xì)節(jié)差異顯著,說明在最優(yōu)取值的邊緣,模型反演質(zhì)量會極速下滑。此外,后十次測試設(shè)置的形變參數(shù)速率相同,僅不斷調(diào)試殘余高度的取值,但形變速率圖也發(fā)生了一些變化,說明二者是互相影響的。
本文首先提出用SBAS-InSAR的處理結(jié)果作為PS-InSAR的形變參數(shù)估計初值是可行的,之后從初值精度不高出發(fā),引出了基于中間結(jié)果圖像質(zhì)量改進(jìn)參數(shù)精度所面臨的三個問題:結(jié)果變化不明顯、缺乏評價標(biāo)準(zhǔn)以及如何快速收斂。給出的示例,直接或間接地回答了以上問題:增加調(diào)參幅度以放大結(jié)果的差異、從結(jié)果的差異中總結(jié)一些標(biāo)準(zhǔn)、據(jù)此標(biāo)準(zhǔn)利用改進(jìn)的二分法達(dá)到快速收斂。
這里對二分法的改進(jìn)之處是:通過中間結(jié)果圖像質(zhì)量與兩端取值的靠近程度,來判斷最優(yōu)取值在前半段還是后半段,之后對該段再次取中間值以同樣的方法處理,直至圖像質(zhì)量穩(wěn)定。
本實(shí)驗通過增加調(diào)參幅度放大 PS-InSAR模型反演結(jié)果的差異,進(jìn)而總結(jié)了一些質(zhì)量評價標(biāo)準(zhǔn),用以在調(diào)整估計參數(shù)時,判斷調(diào)整結(jié)果的優(yōu)劣。這些標(biāo)準(zhǔn)取材于參數(shù)的實(shí)際物理意義及常見的圖像質(zhì)量指標(biāo),并無過于模糊之處,但不可否認(rèn),該標(biāo)準(zhǔn)仍有很大改進(jìn)及完善的余地。
同時還注意到,本方法效率不高,工作量大,每一幅結(jié)果都要做一次處理,費(fèi)時費(fèi)力。若純粹為了PS-InSAR輸入的估計參數(shù)精確一些,顯然沒有必要。由于SqueeSAR本質(zhì)上仍是一種PS-InSAR,也就意味著,其同樣面臨著參數(shù)估計的問題。不同于滑坡探測,滑坡監(jiān)測以及滑坡預(yù)測對測量精度都有著極高的要求,而 SqueeSAR在滑坡監(jiān)測以及滑坡預(yù)測方面展現(xiàn)出的巨大應(yīng)用潛力[5][6]以及當(dāng)前并不足夠成熟、開放的技術(shù)細(xì)節(jié),都要求我們盡可能地提高每一步流程的精度,因而本文的工作對后人開展 SqueeSAR在滑坡預(yù)測方面的應(yīng)用研究亦可提供一定的參考價值。
最后需要指出的是,本文關(guān)注的重點(diǎn),是因應(yīng)估計參數(shù)設(shè)置有偏,之后根據(jù)怎樣的標(biāo)準(zhǔn)如何一步步修正的問題,與文獻(xiàn)[7]參數(shù)設(shè)置無誤但區(qū)域參考點(diǎn)落到失相干區(qū)域時造成局部圖像中間結(jié)果質(zhì)量變差是有本質(zhì)區(qū)別的。