基于路徑柵格的機場噪聲等值線追蹤算法

徐 濤，曹枝東

(1. 中國民航大學計算機科學與技術學院 天津 東麗區 300300; 2. 中國民航大學中國民航信息技術科研基地 天津 東麗區 300300)

隨著民航業的發展以及人們環境意識的逐步提高，機場噪聲問題日益突出。合理規劃使用機場周圍的土地是解決機場噪聲問題的有效手段，而機場噪聲等值線圖則是確定機場噪聲對居民的影響范圍、控制機場噪聲以及合理規劃機場周圍土地使用的重要依據。因此，繪制機場噪聲等值線圖是機場噪聲控制工作和機場規劃設計的重要環節[1]。

等值線圖是數字高程模型的表示形式之一[2]，其繪制一般分為以下步驟：離散數據點網格化、開放等值線起始網格單元的確定及追蹤、封閉等值線起始網格單元的確定及追蹤、等值線光滑、等值線填充[3]。文獻[4]介紹了一種柵格右關聯標志的規則網格數據等值線繪制算法，該算法雖避免了等值線交叉，但在等值點追蹤上對某個特定網格點需查看與其相鄰的4個網格單元，在存儲空間和搜索方式上具有一定的缺陷，算法效率不高。文獻[5]提出了一種基于搜索圓法的等值線追蹤算法，該算法具有快速掃描、追蹤效率高等優點，但仍需要先追蹤開放等值線，后追蹤封閉等值線，如網格邊界不規則，該方法效率會下降。文獻[6]針對數字地形圖中等值線的提取提出了一種改進的摩爾相鄰追蹤算法，可以完整地追蹤到整個地形圖等值線，但對同一個網格單元的重復追蹤需回溯到前一個網格單元，一旦很多網格單元被回溯時，算法效率急劇降低。機場噪聲問題過去幾十年受到極少關注，傳統的環境噪聲等值線圖大多結合現有軟件Surfer加載噪聲數據繪制而成[7]，而機場噪聲具有影響范圍廣以及因航跡覆蓋范圍不同而噪聲分布差異性大等特點[8]。

基于上述原因，本文提出了一種基于路徑柵格的機場噪聲等值線追蹤算法，通過構造有效網格、生成路徑柵格，利用入口方向和路徑柵格頂點編碼和準確判斷等值線走向趨勢，唯一確定下一個等值點，避免封閉等值線和開放等值線[9]的分別追蹤，可以快速生成機場噪聲等值線圖。

1 基本概念

1.1 等值線基本性質

等值線需滿足以下要求[10]：

1) 通常為一條光滑連續曲線；

2) 給定高程值相應域上的等值線不限于一條；

3) 等值線可以是封閉或者開放的；

4) 不能相互交叉。

從這些特點可以看出，對于特定高程值的等值線(不限于一條)，與特定網格單元的交點個數可能為0、2、4[11]，與特定網格單元的一條邊上的交點個數為0或1。

1.2 頂點編碼

不失一般性，以規則地形網格為例，一個規則地形網格水平和垂直方向分別由M和N個等距離排列的點組成，其中整個網格用G表示，單個網格單元用cells[i,j]表示，網格單元頂點用pts[i,j]表示，其中網格單元的坐標為該網格單元右下角頂點的坐標，如圖1所示。當某條等值線的高程值與網格單元的頂點屬性值相等時，等值線剛好經過該點。為避免這種情況對后續算法造成不便，沿襲一貫做法，在該網格單元頂點屬性值上加一個微小的數e，使等值線不直接經過網格單元頂點。

頂點編碼是特定高程值的等值線繞過網格單元的頂點而得到約定的固定數值。任意網格單元只有4條邊(不包括頂點)與等值線相交，以網格單元的4個頂點為參考，如圖2a～圖2d所示。當網格單元頂點A、B、C、D的相鄰兩邊與同一條等值線相交時，約定該網格單元的頂點編碼為1、2、4、8。當一個網格單元的兩個頂點被特定高程值的等值線繞過時，將網格單元的頂點編碼相加得到網格單元的路徑柵格頂點編碼和。

圖1 規則網格

1.3 路徑柵格頂點編碼和

相對于某個特定高程值，任意網格單元與等值線的交叉關系只有3種情況：

1) 沒有等值線經過：約定路徑柵格頂點編碼和PointCodeSum= 0。

2) 有一條等值線經過：每個網格單元中只有一個頂點被等值線環繞，如圖2a～圖2d所示，于是網格單元的路徑柵格頂點編碼和分別為1、2、4、8。

