基于模擬退火和區域劃分的海表流場反演算法研究
2016-07-20 11:28:56張夏榮黃云仙嚴衛趙現斌
海洋學報 2016年7期
張夏榮，黃云仙*，嚴衛，趙現斌
(1.解放軍理工大學 氣象海洋學院，江蘇 南京 211101)
?
基于模擬退火和區域劃分的海表流場反演算法研究
張夏榮1，黃云仙1*，嚴衛1，趙現斌1
(1.解放軍理工大學 氣象海洋學院，江蘇 南京 211101)
摘要：經典的海洋表面流場迭代反演算法采用固定的校正系數對探測區域進行整體分析和計算，因此耗時較長且反演精度也受到了一定程度的限制。本文提出利用模擬退火算法對反演過程中的校正系數進行優化，使其能根據每一次反演的結果進行自適應的調整和改變，從而減少迭代次數；同時根據流場分布的特點，提出劃分探測區域的反演方法，最終，在經典算法的基礎上給出了改進的海表流場反演算法。仿真反演實驗結果表明，改進后的反演算法，其時效性和精度都有了明顯提高。
關鍵詞：順軌干涉SAR；海洋表面流場；M4S模型；模擬退火算法
1引言
海洋表面流場是一種重要的海洋動力學現象，無論是在輸送海洋熱量、鹽度等[1]方面，還是對研究海洋、陸地和大氣之間的相互作用[2]以及海洋災害預警方面都具有十分重大的科學意義[3]。
目前針對海表流場的測量方法主要有兩大類，分別為傳統手段測量和遙感測量。其中傳統的海流測量包括浮標和電磁式海流計等[4]，主要存在采樣點少、成本較高的問題；而通過遙感手段雖然可以獲取空間覆蓋范圍較廣的海表流場數據，但無論是星載高度計還是合成孔徑雷達(Synthetic Aperture Radar，SAR)都有各自的局限性[5—6]。其中星載高度計由于受星下點觀測幾何的限制，其獲得的流場距離向分辨率較為粗糙，約為幾十到幾百千米量級；而合成孔徑雷達的強度圖像并不能直接反映流場速度的變化，且其在成像過程中存在多種非線性調制機制[7]，這也增大了流場反演的難度。
順軌干涉SAR(Along-Track Interferometric SAR，ATI-SAR)是1987年由美國Goldstein和Zebker首先提出來的一種技術，其探測的干涉相位與表面流場的徑向速度分量成正比，據此可以反演高分辨率的海表流場信息[8—9]?；贛4S微波成像仿真模型，Romeiser等[10—12]提出了順軌干涉SAR流場迭代反演算法，于祥禎等[13]給出了詳細的海表流場迭代校正方案，并利用美國噴氣推進實驗室JPL(Jet Propulsion Laboratory)的機載SAR數據開展了反演實驗，驗證了反演方法的有效性。但是，該反演算法對于校正系數的選取未作詳細的說明，而合理選擇校正系數可以減少迭代次數，縮短M4S模型仿真計算的時長，從而提高反演速度。另外，在計算過程中，當目標區域的海表流場變化較大時，利用該算法反演的流場與真實流場存在較大的誤差。本文擬在該經典的反演算法基礎上，對校正算法進行優化，從而有效提高全局流場的反演精度并縮短整個反演過程的計算時間。
2經典的海表流場反演算法
2.1M4S微波成像仿真模型
M4S模型是一個可以仿真計算海表流場、海面風場以及海浪等海洋活動的工具包[14]，它由漢堡大學的Roland Romeiser開發并首次應用于海表流場的仿真和反演。
M4S模型主要包括M4Sw320、M4Sr320和M4Sq320等計算模塊。其中，M4Sw320為海浪譜計算模塊，用于計算給定海面風場、海表流場條件下海面微尺度波的波高譜；M4Sr320為雷達成像計算模塊，利用生成的海浪譜，以及設置的雷達參數和平臺參數等計算多普勒譜、SAR強度和干涉相位圖像；M4Sq320為單點計算模塊，計算單點雷達散射信號或者多普勒譜。
采用M4S模型開展海表流場反演研究時，給定海面風場、海表流場、雷達參數和平臺參數，可以計算順軌干涉SAR的強度和干涉相位圖像，M4S模型仿真計算的主要流程如圖1所示。
圖1 M4S模型計算流程圖Fig.1 The flow chart of M4S model
2.2經典海表流場迭代反演算法
基于M4S模型，Romeriser和Hartmut[9]提出了順軌干涉SAR海表流場迭代反演算法，于祥禎等[7，13]給出了具體的海表流場校正方案，如圖2所示。
圖2 經典的海表流場迭代反演算法流程圖[7，13]Fig.2 Classical iterative inversion algorithm flow chart of the ocean surface current[7，13]
圖2中，un表示最終反演流場，n表示迭代次數，T表示順軌干涉相位的均方根誤差閾值。u0代表初猜流場，通過式(1)確定：
(1)
式中，λ表示雷達波長，V表示平臺飛行速度，L表示有效基線長度，θ表示雷達入射角，φ0表示實測順軌干涉相位。
具體流程為：
(1)首先，利用式(1)通過實測順軌干涉相位計算得到初猜流場u0。
(2)初始化海表流場迭代反演所涉及的參數。
(3)將風場以及得到的海表流場、平臺參數、雷達參數等輸入M4S模型，計算順軌干涉相位圖像。
(4)由于流場方位向梯度的存在[7，13]，需對得到的干涉相位圖像進行方向偏移的校正和均值濾波處理。
(5)然后將仿真的干涉相位圖像與實際干涉相位圖像進行對比，如果干涉相位的均方根誤差rmsen
(6)將仿真和實際的干涉相位進行逐個比對，通過設置合理的校正系數對干涉相位偏差大于閾值T的相位對應的流場值進行修正，將修正后的流場重新輸入M4S模型，進行下一次迭代校正。
通過實際的迭代計算發現，該方法雖能夠有效的反演出海表流場，但是其對流場校正系數所涉及的比例因子始終設置為一個固定值，這樣就不能很好的根據流場迭代的次數和當前精度實現系數的自適應調整，從而減緩了迭代的收斂速度。另外，算法針對的是目標區域所有的干涉相位值和流場值，沒有考慮到各區域流場的大小和方向不同時對校正系數有不同的要求，這也會影響海表流場反演的精度。
因此，本文在該算法的基礎上，針對校正系數和區域劃分，開展兩個方面的研究和改進。
3改進的海表流場反演算法
3.1基于模擬退火算法優化校正系數
模擬退火算法(Simulated Annealing，SA)的思想最早是由Metropolis于1953年提出的[15]，它來源于固體退火降溫的原理，即將固體加溫至充分高，再讓其慢慢冷卻這一過程。該過程可以解釋為：當對固體加熱時，其內部粒子的熱運動劇烈程度會隨溫度的增大而增大，此時粒子呈現出無序狀，內能也隨之增大；若停止加熱，則固體內部粒子會隨溫度的降低逐漸趨于有序，其熱運動的劇烈程度也逐漸降低，最后在常溫時達到基態，此時內能減為最小[16]。模擬退火算法在實際應用中的基本思路是：首先在高溫下進行搜索，由于這個時候各個狀態出現的幾率相差不大，可通過粗搜索大致找到低能區域；隨著溫度的降低，各狀態出現概率的差距逐漸拉大，搜索精度也不斷提高，進而越來越準確的找到目標函數的全局極小點，也即問題的最優解[17]。
要將模擬退火算法具體應用到流場的校正過程中，需要找到合適的參量來代替固體退火降溫這一過程中的內能和溫度。本文用目標函數值代替固體退火過程中的內能，用校正系數中的比例因子β代替固體退火過程中的溫度。具體思路為：首先，根據校正區域的流場特點給定一個初始的比例因子β0，并在其鄰域中隨機產生一個新的值；然后擬定一個接受準則，該準則具有一定容錯性，即允許目標函數在一定范圍ζ內接受使目標函數惡化的β，算法持續進行“產生新β值—計算目標函數值—判斷是否接受新的β—接受或舍棄”的迭代過程[16]。經過多次β值的變化后，可以求得流場校正問題的相對最優解。接著，逐步減小比例因子β的變化范圍，重復執行上述迭代過程，當這個范圍趨于零時，系統也越來越趨于平衡狀態，也即反演流場趨于最優流場，此時停止迭代并輸出該流場，如圖3所示。
圖3 基于模擬退火算法優化校正系數流程圖Fig.3 The flow chart of optimizing correction factor based on simulated annealing algorithm
具體步驟如下：
(2)
(3)
(2)初始化參數。除了給定初始的比例因子β0以外，還包括比例因子的變化范圍Γ0([-x0,x0])、衰減因子τ(負值)、最大迭代次數N、終止條件S和目標函數可以接受的范圍ζ(包括ζ1和ζ2)。
(3)判斷是否接受新的β值。當目標函數滿足下面的任意一個條件時，接受此時的β值；否則，舍棄該β值并在變化范圍Γ0內隨機產生新的β并重新進行判斷。
(4)
(5)
(4)如果滿足終止條件S或達到最大迭代次數N，則迭代終止；否則，令x0=x0+τ，并重復步驟(3)。其中終止條件S為同時滿足式(4)和式(5)。
3.2探測區域劃分的流場優化方案
在進行流場迭代校正前，根據流場大小和方向的不同，將目標區域進行劃分，結合前文提出的優化校正系數的思想，對每一塊區域分別進行迭代計算，最后再將各區域的輸出流場合成為最終輸出流場，如圖4所示。
圖4 劃分探測區域的流場校正方案Fig.4 Current correction scheme for the division of detection area
在對目標區域進行劃分時，主要涉及到流向和流速兩個參量，風場的信息未作考慮。
由試驗結果可知，脆口蘿卜的各包裝材質產品在各溫度梯度儲存期間理化指標pH和總酸都基本穩定，產品中細菌總數均<10 cfu/g，符合產品設計標準，但是隨著產品儲存期的延長，產品的色澤、香氣、口感都呈下降趨勢，通過Q10方法測算同種包裝材質，不同儲存溫度影響產品貨架期，高溫縮短了產品貨架期；不同種類的包裝材質對脆口蘿卜產品貨架期影響很大，鋁箔袋產品能夠達到324天，鍍鋁袋產品229天，透明袋產品最少，只有95天。不同的包裝材質，脆口蘿卜貨架期差異大，與包裝材質的氧氣透過率、吸濕性以及產品配料特性有關。
如圖5所示，將目標區域根據流場分布特點[18]劃分為9塊區域(這里為方便說明，將其等分，實際計算過程中并非完全均等劃分)，每塊區域的流場大小和方向閾值依具體的流向和流速確定(實驗中以流速差的絕對值不大于0.5 m/s，流向差的絕對值不大于6°為劃分標準)。
圖5 對目標區域進行劃分Fig.5 Division of detection area
3.3改進的海表流場反演算法
在經典的海表流場反演算法基礎上，作上述兩處優化，從而得到了改進后的海表流場反演算法，如圖6所示。具體步驟如下：
(1)首先利用式(1)通過實測順軌干涉相位計算得到初猜流場u0。
(2)根據劃分標準將探測區域劃分為若干個子區域并對每一個子區域進行下面的反演迭代。
(3)輸入風場、平臺參數和雷達參數等，進行海表流場的迭代反演。
(4)將仿真的干涉相位圖像與實際干涉相位圖像進行對比，如果干涉相位的均方根誤差rmsen
(5)將仿真和實際的干涉相位進行逐個比對，在對校正系數進行優化的基礎上，修正干涉相位偏差大于閾值T的相位對應的流場值，將修正后的流場重新輸入M4S模型，進行下一次迭代校正。
(6)將得到的每一個子區域的“最優流場”進行合成，得到最終的反演流場。
圖6 改進后的順軌干涉SAR海表流場迭代反演算法Fig.6 Improved iterative inversion algorithm flow chart of the ocean surface current
4結果與分析
4.1仿真數據及實驗
為了驗證改進后算法的性能，進行如下仿真和反演實驗。
首先，模擬一個二維海表流場(100×100，空間間隔為50 m)作為“真實流場”，海面背景風向為35°，風速為3 m/s。為了獲取二維的海表流場，假設星載SAR(右側視)沿垂直方向探測(按照雷達視向分別記為x和y方向)，探測方案如圖7所示，紅色箭頭表示海表流場。
圖7 海表流場及星載SAR探測方案Fig.7 Ocean surface current and space-borne SAR detection scheme
設置基于星載順軌干涉SAR的參數，包括中心頻率、雷達波長、極化方式、雷達入射角、飛行高度和速度、有效基線長度等，如表1所示。
然后，利用M4S模型計算初始相位φ0，假設其為“真實相位”(實際探測中，初始相位應為干涉SAR實測結果)，M4S模型計算的順軌干涉相位圖像如圖8所示。接著，利用式(1)，通過初始相位φ0計算初猜流場u0，并分別利用經典算法和本文改進的海表流場反演算法反演流場。最后，將反演流場與模擬的海表流場(真實流場)進行比對，并對兩種反演算法進行分析。
表1 星載順軌干涉SAR參數設置
圖8 M4S模型計算的初始順軌干涉相位Fig.8 The initial phase of along-track interferometric SAR calculated by M4S model
4.2反演結果分析
經過計算，海表流場的反演結果如圖9、圖10和圖11所示。其中，圖9描述了反演流場與給定流場(真實流場)的比對分析情況，圖10a和b分別描述了經典算法和改進后算法海表流速的反演結果與真實值比對分析情況，而圖11a和b則分別描述了經典算法和改進后算法海表流向的反演結果與真實值比對分析情況。
圖9 反演流場(a)與真實流場(b)Fig.9 Inversion of flow field (a) and real flow field (b)
圖10 經典算法(a)和改進算法(b)的海表流速反演值與真實值比對分析情況Fig.10 Sea surface velocity inversion and real value analy-sis of classical algorithm (a) and improved algorithm (b)
從圖10可知，經典算法海表流速反演結果的均方根誤差為0.048 m/s，偏差為-0.015 m/s；改進算法的海表流速反演結果的均方根誤差為0.016 m/s，偏差為-0.007 m/s。
從圖11可知，經典算法海表流向反演結果的均方根誤差為4.616°，偏差為1.670°；改進算法的海表流向反演結果的均方根誤差為0.486°，偏差為-0.334°。
分析反演結果可知，改進后的算法較經典算法在流場速度的精度上有較大提高，且反演流場與真實流場的符合程度也更高。其中，海表流速的均方根誤差優于0.02 m/s；海表流向的均方根誤差優于1.0°。
此外，采用經典海表流場反演算法反演海表流場時，迭代次數通常為6～8次，而采用模擬退火算法優化校正系數后，通常迭代2～3次。因此，改進的海表流場反演算法也加快了海表流場的收斂速度，減少了處理時間。
5小結
本文首先詳細介紹了于祥禎等提出的順軌干涉SAR海洋表面流場迭代反演算法，并簡要的分析了該方法存在的不足。然后針對這些不足給出了基于模擬退火算法優化校正系數的具體流程和區域劃分的流場優化方案。經過實驗對比分析可以看出，改進后的流場迭代反演算法較經典的反演算法在迭代的收斂速度和反演的精度上都有所提高，其中本次實驗的迭代次數較之前減少了約65%，流速和流向的偏差較之前分別減少了約53%和79%，證明該改進算法可行有效。
隨著順軌干涉SAR海洋觀測數據的積累以及利用M4S模型反演海表流場算法研究的深入，海洋表面流場的測量精度將進一步提高，這對監測海流和通過海流計算其他海洋參量如海浪譜等具有重要意義。
圖11 經典算法(a)和改進算法(b)的海表流向反演值與真實值比對分析情況Fig.11 Sea surface flow direction and real value analysis of classical algorithm (a) and improved algorithm (b)
參考文獻：
[1]劉巍， 張韌， 王輝贊，等. 基于衛星遙感資料的海洋表層流場反演與估算[J]. 地球物理學進展， 2012， 27(5): 1-4．
Liu Wei， Zhang Ren， Wang Huizan， et al. Sea surface flow field retrieval and estimation based on satellite remote sensing data[J]. Progress in Geophysics，2012，27(5):1-4．
[2]馮貴平， 金雙根， Jose M，等. 利用衛星測高、GRACE和GOCE資料估計全球海洋表面地轉流[J]. 海洋學報， 2014， 36(9):45-55．
Feng Guiping， Jin Shuanggen， Jose M， et al. Global ocean surface geostrophic currents estimated from satellite altimetry， GRACE and GOCE[J]. Haiyang Xuebao， 2014，36(9):45-55．
[3]常亮， 高郭平， 郭立新. 星載SAR海洋表層流場反演綜述[J]. 海洋科學進展， 2015(1):1-4.
Chang Liang， Gao Guoping， Guo Lixin. Review on ocean surface current field measurement by space-borne SAR[J]. Advances in Marine Science，2015(1):1-4．
[4]陳大可， 許建平， 馬繼瑞，等. 全球實時海洋觀測網(Argo)與上層海洋結構、變異及預測研究[J]. 地球科學進展， 2008， 23(1):1-7．
Chen Dake， Xu Jianping， Ma Jirui， et al. Argo global observation network and studies of upper ocean structure， variability and predictability[J]. Advances in Earth Science，2008，23(1):1-7．
[5]Paduan J， O’Donnell J， Allen A， et al. Surface current mapping in U.S. coastal waters: Implementation of a national system[R]. Arlington: Ocean， US， 2004．
[6]Barrick D E,Evans M W,Weber B L. Ocean surface current mapped by radar[J]. Science， 1977， 198(4313):138-144．
[7]于祥禎. 順軌干涉SAR對海洋表面流場監測的若干問題研究[D]. 北京：中國科學院研究生院， 2012．
Yu Xiangzhen. Study on some problems of ocean surface current observation by along-track interferometric SAR[D]. Beijing：Graduate University of Chinese Academy of Sciences，2012．
[8]Goldstein R M， Zebker H A. Interferometric radar measurement of ocean surface currents[J]. Nature， 1987， 328(6132)：707-709．
[9]Romeiser R， Hartmut R. Theoretical evalution of several possible along-track InSAR modes of Terra SAR-X for ocean current measurements[J]. IEEE Trans Geosci Remote Sens， 2007， 45(1):21-35．
[10]Romeiser R. Current measurements by airborne along-track InSAR: measuring technique and experimental results[J]. IEEE Journal of Oceanic Engineering， 2005， 30(3):552-569．
[11]Romeiser R， Schwabish M， Schulz-Stellenfleth J， et al. Study on concepts for radar interferometric from satellites for ocean (and land) applications (KoRIOLiS)[R].Germany：University of Hamburg， 2002．
[12]Romeiser R， Suchandt S， Runge H， et al. First analysis of TerrSAR-X along-track InSAR-derived currents fields[J]. IEEE Trans Geosci Remote Sens， 2010， 48(2):820-829．
[13]于祥禎， 種勁松， 洪文. 順軌干涉SAR海洋表面流場迭代反演算法[J]. 電子與信息學報， 2012(11):2660-2665．
Yu Xiangzhen， Zhong Jinsong， Hong Wen. An iterative method for ocean surface current retrieval by along-track interferometric SAR[J].Journal of Electronics and Information Technology， 2012(11):2660-2665．
[14]Romeiser R. M4S 3.2.0 user’s manual[S].Germany:University of Hambury, 2008．
[15]Metropolis N， Rosenbluth A W, Rosenbluth M N， et al.Equation of state calculations by fast computing machines[J]. Journal of Chemical Physics， 1953， 56(21):1087-1092.
[16]龐峰. 模擬退火算法的原理及算法在優化問題上的應用[D]. 長春：吉林大學， 2006．
Pang Feng. The principle of SA algorithm and algorithm’s application on optimization problem[D]. Changchun：Jilin University，2006．
[17]朱凱. 精通MATLAB神經網絡[M]. 北京:電子工業出版社， 2010．
Zhu Kai. Proficient in MATLAB neural network[M]. Beijing：Publishing House of Electronics Industry， 2010．
[18]陳麗娜. 基于海洋流場的拓撲區域劃分的研究[J]. 電子設計工程， 2011(8): 10-12.
Chen Li’na. Research on topological region division based on ocean flow field[J]. Electronic Design Engineering， 2011(8): 10-12．
A study of ocean surface current inversion algorithm based on simulated annealing and region partition
Zhang Xiarong1，Huang Yunxian1，Yan Wei1，Zhao Xianbin1
(1.InstituteofMeteorologyandOceanography，PLAUniversityofScienceandTechnology，Nanjing211101，China)
Abstract:The classical ocean surface current iterative inversion algorithm used constant correction factor to analyze and calculate the whole detection area， so it is time-consuming and the accuracy is limited. In order to improve the efficiency and accuracy of the inversion calculation， this paper used simulated annealing algorithm to optimize the correction coefficient， which can be adaptively changed according to the result of each inversion. Moreover， detection region was divided according to the characteristics of the current distribution. Finally， a modified ocean surface current inversion algorithm is proposed based on the classical method. The simulation results show that the efficiency and accuracy of the improved inversion algorithm are obviously improved．
Key words:Along-Track Interferometric SAR; ATI-SAR; Ocean surface current; M4S model; simulated annealing algorithm
收稿日期：2015-08-17；
修訂日期：2016-04-05。
基金項目：國家自然科學基金面上項目(41375029，41575028)；國家自然科學基金青年基金項目(41505016)；裝備預研共用技術基金(9140A21070114JB25344)。
作者簡介：張夏榮(1991—)，男，甘肅省定西市人，研究方向為微波海洋遙感。 *通信作者：黃云仙，女，上海市人，副教授，從事氣象信息處理技術研究。E-mail:gsdxzxr@163.com
中圖分類號:P731.21
文獻標志碼：A
文章編號：0253-4193(2016)07-0014-08
張夏榮，黃云仙，嚴衛，等. 基于模擬退火和區域劃分的海表流場反演算法研究[J]. 海洋學報， 2016， 38(7): 14-21， doi:10.3969/j.issn.0253-4193.2016.07.002
Zhang Xiarong， Huang Yunxian， Yan Wei， et al. A study of ocean surface current inversion algorithm based on simulated annealing and region partition[J]. Haiyang Xuebao， 2016， 38(7): 14-21， doi:10.3969/j.issn.0253-4193.2016.07.002
