高永偉，周曉光，任常青

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083；2.國家測繪地理信息局 第一航測遙感院，陜西 西安 710054)

0 引 言

地表覆蓋是指地球表面各種物質(zhì)類型及其自然屬性與特征的綜合體，包括自然形成物及人為建造物等[1-3]。地表覆蓋的分布對(duì)地球水循環(huán)、大氣循環(huán)、碳循環(huán)以及熱循環(huán)有重要的影響。我國于2014年研制出2000年、2010年兩期全球30 m地表覆蓋分類數(shù)據(jù)集,將地表覆蓋分為耕地、水體、人造地表等10大類。為保持其現(xiàn)勢性需對(duì)地表覆蓋數(shù)據(jù)進(jìn)行持續(xù)更新。目前，地表覆蓋更新數(shù)據(jù)主要來源于遙感影像、實(shí)測數(shù)據(jù)等[4-5]。而用此方法提取全球地表覆蓋增量數(shù)據(jù)，存在費(fèi)用高、不及時(shí)等問題。隨著志愿者地理信息的發(fā)展，眾源數(shù)據(jù)成為地理信息的重要來源之一。

眾源地理數(shù)據(jù)是由大量非專業(yè)人員自愿獲取并通過互聯(lián)網(wǎng)向大眾或相關(guān)機(jī)構(gòu)提供的一種開放地理空間數(shù)據(jù)[6-8]。眾源地理數(shù)據(jù)主要來源于公共數(shù)據(jù)、GPS數(shù)據(jù)、VGI網(wǎng)站數(shù)據(jù)以及Web2.0技術(shù)催生下的其他地理數(shù)據(jù)。如OpenStreetMap（OSM）是一個(gè)由地圖制作愛好者組成的社區(qū)網(wǎng)站，用戶可以編輯地圖數(shù)據(jù)等；Wikiloc是GPS軌跡分享的網(wǎng)站。許多學(xué)者對(duì)眾源地理數(shù)據(jù)進(jìn)行了研究，發(fā)現(xiàn)其具有現(xiàn)勢性強(qiáng)、信息豐富、成本低、數(shù)據(jù)量大、質(zhì)量差異大等特點(diǎn)[8-10]；也有學(xué)者提出了對(duì)VGI數(shù)據(jù)質(zhì)量評(píng)價(jià)方法[11-14]及模型轉(zhuǎn)換的規(guī)則[15]。眾源軌跡數(shù)據(jù)是通過裝載GPS記錄儀的設(shè)備記錄的用戶位置數(shù)據(jù)。許多網(wǎng)站提供了大量的軌跡文件，如Wikiloc、OSM網(wǎng)站，軌跡文件數(shù)量達(dá)到5 512 811、2 337 470條（截止到2017年1月14日）；紐約市出租車及轎車委員會(huì)（New York City Taxi and Limousine Commission,簡稱NYC）提供了2009年至2015年全市出租車上下車位置數(shù)據(jù)等。國內(nèi)外眾多學(xué)者對(duì)軌跡數(shù)據(jù)進(jìn)行了研究，主要是路網(wǎng)提取、用戶行為等方面[16]。Zhang、歐陽鴻等基于GPS軌跡數(shù)據(jù)提取了城市道路、支路與小路網(wǎng)等[17-18]；孔慶杰、蔣益娟等通過對(duì)軌跡數(shù)據(jù)柵格化后提取道路網(wǎng)[19-20]；這些可以對(duì)現(xiàn)有路網(wǎng)進(jìn)行補(bǔ)充與更新。Zhou C等研究了手機(jī)數(shù)據(jù)與新加坡土地分類的關(guān)系[23]。而對(duì)于如何從軌跡數(shù)據(jù)中提取人造地表覆蓋數(shù)據(jù)的研究較少。鑒于軌跡數(shù)據(jù)與人造地表覆蓋分布有較高的重合度，本文在借鑒路網(wǎng)密度概念的基礎(chǔ)上，研究了采用眾源軌跡數(shù)據(jù)提取人造地表覆蓋增量信息的方法。

