利用概率松弛法的城市路網(wǎng)自動匹配

張云菲，楊必勝，欒學(xué)晨

1.武漢大學(xué) 測繪遙感信息工程國家重點實驗室，湖北 武漢430079；2.武漢大學(xué) 時空數(shù)據(jù)智能獲取技術(shù)與應(yīng)用教育部工程研究中心，湖北 武漢430079

1 引 言

隨著傳感器技術(shù)與地球空間信息科學(xué)的發(fā)展，空間數(shù)據(jù)的采集、建模、融合、更新等正處于快速發(fā)展時期。實際應(yīng)用中常需要對不同設(shè)備、不同部門采集的同一地區(qū)的多源、多尺度空間數(shù)據(jù)進行有效集成，以用于信息共享與變化檢測。空間數(shù)據(jù)匹配是實現(xiàn)空間數(shù)據(jù)變化檢測和更新的前提條件，其目的是建立同一區(qū)域內(nèi)兩份或多份數(shù)據(jù)間同名目標(biāo)的對應(yīng)關(guān)系，通過這種對應(yīng)關(guān)系可以實現(xiàn)兩者的屬性共享和幾何精度改進［1－2］，亦可對無對應(yīng)關(guān)系的目標(biāo)進行變化檢測和數(shù)據(jù)更新［3］。眾源（crowdsourcing）地理數(shù)據(jù)，如歐洲的“開放式道路地圖”（open street map，OSM）的出現(xiàn)，為空間數(shù)據(jù)的更新提供新的、免費的數(shù)據(jù)源。然而，由于其生產(chǎn)過程非標(biāo)準(zhǔn)化，眾源數(shù)據(jù)匹配給傳統(tǒng)的空間數(shù)據(jù)匹配方法帶來新的挑戰(zhàn)。

文獻［4］對空間數(shù)據(jù)匹配的研究現(xiàn)狀進行了系統(tǒng)的總結(jié)。對路網(wǎng)匹配而言，部分學(xué)者通過計算與候選路段間的幾何、拓?fù)洹⒄Z義等相似性特征選取最優(yōu)的匹配路段［5－9］。多源空間數(shù)據(jù)之間通常存在非剛性偏差，通過幾何、拓?fù)洹⒄Z義等相似性可以確定多條可能的候選路段，但很難進一步確定實際的匹配路段，因此，需要考慮在整個路網(wǎng)中與其他匹配對之間的兼容性。一些學(xué)者將相似性指標(biāo)融合成概率或其他綜合指標(biāo)，考慮鄰接目標(biāo)的匹配關(guān)系更新自身匹配概率［10－13］。文獻［11—12］嘗試?yán)酶怕仕沙诜ń鉀Q矢量道路匹配問題，實現(xiàn)了1∶1的道路節(jié)點對應(yīng)，但需要進一步建立路段間的對應(yīng)關(guān)系，文獻［13］將概率松弛法直接應(yīng)用于路段，需要設(shè)定多種閾值，且匹配概率受路網(wǎng)密度影響較大。文獻［11—13］忽略1∶0匹配，即無對應(yīng)關(guān)系的獨立路段，這些獨立路段是潛在的變化檢測目標(biāo)。另外，由于OSM等眾源數(shù)據(jù)的非標(biāo)準(zhǔn)化表達方式以及內(nèi)部尺度不一性，以不同數(shù)據(jù)作為參考得到的結(jié)果不一（即存在匹配方向性問題），且1∶N、M∶N的對應(yīng)關(guān)系很難有效識別。

概率松弛匹配方法在計算機視覺和模式識別領(lǐng)域是一種經(jīng)典的圖匹配方法［14］，該方法首先對兩份圖結(jié)構(gòu)建立初始匹配矩陣，再結(jié)合鄰域目標(biāo)的結(jié)構(gòu)關(guān)系啟發(fā)式地更新初始的匹配矩陣使之達到全局最優(yōu)。城市道路網(wǎng)是一種復(fù)雜的圖結(jié)構(gòu)，實現(xiàn)路網(wǎng)匹配需要綜合考慮位置、形狀、結(jié)構(gòu)等特征，采用概率松弛法能夠?qū)⒍喾N約束有效結(jié)合，實現(xiàn)路網(wǎng)的整體匹配，而不是局部映射。因此，本文在改進文獻［13］的概率和兼容函數(shù)計算模型基礎(chǔ)上提出一種基于概率松弛法的城市路網(wǎng)自動匹配方法，該方法能夠建立1∶0的匹配關(guān)系，減少閾值和路網(wǎng)密度的影響。最后，本文選取不同地區(qū)的兩份路網(wǎng)數(shù)據(jù)進行算法驗證。

