基于對稱非負矩陣分解的終端區扇區劃分方法

張兆寧， 柯智舟

(1.中國民航大學空中交通管理學院， 天津 300300; 2.中國民航大學安全科學與工程學院， 天津 300300)

近年來，隨著航班數量的激增，終端區資源十分有限，一旦終端區扇區劃分不合理，極易造成管制超工作負荷，進而影響飛行的安全性。因此，若要使得管制負荷處于合理的區間，有必要對終端區扇區劃分方法進行研究，便于管制部門對扇區單元中航空器保障工作的整體把控。

許多國內外學者對扇區劃分方法進行了研究。在國外，2012年，Kumar[1]為了對相鄰扇區進行劃分，運用ART1神經網絡方法均衡了管制員工作負荷的目標需求。同年，Xue[2]基于利用遺傳算法對Voronoi
圖的扇區邊界進行優化，多次迭代確定最終的扇區劃分方案。2013年， Flener等[3]提出以管制負荷為約束條件，進行幾何扇區優化的方法。2014年，Flener等[4]運用遺傳算法對3D空域扇區進行劃設優化。2015年，Sergeeva等[5]歸納總結扇區劃分的分類標準和限制因素，提出了自動扇區劃分的概念，對自動扇區劃分做了調查研究。2019年，Granberg等[6]提出了一種譜聚類算法，利用 Voronoi剖分，結合整數編程，對航路和空域扇區進行了劃分及優化。2020年，Oktal等[7]結合地理信息系統分析空域，構建空域的區域劃分多目標整數規劃數學模型，對土耳其空域進行了扇區邊界的優化。

在國內，2009年，張明等[8]使用可變精度粗集理論，分析飛行數據，提出了不同航班數量下的動態扇區數目劃分規則。2013年，王超等[9]提出頂點間聯系程度模型，運用譜聚類算法實現扇區凸形有效分割。同年，尹文杰等[10]將終端區扇區劃分當作整數規劃問題，解得最優解來生成終端區扇區。2014年，羅軍等[11]通過構建拓撲結構模型，并使用粒子群算法解決了扇區優化問題。2015年，王超等[12]以均衡扇區的飛行流量為目的，運用動態規劃的方法，進行空域的劃分。2017年，王莉莉等[13]基于復雜度因子影響程度來定義管制員工作負荷，利用生長算法對扇區進行劃分優化。2019年，王興隆等[14]把扇區移交點作為研究對象，以最小航班延誤成本為目標，利用遺傳算法進行延誤成本的優化。2021年，張兆寧等[15]提出了一種基于態勢的動態優化算法。同年，朱承元等[16]定義交通管制復雜度，利用反向傳播(back propagation，BP)神經網絡對管制員工作負荷進行預測。同年，高偉等[17]根據空域關鍵點之間的靠近程度，建立帶權重的空域網絡結構模型，利用譜聚類實現了扇區的劃分。

可以看出，前人研究主要集中于利用幾何模型以及啟發式算法對扇區單元的劃設方法研究，然而使用聚類方法的研究較少。現通過增加惡劣天氣和航行情報因素，建立終端區劃分模型。以新的組合特征參數生成的 “snake”扇區序列構建對稱且非負的相似矩陣，實現聚類效果。為此，使用對稱非負矩陣分解(symmetric nonnegative matrix factorization， SNMF)應用于終端區劃分，以期實現特征提取和快速降維，提高時效性。

1 終端區劃分模型

分析終端區網絡加權圖，利用節點和邊集的權重計算管制員的工作負荷。在原管制負荷定義的基礎上增加考慮惡劣天氣和航行情報產生的負荷，建立以滿足各扇區間管制負荷均衡為目標的終端區劃分模型。

1.1 管制負荷和終端區描述

管制員工作負荷，即航空器運行過程對管制員造成了客觀的任務需求，管制員為完成監視、沖突、協調等任務對身體和精神上產生的消耗。扇區的劃分依據主要是對管制負荷進行量化，因此，定義第k扇區的管制負荷Wk為

(1)

假設終端區模型為二維平面進行劃分研究，不考慮高度的因素。在空間二維平面上可以用二元組G(V，E)表示，其中V和E分別為節點和邊的集合，以自然航路點為節點，航段為邊。終端區網絡圖則分為有向圖和無向圖，無向圖指的是兩個節點之間只有相連接和不相鏈接的兩種關系，即為至多存在一條邊，而有向圖則指的是兩點之間有兩個不同方向的連接線組成。采用無向圖來進行終端區的描述，無向圖中終端區結構的節點集合為V={V1，V2，…，VN}，邊集為E={E12，E13，…，Eij}。

