簡(jiǎn)介本資源是一套面向遙感與InSAR研究者的MATLAB相位解纏實(shí)踐代碼包聚焦干涉SAR數(shù)據(jù)處理中的核心難點(diǎn)——2π周期性相位展開(kāi)問(wèn)題適用于地表形變監(jiān)測(cè)、地質(zhì)災(zāi)害評(píng)估等科研與工程場(chǎng)景適合具備基礎(chǔ)SAR知識(shí)和MATLAB編程能力的研究生、科研人員及工程師。壓縮包共10個(gè)文件41KB含7個(gè)核心.m函數(shù)如QualityGuidedUnwrap2D、BranchCuts、GoldsteinUnwrap2D等、2個(gè)說(shuō)明類(lèi)txt文檔及1個(gè)示例干涉相位數(shù)據(jù).mat文件分別實(shí)現(xiàn)質(zhì)量圖指導(dǎo)法與枝切法兩種主流解纏策略并集成相位殘差檢測(cè)、質(zhì)量圖計(jì)算、洪水填充等關(guān)鍵子模塊。已有2285人學(xué)習(xí)下載代碼結(jié)構(gòu)清晰、注釋完整可直接運(yùn)行驗(yàn)證算法效果支持參數(shù)調(diào)優(yōu)與結(jié)果可視化是理解InSAR相位解纏原理、開(kāi)展算法對(duì)比實(shí)驗(yàn)及實(shí)際數(shù)據(jù)處理的實(shí)用工具集。1. 相位解纏為什么是InSAR處理里繞不過(guò)去的坎干過(guò)InSAR處理的人都有體會(huì)一幅干涉圖生成之后最讓人頭疼的往往不是濾波不是配準(zhǔn)而是相位解纏。原因很簡(jiǎn)單——干涉相位本身是纏繞的它的取值范圍被限制在(-π, π]之間而真實(shí)的地形形變相位其實(shí)是連續(xù)變化的可能遠(yuǎn)超這個(gè)范圍。換句話說(shuō)我們手里拿到的干涉圖是一張“折疊”過(guò)的圖解纏就是把它重新“展開(kāi)”。先想明白一個(gè)問(wèn)題為什么相位一定是纏繞的這和SAR系統(tǒng)的成像機(jī)制有關(guān)。InSAR通過(guò)兩幅SAR影像的干涉獲取相位差這個(gè)相位差是由往返路徑差決定的而路徑差又以波長(zhǎng)為周期。以C波段波長(zhǎng)約5.6 cm為例視線向形變僅2.8 cm就會(huì)引起一個(gè)完整的2π相位周期。如果形變量超過(guò)這個(gè)量級(jí)相位就會(huì)發(fā)生“卷繞”。更麻煩的是地面上的陡峭地形、大氣延遲、噪聲都會(huì)疊加在相位上讓纏繞模式變得極其復(fù)雜。所以解纏不是簡(jiǎn)單地加上2π的整數(shù)倍而是要估計(jì)出每一個(gè)像素對(duì)應(yīng)的整周模糊度integer ambiguity——也就是那個(gè)“2π的倍數(shù)”。這也就解釋了為什么干涉圖質(zhì)量會(huì)直接決定解纏的成敗。如果圖像里噪聲太大或者存在大范圍的去相干區(qū)域殘差點(diǎn)residue point就會(huì)密集分布解纏路徑會(huì)被切斷結(jié)果產(chǎn)生所謂的“跳變”phase jump導(dǎo)致最終形變圖出現(xiàn)明顯的分塊或條紋斷裂。我自己早期的教訓(xùn)是拿到一幅質(zhì)量很差的干涉圖不做任何檢查就直接跑解纏結(jié)果出來(lái)的圖看起來(lái)“很平滑”實(shí)際上是算法把噪聲也一起解纏了形變量完全失真。所以在做任何解纏之前第一步永遠(yuǎn)不是寫(xiě)代碼而是評(píng)估干涉圖質(zhì)量。相干性圖、殘差點(diǎn)密度、條紋清晰度這三個(gè)指標(biāo)能直接告訴你這幅圖值不值得解纏或者應(yīng)該在解纏前做哪些預(yù)處理。2. 在寫(xiě)解纏代碼之前先花10分鐘檢查干涉圖質(zhì)量2.1 相干性圖最直觀的質(zhì)量標(biāo)尺相干性coherence是評(píng)估干涉圖質(zhì)量最核心的指標(biāo)它反映了兩次成像期間地面散射特性的一致性。數(shù)學(xué)上相干性估計(jì)通常用空間窗口內(nèi)的樣本協(xié)方差來(lái)計(jì)算% 假設(shè) intf 是復(fù)數(shù)干涉圖float32尺寸為 [rows, cols] % 窗口大小一般取 5x5 或 9x9視分辨率而定 win 5; kernel ones(win) / (win^2); intensity1 abs(slc1).^2; % SLC1強(qiáng)度 intensity2 abs(slc2).^2; % SLC2強(qiáng)度 cross slc1 .* conj(slc2); % 分別做空間平均 mean_int1 conv2(intensity1, kernel, same); mean_int2 conv2(intensity2, kernel, same); mean_cross conv2(cross, kernel, same); coherence abs(mean_cross) ./ sqrt(mean_int1 .* mean_int2);注意如果手里沒(méi)有原始SLC只有復(fù)數(shù)干涉圖也可以用干涉圖的幅度即相干斑的統(tǒng)計(jì)特性來(lái)間接判斷但不如上述方法準(zhǔn)確。相干性值域在0到1之間經(jīng)驗(yàn)上大于0.5的區(qū)域?qū)儆诟哔|(zhì)量區(qū)解纏基本沒(méi)有問(wèn)題0.3到0.5屬于勉強(qiáng)可用區(qū)需要配合濾波或質(zhì)量引導(dǎo)小于0.3的區(qū)域基本不可信強(qiáng)行解纏只會(huì)帶來(lái)噪聲。這就像看一張照片的清晰度模糊的區(qū)域你強(qiáng)行去識(shí)別上面的文字猜出來(lái)的東西大概率是錯(cuò)的。2.2 殘差點(diǎn)分布解纏的“地雷圖”殘差點(diǎn)是InSAR相位解纏里最核心的概念之一。它的定義是沿著一個(gè)2×2像素的閉環(huán)依次計(jì)算相鄰像素間的纏繞相位差然后求和。如果這個(gè)和不為0說(shuō)明這個(gè)閉環(huán)處的相位場(chǎng)是不守恒的——就像地形圖里出現(xiàn)了“高差對(duì)不攏”的地方。這里給出MATLAB計(jì)算殘差點(diǎn)的核心邏輯function residues calculate_residues(phase) % phase: 纏繞相位矩陣 [rows, cols] [rows, cols] size(phase); residues zeros(rows, cols); % 2x2閉環(huán) for i 1:rows-1 for j 1:cols-1 d1 wrapToPi(phase(i, j1) - phase(i, j)); d2 wrapToPi(phase(i1, j1) - phase(i, j1)); d3 wrapToPi(phase(i1, j) - phase(i1, j1)); d4 wrapToPi(phase(i, j) - phase(i1, j)); sum_d d1 d2 d3 d4; if sum_d 0.5 % 正殘差點(diǎn) residues(i, j) 1; elseif sum_d -0.5 % 負(fù)殘差點(diǎn) residues(i, j) -1; end end end end這里面有個(gè)細(xì)節(jié)容易踩坑wrapToPi函數(shù)在MATLAB里處理的是弧度制輸入輸出范圍是(-π, π]。如果你用的是角度制的相位圖必須先轉(zhuǎn)換而且閉環(huán)求和后的閾值判斷也要相應(yīng)調(diào)整。我見(jiàn)過(guò)不少人在這一步直接把角度制的相位差值套進(jìn)這個(gè)函數(shù)結(jié)果殘差點(diǎn)分布完全不對(duì)后面解纏自然全盤(pán)皆輸。殘差點(diǎn)密度的經(jīng)驗(yàn)判定標(biāo)準(zhǔn)如果一幅干涉圖里殘差點(diǎn)的數(shù)量占總像素?cái)?shù)的比例低于1%屬于優(yōu)秀質(zhì)量如果超過(guò)5%解纏會(huì)比較棘手需要重點(diǎn)考慮濾波或裁剪低相干區(qū)域。殘差點(diǎn)就像地雷枝切法branch cut就是在這些地雷之間搭橋把正負(fù)殘差點(diǎn)連接起來(lái)讓積分路徑避開(kāi)這些不連續(xù)區(qū)域。2.3 濾波解纏前的最后一道防線濾波選什么、怎么選直接影響解纏結(jié)果。常用的有兩種方向空間域?yàn)V波和頻域?yàn)V波?？臻g域里最常見(jiàn)的是Goldstein濾波實(shí)際上是對(duì)干涉條紋頻譜進(jìn)行自適應(yīng)濾波頻域里則是經(jīng)典的boxcar濾波或自適應(yīng)窗口濾波。我個(gè)人的經(jīng)驗(yàn)是不要一上來(lái)就重度濾波。濾波本質(zhì)上是平滑會(huì)損失相位細(xì)節(jié)。對(duì)于質(zhì)量尚可的區(qū)域輕度濾波甚至不濾波反而能保留更多形變細(xì)節(jié)對(duì)于低相干區(qū)域重度濾波也無(wú)法挽回本質(zhì)上去相干的區(qū)域只會(huì)抹平邊界讓解纏結(jié)果看起來(lái)“平滑”但失真。一個(gè)穩(wěn)健的做法是先用相干性圖生成一個(gè)掩膜mask把相干性低于0.25的區(qū)域直接剔除不參與后續(xù)解纏。然后再對(duì)剩余區(qū)域做輕度Goldstein濾波窗口大小取32或64。這比全圖一刀切的濾波方式要靠譜得多因?yàn)樗选耙呀?jīng)死掉”的像素隔離出了處理流程而不是強(qiáng)行去“修復(fù)”它們。mask coherence 0.25; filtered_phase goldstein_filter(phase, coherence, 32); filtered_phase(~mask) 0;3. MATLAB里實(shí)現(xiàn)三種主流解纏算法3.1 枝切法經(jīng)典中的經(jīng)典但別指望它處理爛圖Goldstein枝切法的核心思路是識(shí)別殘差點(diǎn)然后用“樹(shù)枝”連接正負(fù)殘差點(diǎn)使得積分路徑上不會(huì)遇到不成對(duì)的殘差。說(shuō)白了就是先把“地雷”排掉再放心地走路。MATLAB實(shí)現(xiàn)枝切法的完整流程大致包括四步提取殘差點(diǎn)見(jiàn)上文calculate_residues函數(shù)生成枝切線把鄰近的正負(fù)殘差點(diǎn)連接起來(lái)原則是總長(zhǎng)度最短。這一步本質(zhì)上是組合優(yōu)化問(wèn)題常用的是最近鄰匹配或Delaunay三角網(wǎng)搜索設(shè)置障礙枝切線經(jīng)過(guò)的像素在積分時(shí)被跳過(guò)沿路徑積分從參考點(diǎn)開(kāi)始對(duì)不穿過(guò)枝切線的像素逐點(diǎn)解纏% 枝切法主流程偽代碼示意邏輯 % residues: 殘差點(diǎn)圖1為正-1為負(fù)0為正常 branches generate_branch_cuts(residues, max_branch_length); unwrap_phase integrate_along_path(phase, branches, ref_point);這里面最容易出問(wèn)題的是第二步——枝切線的生成策略。如果兩個(gè)殘差點(diǎn)距離過(guò)遠(yuǎn)強(qiáng)行連成一條長(zhǎng)樹(shù)枝反而會(huì)切斷大片有效區(qū)域。一般會(huì)設(shè)置一個(gè)最大枝切長(zhǎng)度閾值比如20個(gè)像素。超過(guò)這個(gè)距離的殘差點(diǎn)寧可留在那里或者直接裁掉也不要連出超長(zhǎng)的樹(shù)枝。枝切法的優(yōu)點(diǎn)是解纏結(jié)果保留了相位的“硬邊界”不會(huì)像最小二乘法那樣把突變區(qū)域抹平。缺點(diǎn)也很明確殘差過(guò)多時(shí)樹(shù)枝會(huì)密集到把有效區(qū)域切割得支離破碎導(dǎo)致大片區(qū)域的解纏值缺失或者出現(xiàn)明顯跳變。所以枝切法更適合高質(zhì)量干涉圖——那種殘差點(diǎn)稀疏、噪聲少的圖。3.2 最小二乘解纏全局優(yōu)化的穩(wěn)健選擇最小二乘法的思路是找一個(gè)“全局最優(yōu)”的解纏相位使得它的梯度相鄰像素差在最小二乘意義下最接近觀測(cè)到的纏繞相位梯度。它不追求每個(gè)像素的精確整周模糊度而是從全局讓誤差最小化。在MATLAB里經(jīng)典的實(shí)現(xiàn)方式是帶權(quán)重的最小二乘解纏通常配合快速離散余弦變換DCT來(lái)求解function unwrapped phase_unwrap_LS(phase, weight) % phase: 纏繞相位 [rows, cols] % weight: 權(quán)重矩陣一般用相干性 [rows, cols] [rows, cols] size(phase); % 計(jì)算梯度x方向和y方向 dx wrapToPi(diff(phase, 1, 2)); dy wrapToPi(diff(phase, 1, 1)); % 構(gòu)建泊松方程右側(cè) rho zeros(rows, cols); rho(:, 2:cols) rho(:, 2:cols) weight(:, 2:cols) .* dx; rho(:, 1:cols-1) rho(:, 1:cols-1) - weight(:, 1:cols-1) .* dx; rho(2:rows, :) rho(2:rows, :) weight(2:rows, :) .* dy; rho(1:rows-1, :) rho(1:rows-1, :) - weight(1:rows-1, :) .* dy; % DCT求解 unwrapped solve_poisson_dct(rho); end這里有幾個(gè)關(guān)鍵點(diǎn)梯度計(jì)算必須用wrapToPi否則差分值仍然會(huì)纏繞求解結(jié)果還是纏繞的這一條最容易犯錯(cuò)權(quán)重矩陣的作用不能省。低相干區(qū)域權(quán)重小解纏時(shí)對(duì)全局優(yōu)化的影響就小可以有效抑制噪聲傳導(dǎo)DCT求解的前提是假設(shè)邊界處梯度為零Neumann邊界條件這對(duì)InSAR數(shù)據(jù)基本是合理的最小二乘法的最大優(yōu)勢(shì)是穩(wěn)健即使殘差很多它也能給出一個(gè)“整體看起來(lái)合理”的結(jié)果。代價(jià)是真實(shí)形變中如果存在斷層或陡峭的形變梯度比如地震同震形變的斷層處最小二乘會(huì)把這種突變“抹平”導(dǎo)致形變梯度被低估。所以在斷層形變研究里我通常更偏向枝切法或者質(zhì)量引導(dǎo)法。3.3 質(zhì)量引導(dǎo)法把好像素先用起來(lái)質(zhì)量引導(dǎo)法Quality-Guided Phase Unwrapping的核心思想非常直觀先從高質(zhì)量區(qū)域高相干性、低殘差密度開(kāi)始解纏然后像水波擴(kuò)散一樣逐步向低質(zhì)量區(qū)域推進(jìn)。這樣能保證誤差盡可能被“關(guān)”在低質(zhì)量區(qū)域不會(huì)大面積擴(kuò)散。實(shí)現(xiàn)質(zhì)量引導(dǎo)法的關(guān)鍵有兩個(gè)質(zhì)量圖的構(gòu)建和排序策略。質(zhì)量圖可以用相干性圖直接充當(dāng)也可以用相位導(dǎo)數(shù)方差phase derivative variance來(lái)構(gòu)建——后者對(duì)條紋密集區(qū)域更敏感。% 相位導(dǎo)數(shù)方差質(zhì)量圖示意 qual zeros(rows, cols); for i 2:rows-1 for j 2:cols-1 % 計(jì)算4鄰域相位導(dǎo)數(shù)的方差 dzx wrapToPi(phase(i, j) - phase(i, j-1)); dzy wrapToPi(phase(i, j) - phase(i-1, j)); % 實(shí)際實(shí)現(xiàn)需要計(jì)算鄰域內(nèi)的統(tǒng)計(jì)量 qual(i, j) sqrt(var([dzx, dzy, ...])); % 值越小質(zhì)量越高 end end排序策略上最簡(jiǎn)單的方法是堆垛法flood fill priority queue初始選取一個(gè)質(zhì)量最高的種子點(diǎn)將它加入隊(duì)列每次從隊(duì)列中取出質(zhì)量最高的像素解纏它并把它的四個(gè)鄰域如果還沒(méi)解纏加入隊(duì)列。MATLAB里可以用containers.Map配合排序或者直接用sortrows維護(hù)一個(gè)按質(zhì)量值排序的列表數(shù)據(jù)量不大時(shí)效率足夠。質(zhì)量引導(dǎo)法的優(yōu)勢(shì)在于它能充分利用干涉圖里“還不錯(cuò)的”區(qū)域即使整體質(zhì)量一般也能得到連貫的解纏結(jié)果。缺點(diǎn)是對(duì)孤立低質(zhì)量區(qū)域的解纏能力弱如果低相干區(qū)域被高質(zhì)量區(qū)域包圍解纏值會(huì)被“鎖死”可能出現(xiàn)孤島狀錯(cuò)誤。3.4 三種算法的選型建議算法適用場(chǎng)景優(yōu)點(diǎn)缺點(diǎn)MATLAB實(shí)現(xiàn)復(fù)雜度枝切法高質(zhì)量干涉圖、斷層形變保留突變邊界低質(zhì)量圖效果差中等最小二乘法大面積形變、噪聲較多穩(wěn)健、全局最優(yōu)平滑掉突變較低DCT求解質(zhì)量引導(dǎo)法質(zhì)量參差不齊的干涉圖自適應(yīng)、靈活孤立低質(zhì)量區(qū)域易出錯(cuò)較高實(shí)操建議實(shí)際項(xiàng)目里我通常先用最小二乘法跑一遍全圖得到一個(gè)參考解纏結(jié)果再對(duì)重點(diǎn)關(guān)注區(qū)域比如形變梯度大的斷層附近用枝切法或質(zhì)量引導(dǎo)法細(xì)化。兩種結(jié)果對(duì)比可以快速定位潛在的解纏錯(cuò)誤區(qū)域。4. 完整實(shí)操?gòu)囊环缮鎴D到解纏結(jié)果4.1 數(shù)據(jù)準(zhǔn)備和參數(shù)設(shè)定假設(shè)我們手頭有一幅由GAMMA或ISCE生成的復(fù)數(shù)干涉圖intf.float數(shù)據(jù)格式為float32復(fù)數(shù)尺寸為500×500以及對(duì)應(yīng)的SLC1和SLC2。以下幾行代碼是處理流程的基礎(chǔ)% 讀取復(fù)數(shù)干涉圖 fid fopen(intf.float, rb); intf fread(fid, [500, 500], float32); fclose(fid); phase angle(intf); % 纏繞相位值域 [-pi, pi] amp abs(intf); % 幅度信息這里一個(gè)很常見(jiàn)的坑是數(shù)據(jù)字節(jié)序問(wèn)題。GAMMA默認(rèn)輸出的是小端序little-endian但不同版本可能有差異。如果讀出來(lái)的數(shù)據(jù)明顯是“花屏”狀態(tài)先檢查fread是否需要加參數(shù)l或b。另外注意矩陣讀入后是否需要轉(zhuǎn)置——GAMMA輸出是按行優(yōu)先存儲(chǔ)的但MATLAB默認(rèn)按列優(yōu)先讀取所以讀出來(lái)后通常要.T轉(zhuǎn)置一下。4.2 預(yù)處理去平地效應(yīng)在解纏之前如果干涉圖還包含平地相位即由參考橢球面引起的系統(tǒng)性相位變化需要先去掉。常見(jiàn)做法是用軌道信息和成像幾何計(jì)算平地相位并減去或者在頻域里把主頻峰移到中心。% 頻域去平地將干涉圖變換到頻域把零頻移到幅度譜峰值位置 F fft2(intf); [rows, cols] size(F); % 找到幅度譜峰值的位置避開(kāi)零頻附近 shift_x ...; % 通過(guò)尋找峰值計(jì)算 shift_y ...; % 直接在頻域移動(dòng)或者使用相位斜坡擬合均可去平地這一步很多人會(huì)忽略或做錯(cuò)其實(shí)它直接影響后續(xù)的條紋頻率和解纏效果。如果平地沒(méi)去干凈干涉圖里會(huì)出現(xiàn)大量的平行條紋它們的密度很高容易造成殘差點(diǎn)密集分布。4.3 解纏執(zhí)行我默認(rèn)采用質(zhì)量引導(dǎo)法作為主流程因?yàn)樗骖櫫朔€(wěn)健性和邊界保留能力% 1. 構(gòu)建質(zhì)量圖用相干性 coherence estimate_coherence(slc1, slc2, 5); % 5x5窗口 % 2. 低相干掩膜 mask coherence 0.3; % 3. 質(zhì)量引導(dǎo)解纏 unwrapped_phase quality_guided_unwrap(phase, coherence, mask); % 4. 去除參考點(diǎn)通常選一個(gè)高相干、遠(yuǎn)離形變區(qū)的點(diǎn) ref_idx ...; % 參考點(diǎn)像素坐標(biāo) unwrapped_phase unwrapped_phase - unwrapped_phase(ref_idx);這里有一個(gè)容易被忽略的細(xì)節(jié)參考點(diǎn)的選擇會(huì)直接影響最終形變的絕對(duì)量級(jí)。所有解纏結(jié)果都是相對(duì)于參考點(diǎn)的相對(duì)值參考點(diǎn)和形變區(qū)如果在同一幅圖內(nèi)其自身可能也在形變就會(huì)導(dǎo)致全圖的形變被“抬升”或“下沉”。所以參考點(diǎn)一定要選在形變區(qū)之外最好結(jié)合實(shí)際地面情況如基巖、穩(wěn)定建筑區(qū)來(lái)定。4.4 結(jié)果輸出和可視化解纏完成后輸出是最容易忽略卻也很重要的環(huán)節(jié)。因?yàn)楹罄m(xù)往往要用GIS或者其他軟件做進(jìn)一步的形變分析數(shù)據(jù)格式要提前想好。% 轉(zhuǎn)換為形變值以C波段為例單位米 lambda 0.056; % 波長(zhǎng) los_displacement unwrapped_phase * lambda / (4 * pi); % 保存為GeoTIFF需要映射信息 geotiffwrite(los_displacement.tif, los_displacement, R, CoordRefSysCode, 32650);注意這里視線向形變的符號(hào)約定要小心。不同軟件GAMMA、ISCE、SNAP對(duì)形變方向的正負(fù)號(hào)定義不完全一致導(dǎo)出前一定確認(rèn)清楚否則做出來(lái)的形變圖在符號(hào)上是反的明明沉降會(huì)被畫(huà)成抬升。5. 解纏過(guò)程中最常踩的坑和排查方法5.1 解纏結(jié)果出現(xiàn)“跳變”或“條紋斷裂”這個(gè)問(wèn)題的典型表現(xiàn)是解纏后的相位圖在某一區(qū)域出現(xiàn)明顯的高低值突變甚至相差多個(gè)2π周期。原因通常有三個(gè)一是殘差點(diǎn)密度過(guò)高枝切線或質(zhì)量引導(dǎo)路徑繞不過(guò)去二是低相干區(qū)域形成了“通道”噪聲從通道擴(kuò)散到了有效區(qū)域三是濾波窗口不合適把真實(shí)的相位突變也平滑掉了。排查方法很簡(jiǎn)單先把掩膜mask疊加在解纏結(jié)果上看跳變位置是否和低相干區(qū)域?qū)?yīng)。如果是說(shuō)明是掩膜閾值設(shè)置太低把噪聲區(qū)納入了解纏范圍。如果跳變出現(xiàn)在高相干區(qū)域那多半是解纏算法本身的路徑選擇出了問(wèn)題可以嘗試改用質(zhì)量引導(dǎo)法或者調(diào)節(jié)枝切長(zhǎng)度閾值。5.2 解纏結(jié)果非常平滑但總覺(jué)得形變梯度被削弱了這種情況多半出在最小二乘法上。最小二乘解的固有特性就是“能量最小化”它會(huì)盡可能地把相鄰像素的差異拉小所以真實(shí)形變中的陡峭梯度比如斷層會(huì)被弱化。如果研究目標(biāo)是地震形變或滑坡邊界建議改用枝切法或混合方法先在低相干區(qū)域用最小二乘法給一個(gè)初始估計(jì)再在高相干區(qū)域用枝切法修正。5.3 解纏速度慢到無(wú)法忍受對(duì)于大范圍干涉圖比如10萬(wàn)×10萬(wàn)像素即使是MATLAB也需要考慮效率問(wèn)題。優(yōu)化思路有兩個(gè)方向降采樣再解纏先把干涉圖降采樣到1/4或1/16大小解纏得到粗結(jié)果再用粗結(jié)果作為初值在原分辨率下做局部修正。這個(gè)思路類(lèi)似金字塔策略速度快且穩(wěn)定。分塊解纏把干涉圖切成有重疊的小塊分別解纏后拼接。注意要保證塊與塊之間有足夠的重疊區(qū)域推薦不小于256像素并且對(duì)齊時(shí)利用重疊區(qū)域的平均相位差來(lái)消除塊間偏移。% 分塊解纏的邊界對(duì)齊關(guān)鍵步驟 % blk1, blk2: 兩塊解纏結(jié)果overlap_region為重疊區(qū) offset median(unwrapped_blk1(overlap_region) - unwrapped_blk2(overlap_region)); unwrapped_blk2 unwrapped_blk2 offset;5.4 相干性不低、卻解纏錯(cuò)誤的情況這種情況最常見(jiàn)的原因是相位混疊——干涉條紋太密超出了采樣率能夠承載的范圍。當(dāng)天線的空間基線過(guò)長(zhǎng)、地形起伏過(guò)大時(shí)局部干涉條紋頻率可能接近甚至超過(guò)奈奎斯特頻率此時(shí)相位在相鄰像素間本身就跳變了超過(guò)π任何解纏算法都無(wú)法恢復(fù)。如果遇到這種情況處理方向不在解纏算法本身而在干涉圖生成之前縮短空間基線選擇時(shí)間基線更近的影像對(duì)、做外部DEM輔助去除地形相位、或者使用多孔徑InSARMAI等替代技術(shù)。5.5 常見(jiàn)問(wèn)題速查表問(wèn)題表現(xiàn)可能原因應(yīng)對(duì)策略解纏結(jié)果有大面積亂碼掩膜未用低相干區(qū)參與解纏檢查掩膜閾值低于0.3區(qū)域剔除跳變沿特定方向分布?xì)埐铧c(diǎn)成串分布枝切線過(guò)長(zhǎng)減小最大枝切長(zhǎng)度改用質(zhì)量引導(dǎo)法形變梯度明顯偏小最小二乘平滑效應(yīng)改用枝切法或混合解纏策略解纏值出現(xiàn)周期性的“條帶”平地效應(yīng)未去除干凈檢查去平地流程頻域?yàn)V波重新處理參考點(diǎn)區(qū)域形變值不為0參考點(diǎn)自身位于形變區(qū)重新選參考點(diǎn)置于穩(wěn)定區(qū)域運(yùn)行內(nèi)存溢出或速度極慢數(shù)據(jù)量過(guò)大或未降采樣分塊處理或金字塔策略降采樣6. 解纏之外的幾個(gè)延伸方向解纏本身只是InSAR形變測(cè)量鏈條中的一環(huán)但解纏質(zhì)量的好壞直接決定了后續(xù)所有產(chǎn)品的可靠性。解纏結(jié)果如果出了問(wèn)題后面無(wú)論是做形變速率估計(jì)、時(shí)間序列分析還是地球物理反演都會(huì)帶著這個(gè)誤差往下走。我個(gè)人的建議是在項(xiàng)目流程里把解纏質(zhì)量評(píng)估作為一個(gè)獨(dú)立的檢查關(guān)卡。每次解纏完成后不要急著出圖先做三件事一是統(tǒng)計(jì)殘差點(diǎn)在解纏后是否大幅減少二是把解纏結(jié)果和原始干涉條紋疊加顯示目視檢查條紋是否連貫三是如果有外部數(shù)據(jù)如GPS、水準(zhǔn)測(cè)量用獨(dú)立觀測(cè)值驗(yàn)證一個(gè)或幾個(gè)點(diǎn)的形變值。這三步能在早期攔截大部分解纏錯(cuò)誤。另外相位解纏并不是一個(gè)“一勞永逸”的問(wèn)題。對(duì)于不同波段L、C、X、不同地形條件、不同地表覆蓋類(lèi)型最合適的解纏策略可能完全不同。L波段雷達(dá)波長(zhǎng)長(zhǎng)形變相位梯度容易滿足采樣條件解纏相對(duì)容易X波段波長(zhǎng)短對(duì)形變極其敏感但也更容易出現(xiàn)相位混疊。所以做項(xiàng)目的時(shí)候提前根據(jù)波段和區(qū)域特征選擇合適的解纏算法比盲目追求“最強(qiáng)算法”務(wù)實(shí)得多。關(guān)于MATLAB實(shí)現(xiàn)本身我還想多說(shuō)一句?，F(xiàn)在有不少開(kāi)源的解纏工具包比如SNAPHU、snaphu_mex質(zhì)量和效率都很高。如果只是做常規(guī)解纏直接調(diào)用這些工具包完全夠用。自己寫(xiě)MATLAB實(shí)現(xiàn)最大的價(jià)值在于你親手把每一步算了一遍你會(huì)真正理解殘差點(diǎn)、質(zhì)量圖、枝切線這些概念是怎么來(lái)的踩過(guò)坑之后你才不會(huì)把解纏當(dāng)成一個(gè)“黑盒”隨便調(diào)參數(shù)。我自己的習(xí)慣是寫(xiě)完一套解纏代碼之后一定會(huì)用一幅已知的模擬干涉圖正演一個(gè)已知形變場(chǎng)加上纏繞和噪聲來(lái)做驗(yàn)證。如果解纏結(jié)果能精確恢復(fù)出原始形變場(chǎng)說(shuō)明代碼邏輯沒(méi)問(wèn)題才能在真實(shí)數(shù)據(jù)上放心用。這個(gè)驗(yàn)證步驟建議所有剛接觸解纏的人都做一遍。最后分享一點(diǎn)經(jīng)驗(yàn)之談。做InSAR處理尤其是解纏這個(gè)環(huán)節(jié)心態(tài)上要有“誤差終究無(wú)法完全消除只能控制其傳播”的意識(shí)。你不可能讓每一幅干涉圖都解纏得完美無(wú)缺但你可以通過(guò)預(yù)處理、掩膜、算法選型和結(jié)果檢查把解纏誤差控制在一個(gè)可接受的范圍內(nèi)。做好這一步后面的形變分析才會(huì)更可靠。本文還有配套的精品資源點(diǎn)擊獲取