2 基于概率松弛法的自動匹配方法

根據(jù)文獻［11—13］，概率松弛匹配方法包括概率矩陣初始化和迭代更新兩個基本過程。本文考慮路網(wǎng)的復(fù)雜性增加了匹配結(jié)果選取過程，具體算法流程如圖1所示，主要分為3步：①利用緩沖分析檢測候選匹配路段，然后根據(jù)候選路段間的幾何差異指標(biāo)計算不同匹配方向的概率并綜合得到初始概率矩陣；②由于初始概率未考慮鄰接路段對其的兼容程度，即兼容系數(shù)，本文通過鄰接路段間的相對關(guān)系計算兩鄰接匹配對的兼容系數(shù)，并基于各鄰接匹配對的兼容系數(shù)計算所有鄰接路段對其的綜合影響，即支持系數(shù)，然后根據(jù)支持系數(shù)不斷更新原概率矩陣，直到概率矩陣變化量收斂于某一極小值；③計算各候選匹配對的結(jié)構(gòu)相似性，結(jié)構(gòu)相似性反映某一匹配對及其鄰接路段間的整體匹配程度，之后基于結(jié)構(gòu)相似性設(shè)定一定規(guī)則選取穩(wěn)健匹配集，識別1∶1和1∶M，并對1∶M進行匹配增長，綜合得到最后匹配集。

圖1 本文算法流程Fig.1 Flow chart of the proposed method

2.1 概率矩陣初始化

假設(shè)待匹配兩份數(shù)據(jù)data1＝｛N1，E1｝，data2＝｛N2，E2｝，其中N代表端點集，E代表路段集。本文首先通過緩沖分析確定data1與data2之間的候選匹配路段，再計算候選路段間的差異性指標(biāo)并綜合得到初始的概率矩陣。差異性指標(biāo)是描述兩候選路段間差異的指標(biāo)，評價線性目標(biāo)的差異性指標(biāo)有距離、方向、長度、語義、拓?fù)涞龋?－7］，本文采用距離、方向、長度衡量路段間差異特征。對一種差異性指數(shù)t（t＝dis，dir，len，分別表示距離，方向，長度差異性指標(biāo)），根據(jù)文獻［10］按公式（1）計算路段j∈N2匹配到路段i∈N1的概率pti（i，j）。當(dāng)j＝ －1時，pti（i，－1）表示路段i匹配到空（1∶0）的概率

式中，CSi為路段i的候選匹配路段集；dt（i，j）為路段i和j間差異性指標(biāo)t的值；β1t為關(guān)于差異性指標(biāo)t從data2到data1的誤差因子，用于計算1∶0的匹配概率，文獻［13］未考慮1∶0匹配關(guān)系，導(dǎo)致獨立路段誤匹配至最近的路段。本文通過計算兩份數(shù)據(jù)間Hausdorff距離估算誤差因子，公式為

同理，計算從data1到data2的誤差因子βt2和路段i匹配到路段j的概率ptj（i，j）。實際上，匹配概率應(yīng)為兩路段相互匹配的概率，即i匹配到j(luò)且j也匹配到i的聯(lián)合概率。因此，需將不同方向匹配的概率綜合為相互匹配概率，公式如下

候選匹配對（i→j）的概率定義為不同匹配方向概率的調(diào)和平均數(shù)［15］。i沒有匹配路段的概率pt（i，－1）定義為i匹配到空且所有候選路段不匹配到i的概率，同理，j沒有匹配路段的概率為pt（－1，j）。

為綜合多種差異指標(biāo)［6］，本文分別計算3種差異指標(biāo)的匹配概率，加權(quán)平均得到綜合的初始概率。距離、方向、長度的權(quán)值分別為0.3，0.4，0.3。

2.2 概率矩陣迭代更新