在實際運行過程中，除了考慮終端區的物理結構，還需要對終端區的整體流量趨勢進行分析，因此，終端區劃分模型的構建還要考慮流量因素。以通過節點Vi的航空器架次記為權重wi，即圖1中圓圈中的數字2、4、7、8、9、10、11、12，通過邊Eij的航空器架次記為權重wij，即圖1中連線上的數字1、2、3、4、5、6。則終端區網絡加權圖可以用權重來表示流量，以管制過程中的航空器架次為測度標準的表達，用F={wi，wij}表示。若無向圖中具有此屬性，則可用三元組G(V，E，F)來表示，圖1所示為終端區網絡加權圖。

1.2 模型構建

從定義上看，管制負荷主要由監視負荷、沖突負荷和協調負荷三部分參數組成，終端區劃分需要量化各扇區內管制負荷，建立以均衡扇區間管制員工作負荷為目標的模型。

管制員對扇區k的監視負荷和沖突負荷可用在處理該扇區內經過航路點的飛機總架次的監視和沖突任務所需的時間來表示，飛機總架次可用包含在扇區k內的節點權重wi來表示，即

(2)

(3)

式中：TM、TC分別為處理監視和解決沖突任務所需時間；xik為決策變量。

管制員對扇區k的協調負荷指的是進入和離開扇區k所需的時間，即航空器進入和離開扇區k的飛機總架次移交任務所需的時間，飛機總架次可用權重wij來表示，即

(4)

式(4)中：TO為處理協調任務所需時間；xik、xjk均為節點決策變量，可表示為

(5)

終端區內管制員在日常工作中除了監視航空器、防止航空器產生沖突和跨區飛行航空器的協調工作外。若因風向的改變，引起跑道換向；因日常設備維護，導致部分導航設備關閉；因某條跑道維修，改用單跑道運行等情況，是終端區管制員所特殊需要考慮的問題。且終端區內航空器的進近和離場過程中，天氣因素的變化給管制員工作造成的困擾表現較為顯著。

因此在原來的管制員工作負荷中還考慮了惡劣天氣和航行情報兩個因素對管制負荷的影響。則扇區k內的管制負荷的表達式可更新為

(6)

終端區的劃分是以滿足各扇區單元之間的管制負荷均衡為目標。基于此對扇區單元進行聚類劃分，得到劃分目標函數為

(7)

式(7)中：Wi為扇區i的管制負荷；扇區數目|S|=NS。

扇區單元最終劃分個數的確定主要結合英國運籌與分析理事會提出的DORATASK[18]方法分析得出式(8)，扇區單元集合為S={S1，S2，…，SNS}，最終劃分扇區數目為|S|=NS。

(8)

2 模型求解方法

針對模型目標中兩兩扇區間的管制負荷均衡問題，可以考慮將不同管制負荷之間的關系組合成新的特征參數并用于構建“snake”扇區序列，通過不同扇區單元在“snake”中的不同位置來判斷兩兩扇區單元之間的相似度，構建相似度矩陣。使用對稱非負矩陣分解對相似度矩陣進行求解，最終聚類劃分的扇區滿足目標函數的最優解。

2.1 層次聚類分析

扇區的劃分過程中構建扇區序列需要對參數進行加權融合處理。利用專家打分法確定監視、沖突、協調負荷之間的權重大小，生成新的組合特征參數，計算公式為

(9)

(10)

式(10)中：σ為方差。

在組合特征參數的基礎上，參考文獻[19]中“snake”的概念。首先，對所有扇區單元運用層次迭代聚類算法進行分析。同時，管制負荷具有相關性，如果終端區空域中的一塊扇區在某一時刻發生了負荷過高情況，則會導致其相鄰扇區在同一時刻很快也發生同類情況。參照此想法，為了方便在第二步驟完成相似矩陣相似度的計算，在第一步驟聚類過程中，使用層次迭代聚類算法找到各個扇區單元之間相鄰且相似的高負荷情況的扇區。其中層次迭代聚類算法的順序流程圖如圖2所示。

由圖2可看出，終端區中每塊特定的扇區模塊，都會產生對應的一個扇區序列，由終端區中所有扇區組成，被稱為“snake”，其中第一個元素的扇區即為該特定扇區。然而每一條“snake”扇區序列的產生都是由特定扇區開始的，逐一吸收當前已加入的扇區的相鄰扇區，直到包含了整個終端區中所有的扇區模塊。

在層次迭代聚類算法分析各個扇區單元時，為了得到扇區單元所對應的組合特征參數，“snake”扇區序列首先要識別自己相鄰的所有扇區。同時，將相鄰的扇區與已經聚類的扇區平均值進行比較。為了在增加扇區單元時，確保扇區與扇區單元之間的低方差和空間相鄰性質，選擇一個與該扇區最為相似的加入當前扇區系列中。扇區間方差的計算公式為

圖2 層次迭代聚類算法Fig.2 Hierarchical iterative clustering algorithm

(11)