1 總體思路

眾源軌跡數(shù)據(jù)因采集的設(shè)備、人群、公司不同提供的軌跡文件的內(nèi)容格式也不相同。在軌跡數(shù)據(jù)研究中，軌跡點(diǎn)包含的信息由專業(yè)人員根據(jù)需要而定。如在Schroedl等基于K-means 算法以及最小二乘法提取道路網(wǎng)及道路中心線[16]、Zhang 等通過融合大量的軌跡點(diǎn)來更新已有路網(wǎng)信息[12]，采集的軌跡點(diǎn)有經(jīng)度、緯度、時(shí)間、速度、方向等屬性。而眾源軌跡數(shù)據(jù)，軌跡點(diǎn)的信息并不如此完善。表1展示了OSM、Wikiloc、NYC3個(gè)網(wǎng)站提供的軌跡點(diǎn)文件及屬性信息。

表1 眾源軌跡數(shù)據(jù)軌跡點(diǎn)文件及屬性表Tab.1 Crowd sourcing trajectory data and it's attribute

為適用于眾源軌跡數(shù)據(jù)處理，本文采用圖1所示思路對(duì)人造地表覆蓋增量數(shù)據(jù)進(jìn)行提取。即先經(jīng)過數(shù)據(jù)預(yù)處理：對(duì)共有屬性如經(jīng)度、緯度、時(shí)間等提取，并對(duì)軌跡點(diǎn)進(jìn)行投影；為簡化處理及計(jì)算，按照一定規(guī)則將軌跡點(diǎn)柵格化處理；接著根據(jù)掃描半徑及軌跡點(diǎn)柵格覆蓋比率提取人造地表覆蓋數(shù)據(jù)；去除粗提取數(shù)據(jù)中的微小空洞及獨(dú)立柵格單元后，經(jīng)過重采樣與柵格運(yùn)算提取人造地表覆蓋增量數(shù)據(jù)。

圖1 基于眾源軌跡數(shù)據(jù)的人造地表覆蓋增量信息提取總體思路Fig.1 The general idea of change-only artificial surface information extraction using crowd sourcing trajectory data

2 人造地表覆蓋增量數(shù)據(jù)提取方法

由于眾源軌跡數(shù)據(jù)來源不同，軌跡點(diǎn)屬性、格式不盡相同，圖2是需要進(jìn)行預(yù)處理提取的經(jīng)度、緯度、時(shí)間等共有屬性。軌跡數(shù)據(jù)大部分采用WGS-84坐標(biāo)系。為方便后續(xù)距離及面積運(yùn)算，將地理坐標(biāo)轉(zhuǎn)為投影坐標(biāo)，采用Web墨卡托投影坐標(biāo)系。

圖2 眾源軌跡數(shù)據(jù)格式圖Fig.2 Examples for crowd sourcing trajectory data

2.1 軌跡點(diǎn)柵格化

孔慶杰等[14]在基于GPS軌跡點(diǎn)的路網(wǎng)提取中將軌跡點(diǎn)柵格化生成柵格路網(wǎng)，然后進(jìn)行矢量化路網(wǎng)提取。本文借鑒其柵格化思想，增加軌跡點(diǎn)過濾規(guī)則進(jìn)行柵格化處理，具體如下：

首先根據(jù)處理軌跡點(diǎn)的范圍確定柵格化圖層范圍，依據(jù)柵格單位分辨率將柵格圖層分割成M×N個(gè)單元；接著確定柵格單元屬性值。柵格單元屬性值的確定原則：對(duì)于柵格單元(i,j)，統(tǒng)計(jì)其軌跡點(diǎn)數(shù)量Pcount，若Pcount大于最少軌跡點(diǎn)數(shù)Pmin,則將其單元格的屬性值設(shè)為1；否則設(shè)為0.如圖3所示，當(dāng)Pmin=1時(shí)，Pcount＞Pmin=1的單元有左上、左下和右下角3個(gè)，將其值均設(shè)置為1，其余單元格值設(shè)為0。