初始概率矩陣確定以后，需計算各匹配對的鄰接路段對其的兼容系數(shù)以不斷更新原概率矩陣。兼容系數(shù)反映了鄰接匹配對間的兼容程度，而通常鄰接匹配對不止一個，需要綜合這些鄰接匹配對的影響，即支持系數(shù)來改進初始概率矩陣，使之更準(zhǔn)確地反映整個路網(wǎng)的匹配關(guān)系。另外，考慮到路段匹配可能存在部分匹配，即1∶M，而根據(jù)一定的距離或長度比閾值很難有效識別這類型匹配。如圖2（a），路段a1和b1的長度，起終點距離皆相同，根據(jù)文獻［13］將被誤識別為1∶1匹配，而實際匹配應(yīng)為（a1→b1，b2）。因此，需要考慮部分匹配對自身的兼容性。

圖2 部分匹配與兼容系數(shù)計算Fig.2 Partial matching and calculation of the compatibility coefficient

基于文獻［12］，本文考慮鄰接路段的相對位置，方向計算兼容系數(shù)，一定程度上克服數(shù)據(jù)本身的旋轉(zhuǎn)和縮放誤差。如圖2（b）所示，計算鄰接匹配（h→k）對（i→j）的兼容系數(shù)C（i，j；h，k），公式為

式中，ratio為點o到e的路段長度與o′到e′之間長度的比值；δdis和δdir分別為路段i、h與j、k之間的位置和方向一致性指標(biāo)，計算公式如下

d1和d2分別為點o、o′和e、e′的距離；α（β）為路段i與h（j與k）的夾角，β1dir（β2dir）與β1dis（β2dis）為公式（2）計算的位置和方向誤差因子。對于自兼容性，如C（i，j；i，k），首先在路段i上尋找中間節(jié)點o′的對應(yīng)點o，再按公式（5）計算兼容性，同理計算C（i，j；h，j）

而某一候選匹配對的概率受到多個鄰接匹配對的影響，需要綜合所有鄰接匹配對的影響以計算支持系數(shù)。如圖3所示，對于候選匹配對（i→j），F(xiàn)1，F(xiàn)2和T1，T2分別為對應(yīng)起點和終點，相應(yīng)的，SF、QF和ST、QT為i和j在起點和終點處的鄰接路段集。如表1所示，不同匹配方向的支持系數(shù)皆由起點、終點處支持系數(shù)構(gòu)成，q1和q2為（i→j）是完全匹配和部分匹配假設(shè)下的支持系數(shù)，本文取兩者中較大值。

圖3 計算支持系數(shù)Fig.3 Calculating the support values

表1 支持系數(shù)計算Tab.1 Calculation of sub－support values

例如，當(dāng)j匹配到i時，在起點的支持系數(shù)qi，F(xiàn)為

式中，r為迭代次數(shù)；C為計算的兼容系數(shù)；p（r）為第r次迭代時概率。同理可計算在終點處的支持系數(shù)qi，T以及i匹配到j(luò)時在起點、終點處的支持系數(shù)qj，F(xiàn)、qj，T。將起終點處支持系數(shù)相加得到j(luò)匹配到i以及j匹配到i的支持系數(shù)q（r）i和q（r）j，再按公式（7）綜合得到最后的支持系數(shù)q（r）。

ηi為i匹配到空的先驗值，計算公式為

最后根據(jù)文獻［11］提出的迭代機制重新計算概率，以路段j匹配到i為例，公式為

根據(jù)不同匹配方向計算的概率按類似公式（3）綜合得到新的概率矩陣。如此迭代，直到概率矩陣變化量小于某一給定閾值ε。

2.3 匹配對選取

為有效克服誤匹配和漏匹配，本文設(shè)定以下3步用于匹配對的選取。

第1步：計算結(jié)構(gòu)相似性。結(jié)構(gòu)相似性反映某一匹配對及其鄰接路段間的整體匹配程度，本文定義候選匹配對（i→j）的結(jié)構(gòu)相似性為

式中，strF、strT分別為對應(yīng)起終點處的結(jié)構(gòu)相似性。假設(shè)SF、QF分別為i和j在起點處的鄰接路段集，根據(jù)二分圖匹配算法［16］計算兩者間的最大匹配，strF即為該匹配組合下的概率之和，同理計算strT。為減少計算量，本文對迭代后概率矩陣按大律法［17］對行和列二值化，兩者相加將概率矩陣分為3級，本文排除第3級候選匹配對，只計算前兩級的結(jié)構(gòu)相似性。