(12)

2.2 相似度矩陣計算

通過層次迭代聚類算法使得終端區內每個扇區單元模塊都有相對應的扇區序列“snake”。在所有“snake”扇區序列的基礎上，考慮兩兩扇區之間的負荷均衡性和空間相連性，在此基礎上執行兩兩扇區單元之間的相似度計算。結果得到每個扇區有且僅有一個由該扇區為首具有相連性和均衡性的扇區序列“snake”區域。因此，通過對所有不同長度尺寸下的扇區序列“snake”中包含的相同元素的數量進行統計，即為兩兩扇區單元之間的相似度值，可以用兩兩扇區序列“snake”之間相似度值來表示。計算公式為

(13)

式(13)中：w(i，j)為相似度矩陣W中的元素；sik、sjk分別為扇區i和扇區j之間對應的“snake”的前k個扇區；N為終端區中具有扇區個數的數量；intersect(sik，sjk)為sik和sjk在長度尺寸k的前提下包含的相同扇區的數量，其中，長度尺寸k指的是截取“snake”前k個扇區。

相似度矩陣為兩個不同扇區單元之間的相似度值組成的，是用來評判兩個不同扇區單元之間的相似程度，即為不同的扇區單元所對應的扇區序列“snake”之間的相似程度，其中該扇區在當前扇區序列“snake”中位置越靠后，則與位于“snake”的起始扇區之間越不相似，極易使得兩者之間的相似性不能完全表現出來。因此，為了可以更好地理解扇區序列“snake”中不同的扇區的位置順序，該扇區應該具有不同權值的原因如圖3所示。

圖3(a)中，sg和sh兩個扇區單元的扇區序列“snake”極為相似的原因，是因為{I，J，K，H}扇區單元從首次被累計計算進來后，每次統計相同扇區時，都會再次被統計一次，直到統計完成所有尺寸的扇區序列“snake”中相同扇區數目為止。因此，這些扇區會被統計很多次。比如，扇區I在sg和sh尺寸為3的時候被統計了一次，這時候尺寸為3的sg和sh有公共扇區數為1，即為{I}，下一步統計尺寸增長到4時，sg和sh的公共扇區數同增長到了2，即為{I，J}，加上之前計算得到的數據，得到統計的數值為3。以此類推，首先將之前統計得到的相同扇區的數目，疊加到當前的相似度值上，直到所有長度尺寸下的 “snake”中相同扇區被統計完時，扇區I被統計的次數為N-2次。

圖3(b)中，扇區{a，b，e}被統計的次數較少，由于處于靠后的位置。因此，sa和sz兩個扇區序列“snake”之間的相似程度不同。雖然扇區a出現在sa中第一位置，但是出現在sz中的倒數第二位，最后被統計到的總體次數為3。扇區之間都存在空間相連性，當扇區在不同的“snake”都處于靠前位置，則計算兩個“snake”相似度的時候，該扇區被統計的次數偏多。所以，式(13)可以用式(14)來代替，即

圖3 相似度計算Fig.3 Similarity calculation

(14)

式(14)中：siN和sjN分別為扇區i和扇區j所對應的完整“snake”；L(h，siN)為扇區h在siN完整“snake”的位置。

扇區在“snake”中的位置越靠前，則越表示它與“snake”之間的相似度越高，若是對兩條不同的扇區序列“snake”進行分析，其中所在位置都比較靠前，則可以間接表示兩兩之間的相似程度越高，同時也反映出扇區序列“snake” 對應的兩個扇區之間的相似度也就越高。通過引入了權重系數Φ≥1，可以得到更加緊湊的扇區模塊，權重因位置不同而具有不同的權值，使得位置對權重的影響更加明顯。

(15)

式(15)中：Φ變大，則更加關注空間位置信息，使得到的扇區模塊更加緊湊，從而可以簡化下一步中對稱非負矩陣分解的步驟。

2.3 對稱非負矩陣分解

對稱非負矩陣分解(SNMF)是一種圖聚類的方法，用構建的相似矩陣可以將原始的數據轉化為圖形，然后分解矩陣得到最后劃分結果。對稱非負矩陣分解主要通過對維度的規范，從而降維解決未解(non-deterministic polynomial，NP)難題。在智慧交通領域中，設計區域聚類劃分方法時常選用標準割圖(normalized cut，NC)作為優化目標函數，目標函數為

(16)

式(16)中：w∈RN×N為相似度矩陣；H∈RN×NS為聚類指示矩陣，其中每一行元素的最大值相對應的列名，則為該實例所屬的類別，N和NS分別為扇區個數和預先估計的最后劃分扇區個數，通常NS<