圖3 軌跡點(diǎn)柵格化原理圖Fig.3 Schematic of trajectory points rasterizing

柵格單元邊界和柵格范圍邊界上的軌跡點(diǎn)P，用CellUp, CellDown, CellLeft, CellRight表示P的左上，左下，右上，右下4個(gè)單元格，按如下規(guī)則處理：

圖4 邊界點(diǎn)處理原則Fig.4 Principles for boundary points processing

規(guī)則1：若ECellUp；則P∈CellUp；否則，P∈CellDown；

規(guī)則2：若ECellLeft，則P∈CellLeft；否則，則P∈CellRight.

2.2 人造地表覆蓋數(shù)據(jù)粗提取

在軌跡點(diǎn)柵格化后，我們即可進(jìn)行人造地表覆蓋數(shù)據(jù)的初步提取。在《城市道路交通規(guī)劃設(shè)計(jì)規(guī)范》(GB50220-95)(簡稱《規(guī)范》)中，路網(wǎng)密度等于某一計(jì)算區(qū)域內(nèi)所有道路的總長度與區(qū)域總面積之比。由于通過軌跡點(diǎn)提取道路算法復(fù)雜，且本文目的為提取人造地表覆蓋數(shù)據(jù)，無需提取道路網(wǎng)。因此，本文在不提取道路前提下采用軌跡點(diǎn)柵格覆蓋比率(用γ表示)進(jìn)行粗提取。γ定義如下：

定義1：γ等于計(jì)算區(qū)域內(nèi)擁有軌跡點(diǎn)的所有柵格單元面積與該區(qū)域面積之比。

式（1）為γ的計(jì)算公式，式中，β為柵格單元屬性值為1的面積，Area為掃描面積；式（2）為β的計(jì)算公式，式中，Ru為柵格單元分辨率，Piexcount為屬性值為1的單元格個(gè)數(shù)；式（3）為Area的計(jì)算公式，式中，R為掃描半徑。

按照上述公式計(jì)算圖3a區(qū)域柵格化得到圖3b的γ值：在圖3b中，軌跡點(diǎn)柵格化單元個(gè)數(shù)為3，總柵格單元個(gè)數(shù)為9，所以圖3b的γ值為33.3%。

人造地表覆蓋數(shù)據(jù)粗提取步驟是：首先以邊長為R的正方形掃描模板對(duì)柵格化軌跡點(diǎn)進(jìn)行掃描，統(tǒng)計(jì)屬性值為1的柵格單元數(shù)Piexcount；計(jì)算掃描區(qū)域的γ值，若γ大于等于軌跡點(diǎn)柵格覆蓋比率提取閥值α，則將該范圍內(nèi)所有柵格單位的屬性值變?yōu)?，否則為0。若區(qū)域面積小于掃描區(qū)域面積，則統(tǒng)計(jì)區(qū)域面積內(nèi)柵格單元為1的個(gè)數(shù)與屬性值為0的個(gè)數(shù)，求其比值，并與α比較，再按照以上原則處理；提取過程如圖5所示,對(duì)于α,若其γ＞=α，則粗提取后為b,否則為c。

圖5 人造地表覆蓋數(shù)據(jù)粗提取示例圖Fig.5 Example of rough extraction for artificial surface data

2.3 微小孔洞及獨(dú)立柵格單元去除

在某些區(qū)域，由于單位面積內(nèi)軌跡點(diǎn)數(shù)量較少，可能導(dǎo)致粗提取后出現(xiàn)圖6所示的孔洞及獨(dú)立柵格單元。對(duì)于孔洞，可按照圖7方式處理：α為中心孔洞；b為凹陷噪聲[17]，繞中心依次旋轉(zhuǎn)90，又有3種類型凹陷；圖8為獨(dú)立柵格單元的處理。具體過程如下：

第一步：首先取3×3的模板掃描柵格數(shù)據(jù)；