第2步：選取穩(wěn)健匹配對。若（i→j）滿足等式（11），則（i→j）為穩(wěn)健匹配對，添加到隊列queue中。

之后，需逐一檢測queue中每一匹配對與鄰接匹配對是否兼容。如圖4所示，當(dāng)兩鄰接匹配對Mn，Mp的中間結(jié)點（a2，b2）不同時是Mn或Mp的對應(yīng)端點，則認(rèn)為Mn、Mp相互不兼容。對于不兼容的兩鄰接匹配對，本文保留結(jié)構(gòu)相似性大的，刪除結(jié)構(gòu)相似性小的。

圖4 鄰接匹配對之間的兼容性檢測Fig.4 Detecting the compatibility between two neighbouring matching pairs

第3步：檢測1∶M匹配對。對于queue中每一匹配對，判斷其對應(yīng)起點和終點是否皆匹配，如果匹配，則為1∶1對應(yīng)，否則應(yīng)為1∶M匹配。判斷匹配對（i→j）在某對應(yīng)端點處是否匹配，應(yīng)滿足以下條件：①若（i→j）在對應(yīng)端點處的結(jié)構(gòu)相似性大于i和j的其他匹配對在該對應(yīng)端點處的結(jié)構(gòu)相似性；②（i→j）在該對應(yīng)端點處的鄰接匹配對中至少存在一匹配對M∈queue。

若（i→j）在某對應(yīng)端點處滿足條件①和②，則認(rèn)為（i→j）在該對應(yīng)端點處是匹配的。如果（i→j）的對應(yīng)起點和終點皆匹配，則該匹配對為1∶1匹配，否則為1∶M的匹配。逐一檢測queue中的每一匹配對是否1∶M的匹配對，若是，則添加到臨時隊列temp。

逐一遍歷temp中的每一匹配對，如圖5中（i→j），若起點不匹配，則從短路段i的起點開始，遍歷其鄰接路段集G＝｛g0，…，gk，…，gn｝，從中選擇與j的結(jié)構(gòu)相似性最大的gk，添加新的匹配對（gk→j）到queue，再從gk的鄰接路段中尋找與j的結(jié)構(gòu)相似性最大的路段，組成新的匹配對添加到queue中，如此循環(huán)直到找到j(luò)的起點的匹配端點。若（i→j）的對應(yīng)終點也不匹配，則進行同樣的匹配增長過程。進行完匹配增長后，對新的queue中每一匹配執(zhí)行第2步中的兼容性檢查，同樣保留結(jié)構(gòu)相似性大的，刪除結(jié)構(gòu)相似性小的。

圖5 匹配增長過程Fig.5 Matching growing from unmatched nodes

3 試驗與評價

本文選取中國武漢，瑞士蘇黎世的兩份路網(wǎng)數(shù)據(jù)進行匹配算法的驗證。武漢為四維圖新（NavInfo?）導(dǎo)航數(shù)據(jù)與OSM數(shù)據(jù)匹配，蘇黎世為NAVITEQ?導(dǎo)航數(shù)據(jù)與OSM數(shù)據(jù)匹配，空間參考為WGS－84，各數(shù)據(jù)的節(jié)點、路段統(tǒng)計如表2所示。導(dǎo)航數(shù)據(jù)格式為＊.shp，OSM數(shù)據(jù)格式為＊.osm，首先利用FME將OSM數(shù)據(jù)轉(zhuǎn)換為統(tǒng)一的＊.shp，并通過 Microsoft Visual Studio 2008（C＃）＋ArcGIS Engine 9.3實現(xiàn)本文匹配算法，其中武漢、蘇黎世的緩沖距離分別設(shè)為200m和40m，迭代終止閾值設(shè)為0.000 5。

表2 路網(wǎng)數(shù)據(jù)統(tǒng)計描述Tab.2 The statistical description of road networks