3) 有兩條等值線(不同連續等值線上的片段)經過[11]：網格單元分別有兩個頂點被等值線環繞，如圖2e、圖2f所示，網格單元的路徑柵格頂點編碼和分別為1+4=5和2+8=10，在圖2g、圖2h狀況下，分別約定以AB點和BC點為參考，網格單元的路徑柵格頂點編碼和分別為1+2=3和2+4=6。

圖2 頂點編碼及路徑柵格頂點編碼和

2 基于路徑柵格的機場噪聲等值線追蹤算法

2.1 關鍵步驟

2.1.1 有效網格的構造

機場及其附近區域的噪聲都是由航班的飛行引起，然而機場四周并非所有地方都有航班經過，距航跡較遠的地方噪聲較小(可以不考慮)。圖3為一個100*100的網格(其中淺色部分無噪聲，深色部分有噪聲)，U、V、W、X、Y、Z6個區域無噪聲，剔除這些區域的網格，于是有效網格為S=(G-U-V-W-X-Y-Z)。

圖3 機場噪聲有效網格

2.1.2 路徑柵格的生成

路徑柵格(相對特定高程值)是在有效網格中剔除沒有等值線經過的網格單元后剩余的網格部分[12]，其中每個網格單元都有路徑柵格頂點編碼和(不為0)，它是進行等值線追蹤的依據。路徑柵格的生成步驟如下。

1) 從高程值列表中選取一個高程值Hi，依次掃描有效網格，對每個網格單元做如下處理：分別計算Hi與網格單元cells[i,j]4個頂點(A、B、C、D)屬性值Attribute的差值，記為s1、s2、s3、s4，若有為0的情況，在差值上加一個微小的數e，并令F=s1s2s3s4。

①s1、s2、s3、s4都大于0或都小于0時，該網格單元的路徑柵格頂點編碼和cells[i,j].PointCodeSum=0。

②F<0(如圖2a～圖2d情況)，做如下判斷：

③F>0和s1s3>0(如圖2e、圖2f情況)，為避免兩條等值線交叉情況的出現，這里判斷等值線的走向需取網格單元的中心點P，令P點屬性值為centerAttr，則：

2) 遍歷完畢有效網格，高程值Hi的路徑柵格形成。

2.1.3 等值線的追蹤

基于生成的路徑柵格，構造網格單元出口方向判斷表，如表1所示。通過表1確定網格單元的出口方向。在等值線追蹤判斷出口方向時需要入口方向和路徑柵格頂點編碼和，確定了起始追蹤網格單元，可以得到路徑柵格頂點編碼和。如該網格單元的兩邊或者四邊有等值線穿過，可任取其中一個方向作為入口方向。

表1 網格單元出口方向判斷表

結合表1，給出一個等值線追蹤的示例。給定高程值Hi=30，已知起始追蹤網格單元cells[1,10]，如圖4所示，該網格單元PointCodeSum=3，上下兩條邊有等值線穿過。1) 取inDirection=上，由表1知，該網格單元outDirection=下，等值線進入網格單元cells[2,10]中，該網格單元PointCodeSum=10，inDirection=上，由表1知，該網格單元outDirection=右，LL，等值線進入網格單元cells[1,13]中，該網格單元PointCodeSum=3，inDirection=下，由表1知，該網格單元outDirection=上。等值線追蹤到邊界，Hi=30的一條等值線追蹤完畢。2) 取inDirection=下，則outDirection=上，等值線追蹤到達邊界，而此時開放的等值線并沒有追蹤完畢(兩端頂點都為邊界點的等值線為開放等值線)，需記錄等值線起始追蹤網格單元，記Begin[i,j]=cells[1,10]，當追蹤沒完成，返回起始追蹤網格單元，更換進入方向繼續追蹤，直至到達邊界。

圖4 等值線追蹤示例

在等值線追蹤過程中，需確定等值線與網格單元邊的交點。已知網格單元邊的兩個頂點(P1、P2)和特定高程值Hi，可以通過線性插值[11]來確定交點Pcross坐標，追蹤完畢順次連接交點即可繪制需要的等值線。

2.2 算法描述

綜上討論，對于規則的機場噪聲數據，可以通過路徑柵格等值線追蹤算法進行繪制，算法完整描述如下：

1) 構造有效網格。讀取機場規則噪聲數據，將網格單元頂點保存為二維數組pts，網格單元保存為二維數組cells。剔除明顯沒有等值線經過的網格區域，形成初始遍歷的有效網格S。

2) 生成路徑柵格。對于特定的等值線高程值Hi，遍歷S區域網格單元，確定每個網格單元的路徑柵格頂點編碼和。若PointCodeSum>0，將網格單元插入路徑柵格鏈表CellsArray末尾。