第二步：若第(i,j)處單元屬性值為0，則取其相鄰上下左右4個(gè)單元格的屬性值；更對(duì)于a類型的中心孔洞，其上下左右都為1，則直接更改其屬性值為1；對(duì)于b凹陷類型，若周圍有3個(gè)值為1，則第(i,j)處的單元格為1；否則不變；若第(i,j)處的屬性值為1，其上下左右四周為0，則將其屬性值改為0，否則不變；

第三步：若屬性改變的柵格單元數(shù)量大于某一數(shù)量，重復(fù)以上步驟，否則，處理完畢。

圖6 微小孔洞與獨(dú)立柵格單元示例Fig.6 Examples of micro-hole and independent grid unit

圖7 兩種微小孔洞處理方法Fig.7 Methods for micro-hole processing

2.4 人造地表覆蓋增量數(shù)據(jù)提取

經(jīng)過粗提取及孔洞、獨(dú)立柵格去除處理后，可獲取人造地表覆蓋數(shù)據(jù)。由于柵格化分辨率及掃描半徑的選取不同，易造成提取后數(shù)據(jù)分辨率與2010年全球地表覆蓋數(shù)據(jù)集(以下簡稱“2010年數(shù)據(jù)集”)分辨率不一致，后者分辨率為30m。若分辨率不同，則首先對(duì)提取后數(shù)據(jù)采用最鄰近法重采樣處理；接著將數(shù)據(jù)集坐標(biāo)系轉(zhuǎn)為Web墨卡托投影坐標(biāo)系以統(tǒng)一坐標(biāo)系統(tǒng)。最后將提取數(shù)據(jù)與2010年數(shù)據(jù)集進(jìn)行柵格代數(shù)運(yùn)算，提取人造地表覆蓋增量數(shù)據(jù)，如圖9所示。

圖9 人造地表覆蓋增量數(shù)據(jù)提取示例圖Fig.9 Extraction of change-only artificial surfaces data

圖8 獨(dú)立柵格單元處理方法Fig.8 Methods for independent grid unit processing

3 實(shí) 驗(yàn)

為驗(yàn)證本文方法的可行性，以紐約市2015年1月份Yellow出租車軌跡及OSM軌跡數(shù)據(jù)做驗(yàn)證。Yellow數(shù)據(jù)僅包含乘客上下車位置,OSM軌跡是整條軌跡坐標(biāo)。預(yù)處理后，提取經(jīng)度-74.545 296°至-73.232 216°與緯度40.388 878°至 41.248 829°范圍內(nèi)的25 018 112個(gè)出租車軌跡點(diǎn)和1 141 465個(gè)OSM軌跡點(diǎn)。

3.1 軌跡點(diǎn)柵格化實(shí)驗(yàn)

圖10 軌跡點(diǎn)柵格化結(jié)果圖及局部放大圖Fig.10 Trajectory point rasterized result and partial enlargement

本文選取10m分辨率柵格單元，每個(gè)單元的最少軌跡點(diǎn)數(shù)目Pmin=1,按2.2中的方法對(duì)軌跡點(diǎn)進(jìn)行柵格化。結(jié)果如圖10(范圍：經(jīng)度：-74.056°～-73.7.1°；緯度：40.57°～40.926°，3.3節(jié)中范圍與此一致)所示，可以看出軌跡點(diǎn)基本轉(zhuǎn)為柵格圖。軌跡點(diǎn)柵格化后基本沿道路區(qū)域分布。

3.2 人造地表覆蓋數(shù)據(jù)粗提取實(shí)驗(yàn)