試驗的部分結(jié)果如圖6所示，黑實線為導(dǎo)航數(shù)據(jù)，灰實線為OSM數(shù)據(jù)，虛線為匹配路段的中點連線，其中圖6（a）、（b）為武漢地區(qū)試驗結(jié)果，圖6（c）、（d）為蘇黎世地區(qū)試驗結(jié)果。由圖6（a）可見，對于武漢地區(qū)，OSM數(shù)據(jù)由于配準(zhǔn)精度不高，與專業(yè)的導(dǎo)航電子數(shù)據(jù)的偏差達上百米，而本文方法可以有效地克服這種偏差實現(xiàn)整體匹配最優(yōu)，圖中虛線橢圓中匹配的兩條路段位置不是最近的；圖6（b）說明本文方法能夠較好地匹配復(fù)雜路網(wǎng)結(jié)構(gòu)。圖6（c）中，OSM數(shù)據(jù)在某些區(qū)域要比導(dǎo)航數(shù)據(jù)詳細(xì)，將OSM數(shù)據(jù)匹配到導(dǎo)航數(shù)據(jù)或反過來，得到的匹配結(jié)果不同，本文方法不受匹配方向影響，有效避免獨立道路誤匹配到另一份數(shù)據(jù)中最相似道路。圖6（d）表明本文能有效處理眾源數(shù)據(jù)中可能存在結(jié)構(gòu)不完整的路網(wǎng)，虛線橢圓中OSM路段由于用戶編輯錯誤導(dǎo)致連通道路斷開。

圖6 試驗匹配結(jié)果Fig.6 The matching results

為定量評價本文算法的匹配精度，通過與人工匹配對比，計算匹配的準(zhǔn)確率P和召回率R公式為

式中，TP為正確匹配數(shù)；FP為錯誤匹配數(shù)；AM為人工無法判斷的匹配數(shù)；FN為漏匹配數(shù)。統(tǒng)計結(jié)果如表3所示，兩試驗區(qū)的匹配精度均達到96%以上，召回率也達到于90%以上，但還有5.3%～8.8%的匹配對未能有效被識別。為分析本文算法的效率，表3統(tǒng)計了單次迭代耗時，武漢為4s（8580候選匹配對），蘇黎世為30s（22559候選匹配對），這是由于蘇黎世路網(wǎng)密度比武漢高，落在緩沖區(qū)中的候選路段相對多導(dǎo)致計算量大。實際上，總匹配耗時受如試驗范圍、數(shù)據(jù)詳細(xì)程度、路網(wǎng)密度、緩沖閾值、收斂條件等多種因素的影響。

表3 匹配結(jié)果統(tǒng)計Tab.3 Statistics of matching results using the proposed method

如圖7所示，通過緩沖分析確定灰色的路段a1、a2的候選路段，分別統(tǒng)計迭代過程中路段a1、a2與各候選路段的概率變化，如圖8（a）和（b）所示。由圖8（a）可見，相比a1的其他候選路段，a1與c2的初始概率最高，且在迭代過程中仍保持較高概率，并收斂于某一穩(wěn)定值；由圖8（b）可見，雖然a2與d1、d2的初始概率最高，但通過不斷迭代，a2與d4、d5的概率達到最高并收斂于某一穩(wěn)定值。結(jié)果表明，本文方法結(jié)合鄰接路段的匹配關(guān)系對概率矩陣進行迭代更新，可以幫助尋找正確的匹配路段，并有效降低誤匹配率。

圖7 迭代過程舉例Fig.7 The matching examples for relaxation

為分析本文方法對參數(shù)的敏感度，本文設(shè)定5組緩沖距離（100m，150m，200m，300m，and 400m）對武漢地區(qū)進行多次試驗，試驗結(jié)果中的正確匹配，錯誤匹配，模糊匹配數(shù)目統(tǒng)計如圖9所示。由圖可見，閾值設(shè)為150m至400m之間時，匹配精度大約為95.4%～97.2%，而當(dāng)閾值為100m時，精度降為67.4%，說明本文方法一定程度上不受閾值的影響，可以達到95%以上，但使用較小的閾值導(dǎo)致不能完全檢測可能的匹配路段，使匹配精度降低。為取得較高匹配精度，建議使用稍大的閾值，然而，過大的閾值將導(dǎo)致計算時間的增加。

圖8 迭代過程中概率變化Fig.8 Matching probability change during iteration

圖9 不同緩沖閾值下的精度比較Fig.9 Matching precision of different buffering thresholds

4 結(jié) 論