3) 確定等值點。取CellsArray首元素，以此為等值線起始追蹤網格單元，利用路徑柵格頂點編碼和任取一個可能的入口方向，記錄起始網格單元(Begin[i,j])和起始入口方向，依次確定下一個等值點，通過線性插值計算等值線與網格單元邊的交點并記錄，順次存入等值點鏈表Path中(將起始追蹤網格單元的入口方向邊的等值點作為Path首元素)。

4) 追蹤判斷等值點。若正在追蹤的等值點到達邊界并且起始等值點Path[0]為邊界點，說明等值線為開放等值線，當前等值線追蹤完畢；若正在追蹤的等值點到達邊界但Path[0]不為邊界點，則重新以Begin[i,j]為起始追蹤網格單元，更換入口方向轉入步驟3)繼續追蹤，直至追蹤的等值點到達邊界；若正在追蹤的等值點與Path[0]相同，說明等值線為封閉等值線，當前等值線追蹤完畢。在追蹤過程中，每次遍歷一個網格單元，取當前被等值線繞過的頂點編碼number，計算PointCodeSum= PointCodeSum-number，若PointCodeSum=0，則從CellsArray中剔除網格單元。

5) 依次連接Path中相鄰兩點，高程值Hi的一條等值線繪制完畢。

6) 重復步驟3)～步驟5)，直至路徑柵格鏈表CellsArray長度為0(高程值Hi的所有等值線繪制完畢)。

7) 重復步驟2)～步驟6)，直至預先設定的等值線高程值Hi集合遍歷完畢。

3 實驗結果及分析

實驗使用兩個機場的噪聲數據集，共選取6組實驗數據，分別是規模為100*100、200*200、300*300、300、400*400、1000*1000、1 989*1 989的規則地形網格數據。選取機場噪聲分貝值層數分別為10、30、66。其中10層的噪聲分貝值數據為：{40,45,50,55,60,65,70,75,80, 85}，30層是從30 dB到88 dB范圍以2 dB為間隔的噪聲分貝值數據構成的集合，66層則是從20 dB到85 dB范圍以1 dB為間隔的噪聲分貝值數據構成的集合。實驗機器配置：Intel? Core?2 Duo CPU E4600 2.40 GHz，內存2.00 GB，操作系統為Microsoft Windows XP Professional 2002 Service Pack 3。

本文算法繪制的等值線圖實驗結果如圖5、圖6所示，圖5是數據規模為400*400層數為10的某單跑道機場噪聲等值線圖，圖6是數據規模為100*100層數為30的某雙跑道機場噪聲等值線圖。從圖中可以看出，即使在噪聲分貝值很密集的情況下，繪制的等值線圖仍能夠很好地符合等值線的基本特點，特別是沒有等值線交叉的情況出現。圖5中心的短粗線部分為機場跑道，沿著跑道兩端的圓點曲線為該機場主要航跡，從中可以看出，單跑道機場的噪聲等值線以跑道為中心，依附航跡分散開來，任意兩條相同分貝值間隔的相鄰等值線，距離航跡越近的區域等值線越密集，噪聲分貝值也越大。雙跑道的情況與單跑道情況類似，只是介于兩條跑道航跡之間的區域噪聲分貝值會被疊加，其值大于單跑道航跡上該地的噪聲分貝值(如圖6中的A點)。

此外，本文選取3種算法(傳統的等值線傳播算法，文獻[5]中的路徑柵格邊界追蹤算法和本文算法)對上述6組實驗數據3種層次的噪聲分貝值進行算法時間效率對比，其結果如表2所示。從表中可以看出，數據規模為100*100、200*200、300*300、400*400時，本文算法時間效率總體上優于傳統的等值線傳播算法，明顯優于路徑柵格邊界追蹤算法；當數據規模增大至1 000*1 000和1 989*1 989時，與傳統的等值線傳播算法和路徑柵格邊界追蹤算法相比，本文算法效率優勢更為明顯。

圖5 某單跑道機場噪聲等值線圖

圖6 某雙跑道機場噪聲等值線圖

表2 3種算法效率比較

本文算法時間效率較高的一個重要原因就是路徑柵格鏈表的使用。算法在構造路徑柵格時，將特定高程值等值線經過的所有網格單元以鏈表存儲，形成路徑柵格鏈表，如圖7所示，圖中深色部分為等值線經過的網格單元，等值線追蹤時只需從路徑柵格鏈表中取出首元素(某個深色網格單元)開始追蹤即可，從而避免常規的先掃描邊界網格單元追蹤開放等值線，再掃描內部網格單元追蹤封閉等值線。因此通過預先構造路徑柵格，建立路徑柵格鏈表，避免了開放等值線和封閉等值線的分開追蹤，提高了等值線追蹤效率。