在人造地表覆蓋粗提取實(shí)驗(yàn)中，需確定掃描半徑R和軌跡點(diǎn)柵格覆蓋比率提取閥值α。由于2010年全球地表覆蓋數(shù)據(jù)的分辨率為30 m，按照專家經(jīng)驗(yàn)法[24]，R取100 m。據(jù)《規(guī)范》大城市路網(wǎng)密度對(duì)應(yīng)的道路覆蓋率為0.145km2/km2.當(dāng)α取0.15時(shí)，提取結(jié)果如圖11a，效果并不理想。α再依次取0.1,0.05,0.01。結(jié)果如圖11所示。底圖為谷歌影像，紅色為柵格化軌跡點(diǎn)。從圖中可以看出α值越小，人造地表覆蓋數(shù)據(jù)提取效果越好，即：當(dāng)α= 0.01時(shí)，具有較好的提取效果。本文取R=100 m,α= 0.01進(jìn)行下一步實(shí)驗(yàn)。

圖11 R=100時(shí)，各提取閥值α提取結(jié)果圖Fig.11 R = 100, extracted results according to diあerent threshold α

3.3 微小孔洞及獨(dú)立柵格單元去除

人造地表覆蓋數(shù)據(jù)粗提取后，需要去除微小孔洞及獨(dú)立柵格單元。實(shí)驗(yàn)采用3×3的模板，按照2.3節(jié)中的方法進(jìn)行處理，在單元格變化個(gè)數(shù)設(shè)置為50時(shí)，經(jīng)過6次處理，去除獨(dú)立柵格單元12 482個(gè)，微小孔洞單元數(shù)分別為13 358,1 542,356,142,65,29個(gè)，微小孔洞及獨(dú)立柵格單元大部分被去除，如圖12所示。

圖12 孔洞與獨(dú)立柵格單元去除結(jié)果圖Fig.12 Result graph of removed micro-holes and independent grid unit

3.4 人造地表覆蓋增量數(shù)據(jù)提取實(shí)驗(yàn)

提取的人造地表覆蓋數(shù)據(jù)并重采樣后與2010年紐約市人造覆蓋數(shù)據(jù)集作柵格運(yùn)算，提取人造地表覆蓋增量數(shù)據(jù)，如圖13所示。從圖中可以看出，軌跡提取的人造地表覆蓋數(shù)據(jù)主要集中在城市區(qū)域，在局部地區(qū)可以補(bǔ)充2010年人造地表覆蓋未提取區(qū)域。如圖13e所示,藍(lán)色為提取的人造地表覆蓋增量數(shù)據(jù)。

圖13 人造地表覆蓋增量數(shù)據(jù)提取結(jié)果圖Fig.13 Results of change-only artificial surfaces data extraction

4 結(jié)束語

目前，采用遙感影像、測繪數(shù)據(jù)提取全球地表覆蓋變化信息的傳統(tǒng)方式，存在費(fèi)用高，不及時(shí)等問題。而眾源軌跡數(shù)據(jù)具有數(shù)據(jù)共享、準(zhǔn)實(shí)時(shí)性等特點(diǎn)。本文在分析眾源軌跡數(shù)據(jù)特點(diǎn)后研究了基于眾源軌跡數(shù)據(jù)提取人造地表覆蓋增量數(shù)據(jù)的方法：首先將軌跡數(shù)據(jù)柵格化，再將柵格化后的數(shù)據(jù)進(jìn)行粗提取并去除微小孔洞及獨(dú)立柵格單元，最后將軌跡提取人造地表覆蓋數(shù)據(jù)與2010年人造地表覆蓋數(shù)據(jù)進(jìn)行柵格運(yùn)算提取增量數(shù)據(jù)，并以紐約市出租車及OSM軌跡數(shù)據(jù)驗(yàn)證了方法的可行性。但在以下方面還需進(jìn)一步研究：①不同區(qū)域、不同城市掃描半徑及軌跡點(diǎn)柵格覆蓋比率提取閾值確定問題；②軌跡數(shù)據(jù)精度導(dǎo)致的人造地表覆蓋數(shù)據(jù)提取錯(cuò)誤問題；③比較采用路網(wǎng)密度和軌跡點(diǎn)柵格化率作為從軌跡數(shù)據(jù)中提取人造地表覆蓋指標(biāo)的效率和精度等問題。