本文方法考慮自身的幾何相似性以及鄰接候選匹配路段的兼容性，迭代計算候選路段的匹配概率，最后基于收斂的概率矩陣選取最后的匹配結(jié)果。試驗結(jié)果表明，該方法對偏差較大的路網(wǎng)數(shù)據(jù)能達到較高的匹配精度，不存在匹配方向性問題，且能夠識別1∶0和1∶M的匹配關(guān)系，該方法適用于偏差較大、存在一定尺度差異的路網(wǎng)數(shù)據(jù)，可用于存在不完整拓?fù)涞谋娫绰肪W(wǎng)與專業(yè)數(shù)據(jù)的集成，同時為路網(wǎng)更新和變化發(fā)現(xiàn)提供了關(guān)鍵技術(shù)支撐。未來的研究將側(cè)重于計算效率的改進，基于匹配關(guān)系的變化提取和數(shù)據(jù)更新。

［1］ SAALFELD A J.Conflation：Automated Map Compilation［J］.International Journal of Geographical Information Systems，1988，2（3）：217－228.

［2］ ZHANG Meng.Methods and Implementations of Roadnetwork Matching［D］.Munich：Technical University of Munich，2009.

［3］ HU Yungang，CHEN Jun，ZHAO Renliang，et al.Matching of Roads under Different Scales for Updating Map Data［J］.Geomatics and Information Science of Wuhan University，2010，35（4）：451－456.（胡云崗，陳軍，趙仁亮，等.地圖數(shù)據(jù)縮編更新中道路數(shù)據(jù)匹配方法［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2010，35（4）：451－456.）

［4］ RUIZ J J，ARIZA F J，UREA M A，et al.Digital Map Conflation：A Review of the Process and a Proposal for Classification［J］.International Journal of Geographical Information Science，2011，25（9）：1439－1466.

［5］ COBB M A，CHUNG M J，F(xiàn)OLEY H I，et al.A Rulebased Approach for the Conflation of Attributed Vector Data［J］.GeoInformatica，1998，2（1）：7－35.

［6］ VOLZ S.An Iterative Approach for Matching Multiple Representations of Street Data［C］∥Proceedings of the ISPRS Workshop on Multiple Representation and Interoperability of Spatial Data.Hanover：ISPRS，2006：101－110.

［7］ MUSTIéRE S，DEVOGELE T.Matching Networks with Different Levels of Detail［J］.Geoinformatica，2007，12（4）：435－453.

［8］ XIONG D，SPERLING J.Semiautomated Matching for Network Database Integration ［J］.ISPRS Journal of Photogrammetry ＆Remote Sensing，2004，59（1）：35－46.

［9］ SAMAL A，SETH S，CUETO K.A Feature－based Approach to Conflation of Geospatial Sources［J］.International Journal of Geographical Information Science，2004，18（5）：459－489.

［10］ SAFRA E，KANZA Y，SAGIV Y，et al.Location－based Algorithms for Finding Sets of Corresponding Objects over Several Geo－spatial Data Sets［J］.International Journal of Geographical Information Science，2010，24（1）：69－106.

［11］ ZHAO Dongbao，SHENG Yehua.Research on Automatic Matching of Vector Road Networks Based on Global Optimization［J］.Acta Geodaetica et Cartographica Sinica，2010，39（4）：416－421.（趙東保，盛業(yè)華.全局尋優(yōu)的矢量道路網(wǎng)自動匹配方法研究［J］.測繪學(xué)報，2010，39（4）：416－421.）

［12］ SONG W，KELLER J M，HAITHCOAT T M，et al.Relaxation－based Point Feature Matching for Vector Map Conflation［J］.Transactions in GIS，2011，15（1）：43－60.

［13］ ZHANG Yunfei，YANG Bisheng，LUAN Xuechen.Automated Matching Road Networks Based on Probabilistic Relaxation［C］∥Proceedings of ISPRS Workshop on Dynamic and Multi－dimensional GIS.Shanghai：ISPRS，2011：75－80.

［14］ CHRISTMAS W J，KITTLER J，PETROU M.Structural Matching in Computer Vision Using Probabilistic Relaxation［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1995，17（8）：749－764.

［15］ HUNTINGTON E V.Sets of Independent Postulates for the Arithmetic Mean，the Geometric Mean，the Harmonic Mean，and the Root－mean－square［J］.Transactions of the American Mathematical Society，1927，29（1）：1－22.

［16］ MUNKRES J A.Algorithms for Assignment and Transportation Problem［J］.Society for Industrial and Applied Mathematics，1957，5（1）：32－38.

［17］ OSTU N.Threshold Selection Method from Gray－scale Histogram［J］.IEEE Transactions on System，Man，Cybernet，1979，9（1）：62－66.