圖7 路徑柵格

4 總 結

等值線追蹤算法很多[13-14]，對于不同的應用環境需要考慮不同的具體情形?；诼窂綎鸥竦臋C場噪聲等值線追蹤算法首先根據機場噪聲數據的特殊情況建立有效網格，掃描有效網格確定網格單元的路徑柵格頂點編碼和，從而建立路徑柵格，在初始追蹤網格單元和入口方向確定的情況下，通過路徑柵格頂點編碼和唯一確定下一個等值點，并且路徑柵格的提前建立避免了傳統的先開放后封閉的等值線追蹤方式。對比實驗可以看出，本文算法在時間效率上具有明顯優勢，繪制的等值線圖很理想，沒有兩條等值線交叉情況出現。

[1] 李冉. 機場航空噪聲預測及其影響因素研究[D]. 天津:中國民航大學, 2008.LI Ran. Airport air noise prediction and influencing factors study [D]. Tianjin: Civil Aviation University of China, 2008.

[2] KIDNERA D B. High-order interpolation of regular grid digital elevation models[J]. International Journal of Remote Sensing, 2003, 21(14): 2981-2987.

[3] 苗潤忠. 光滑的等值線生成算法[J]. 長春理工大學學報,2004, 27(1): 16-18.MIAO Run-zhong. A new smoothed contour lines generating algorithm for quadrilateral meshes[J]. Journal of Changchun University of Science and Technology, 2004,27(1): 16-18.

[4] 王結臣, 錢晨暉, 芮一康. 柵格數據生成等值線的一種實用方法[J]. 測繪科學, 2007, 32(6): 88-90.WANG Jie-chen, QIAN Chen-hui, RUI Yi-kang. A practical method for contour generation based on raster data[J].Science of Surveying and Mapping, 2007, 32(6): 208-282.

[5] 常會, 柴華彬, 鄒友峰. 基于搜索圓法的等值線追蹤技術[J]. 測繪科學, 2009, 34(2): 119-121.CHANG Hui, CHAI Hua-bin, ZHOU You-feng. Isoline tracing technique based on search-circle[J]. Science of Surveying and Mapping, 2009, 34(2): 119-121.

[6] PRADHAN R, KUMAR S, AGARWAL R, et al. Contour line tracing algorithm for digital topographic maps[J].International Journal of Image Processing, 2010, 4(2):156-163.

[7] 過春燕, 張邦俊. 基于Surfer的機場噪聲等值線計算機繪制方法[J]. 中國環境科學, 2003, 23(6): 631-634.GUO Chun-yan, ZHANG Bang-jun. Method of drawing airport noise contours on computer based on Surfer[J].China Environmental Science, 2003, 23(6): 631-634.

[8] BERNARD M S, VAN P, BAARSMA B E. Using happiness surveys to value intangibles: the case of airport noise[J]. The Economic Journal, 2005, 115(500): 224-246.

[9] 孫桂茹, 馬亮, 路登平, 等. 等值線生成與圖形填充算法[J]. 天津大學學報, 2000, 33(6): 816-818.SUN Gui-ru, MA Liang, LU Deng-ping, et al. Investigation on the algorithm of making and filling isoline[J]. Journal of Tianjin University, 2000, 33(6): 816-818.

[10] 余明輝, 萬遠揚, 余飛. 一種繪制等值線圖的新方法[J].武漢大學學報(工學版), 2006, 39(3): 52-54.YU Ming-hui, WAN Yuan-yang, YU Fei. A new method of drawing isoline map[J]. Engineering Journal of Wuhan University, 2006, 39(3): 52-54.

[11] 張顯全, 劉忠平. 基于格網模型的等高線算法[J]. 計算機科學, 2005, 32(9): 199-201.ZHANG Xian-quan, LIU Zhong-ping. An algorithm of contour lines based on regular grid[J]. Computer Science,2005, 32(9): 199-201.

[12] JONES N L, KENNARD M J, ZUNDEL A K. Fast algorithm for generating sorted contour strings[J].Computers and Geosciences, 2000, 26(7): 831-837.

[13] CHANG F, CHEN C J, LU C J. A linear time component labeling algorithm using contour tracing technique[J].Computer Vision and Image Understanding, 2004, 93(2):206-220.

[14] ANGELOVA D, MIHAYLOVA L. Contour extraction from ultrasound images viewed as a tracking problem[C]//Proceedings of the 12th International Conference on Information Fusion. [S.l.]: [s.n.], 2009: 284-291.