上述的目標函數的約束條件H≥0和HTH=I使得優化問題變成了NP難題，為了解決NP難題，SNMF可通過松弛約束變量進行近似求解，SNMF則松弛HTH=I保留H≥0，而約束條件中HTH=I的含義為每個對象有且只有一種類別，通過松弛該約束條件，使得每個對象可以屬于多個劃分結果，雖然這樣得到的結果不符合聚類目的，但是，只保留H≥0的約束條件下，同樣能使得最后的劃分結果滿足HTH≈I。松弛過后的目標函數則變為

(17)

使用公式計算兩兩扇區之間的相似度，構成相似矩陣w進行歸一化處理作為新的相似度矩陣。歸一化處理可以很好地限制了終端區扇區面積的大小，防止出現過大或者過小的扇區模塊，計算公式為

(18)

(19)

在進行對稱非負矩陣分解聚類劃分的過程中，確定最終劃分扇區數量NS，利用隨機初始化的方式生成H矩陣，得到滿足最小化目標函數所對應的低階近似矩陣。其中低階近似矩陣H的每一行代表的都是每個扇區屬于此區域的概率，概率最大的類別則為該扇區的歸屬區域。

3 算例分析

以某終端區空域，選取該終端區某日10：00—12：00時間段內的航班數據為依據，研究對象為現實運行過程中的206架航班，43條航段，27個航路自然點及導航臺， 11：16—11：46內有惡劣天氣的真實數據為例。圖4是基于該終端區的重要航路坐標點與終端區結構所繪制的原先終端區扇區單元Voronoi圖，則該終端區初步可由26個扇區單元組成。

以本文模型結合終端區扇區結構數據來計算該終端區內每個扇區單元的管制員工作負荷如圖5所示，并且表1給予了部分的扇區結構統計數據。

圖4 原先終端區扇區單元Voronoi圖Fig.4 Voronoi diagram of sector unit in the original terminal area

圖5 各扇區單元管制員工作負荷Fig.5 Controller workload of each sector unit

表1 部分扇區結構數據統計Table 1 Part of the sector structure data statistics

對于新的組合特征參數中3種管制負荷關系之間權重的確定，從華東空管局、中南空管局、天津空管局以及民航院校，選取10名專家。均具有豐富的空域扇區劃分研究經驗和長年的一線管制經歷。專家結合實際工作過程中的操作情況，以及科研調研過程中的數據分析對兩兩參數之間進行權重的確定。綜合考慮多名專家的判斷，多次反復溝通后對權重求平均值。最終確定α、β、γ分別為0.4、0.4、0.2，生成新的組合特征參數，對扇區進行聚類劃分，進行終端區的劃設。

根據本文方法，依據Voronoi生成的扇區單元之間進行聚類劃分，在基本參數中，統計得到管制總負荷為14 789 s，確定最終扇區個數NS=3，然后進行求解。最終得到如圖6所示可以滿足扇區間管制負荷均衡效果的劃分結果，01號扇區包含1～3、7～10、19、20、22扇區單元，02號扇區包含4～6、11、13、18、21、23、24、27扇區單元，03號扇區包含12、14～17、25、26扇區單元。聚類劃分后該終端區所劃分的3個扇區的管制員總工作負荷如表2所示。

終端區扇區初步劃分結束后，為了避免扇區邊界不符合扇區劃分原則，產生凹形的扇區邊界，將一個扇區邊界光滑技術應用于扇區初步劃分后的邊界，如圖7所示。

為了體現本文方法的優越性，比較遺傳算法、生長算法對于扇區劃分的收斂結果。SNMF劃分方法對于圖聚類相較于其他兩種方法具有良好的收斂性和收斂速度，如表3所示。

圖6 終端區扇區劃分結果Fig.6 Sector division result of terminal area

表2 聚類后各扇區管制負荷和空域利用率結果Table 2 Control load and airspace utilization results of each sector after clustering

圖7 光滑處理后的終端區扇區劃分結果Fig.7 Sector division result of terminal area after smoothing

表3 不同方法的扇區劃分結果對比結果Table 3 Comparison of sector division results of different methods

結果顯示，對于該時段的扇區劃分，以新的組合特征參數以及SNMF劃分方法的結果中，扇區間管制負荷和空域利用率基本均衡，優于劃分前的扇區單元，驗證了組合特征參數的有效特征，體現出本文方法的時效性和合理性。

4 結論

從終端區網絡圖出發，通過量化管制員負荷，建立扇區劃分模型。利用新的組合特征參數構建“snake”扇區序列，計算兩兩“snake”序列之間的相似度量，構建相似度矩陣并進行歸一化處理，下一步利用對稱非負矩陣分解方法對相似度矩陣進行分解，得到滿足扇區單元之間管制負荷均衡的聚類結果。從結果上來看，劃分方法具有良好的效果，結合某終端區的現實數據，展現了模型和方法的有效性，適用于實際應用，為終端區扇區劃分方法提供理論依據，使之更加合理。