基于邊緣檢測法的風(fēng)自記紙圖像數(shù)字化

李亞麗，黃少平，鞠曉慧

(1.陜西省氣象信息中心，陜西 西安 710014;2.陜西省氣象局秦嶺和黃土高原生態(tài)環(huán)境氣象重點實驗室，陜西 西安 710016;3.江西省氣象信息中心，江西 南昌 330046;4.國家氣象信息中心，北京 100081)

0 引 言

氣象資料是大氣科學(xué)學(xué)科發(fā)展的基礎(chǔ)，無論對天氣預(yù)報、科學(xué)研究、防災(zāi)減災(zāi)，還是對全球氣候變化問題等國家方針政策的制定，都至關(guān)重要。自古以來，人們便對氣象信息及規(guī)律進行觀測和記載，形成了寶貴的氣候資料，眾多科學(xué)家應(yīng)用這些龐大而系統(tǒng)的氣候資料開展氣候研究，氣候資料的指導(dǎo)性激發(fā)了許多新發(fā)現(xiàn)和對機理的新認識?？萍嫉倪M步推動著氣象觀測的變革，在地面氣象觀測自動化之前，全國2 400多個氣象站自20世紀60年代使用降水、風(fēng)、氣溫等各類自記儀器進行小時乃至分鐘尺度的氣象要素觀測，共積累了40多年約2 000多萬頁的風(fēng)自記紙。

圖像處理和模式識別技術(shù)的發(fā)展，為降水、風(fēng)等氣象要素自記紙?zhí)峁┝擞行У臄?shù)字化處理手段，高時空分辨率的氣候資料必將為探尋極端天氣事件成因以及理解氣候長期變化趨勢原因提供有利條件。王伯民等人采用圖像掃描和曲線跟蹤技術(shù)，研發(fā)降水自記紙彩色掃描數(shù)字化處理系統(tǒng)，通過人機交互的方式提取全國降水自記紙曲線數(shù)據(jù)，建立了中國地面氣象站長序列、高質(zhì)量的分鐘和小時降水文件數(shù)據(jù)集，數(shù)據(jù)成果應(yīng)用于各地的暴雨強度公式的研制，城市內(nèi)澇防治等方面?；诨叶确e分投影的分段線性插值算法對風(fēng)速曲線進行重構(gòu)，基于方向濾波的風(fēng)向風(fēng)速自記紙數(shù)字化算法實現(xiàn)風(fēng)向風(fēng)速自記紙圖像的坐標(biāo)線以及風(fēng)向風(fēng)速曲線的提取和量化，為風(fēng)自記紙圖像的跡線識別與數(shù)據(jù)提取做出了有益的嘗試和探索。

隨著數(shù)字化技術(shù)的發(fā)展，通過計算圖像的某個特征值(梯度或灰度)，然后依據(jù)給定的特征值控制閾值對圖像信息進行處理，實現(xiàn)圖像目標(biāo)邊緣的精確定位成為數(shù)字圖像處理技術(shù)的主要內(nèi)容。Canny邊緣檢測法因其具有定位精度高、信噪比高以及單邊緣響應(yīng)等優(yōu)點得到了廣泛應(yīng)用。但傳統(tǒng)的Canny算子因高低閾值比例固定，且需要根據(jù)檢測圖像的具體情況手動設(shè)置，在實際應(yīng)用中效率較低，難以實現(xiàn)自動實時的檢測。由日本學(xué)者大津提出的最大類間方差法(后簡稱Otsu算法)是一種自適應(yīng)的閾值確定方法，根據(jù)圖像的灰度特性，將圖像分成背景和目標(biāo)兩部分，通過遍歷尋找某一閾值使得背景和目標(biāo)之間的類間方差達到最大值，即最佳分割閾值。將Otsu算法的思想應(yīng)用于Canny算子中，可實現(xiàn)Canny的高低閾值自適應(yīng)選取，從而實現(xiàn)圖像邊緣的自適應(yīng)提取，可以應(yīng)用于需要快速自動提取圖像邊緣的場合。

該文介紹了應(yīng)用基于Otsu算法改進的Canny邊緣檢測方法完成風(fēng)自記紙圖像的傾斜校正、坐標(biāo)線及風(fēng)速風(fēng)向跡線的識別和提取的算法，并根據(jù)《地面氣象觀測規(guī)范》和國標(biāo)《地面氣象觀測規(guī)范-風(fēng)向和風(fēng)速》中EL型電接風(fēng)儀器觀測原理和自記紙記錄整理要求，以及氣象行業(yè)標(biāo)準(zhǔn)《風(fēng)自記紙記錄數(shù)字化 EL型》中風(fēng)自記紙記錄提取的要求，設(shè)計研發(fā)的EL型電接風(fēng)自記紙跡線數(shù)據(jù)提取軟件系統(tǒng)的設(shè)計思路、處理流程和主要功能。該系統(tǒng)于2020年投入全國氣象部門應(yīng)用，目前已完成500多站EL型電接風(fēng)自記紙跡線數(shù)據(jù)的提取。并通過比較地面氣象月報表A/A6/A7文件中人工整理的小時和日最大風(fēng)速風(fēng)向(后文簡稱“原風(fēng)”)與利用EL型電接風(fēng)自記紙跡線數(shù)據(jù)提取軟件系統(tǒng)提取的風(fēng)向風(fēng)速數(shù)據(jù)(后文簡稱“新風(fēng)”)，分析系統(tǒng)識別提取的風(fēng)數(shù)據(jù)的質(zhì)量和精度，以期為新風(fēng)記錄的使用、歷史長序列分鐘/小時風(fēng)數(shù)據(jù)集的建立提供有價值的參考意見。

1 EL型電接風(fēng)自記紙跡線數(shù)據(jù)提取軟件系統(tǒng)

EL型電接風(fēng)自記紙(后文簡稱“風(fēng)自記紙”)是由橫縱坐標(biāo)線組成的網(wǎng)格圖，時間跨度1 480 min，每10 min一條豎線(風(fēng)向區(qū))或弧線(風(fēng)速區(qū))。一張紙上有4條風(fēng)向曲線，每一條風(fēng)向曲線由一條風(fēng)向橫線作為判定風(fēng)向的參照線，某一時刻風(fēng)向的判別由此時刻風(fēng)向豎線相對于風(fēng)向橫線的位置決定，分別用來判斷南/北、東北/西南、東/西、東南/西北，風(fēng)向記錄以方位為單位，按文獻[20]中表1的要求記錄對應(yīng)的風(fēng)向符號。風(fēng)速測量以m/s為單位，最小刻度范圍為1.0 m/s。

根據(jù)EL型電接風(fēng)向風(fēng)速計測量原理，風(fēng)速是按空氣行程200 m，風(fēng)速自記筆尖相應(yīng)跳動1次來記錄的。如10 min內(nèi)跳動1次，風(fēng)速為0.3 m/s(即200 m/600 s)，跳動2次,風(fēng)速為0.7 m/s(即400 m/600 s)，跳動3次在自記紙上移動一格，表示10 min平均風(fēng)速1.0 m/s(即600 m/600 s)。因此，風(fēng)速保留1位小數(shù)，小數(shù)位只能取0、3、7。靜風(fēng)時，風(fēng)速跡線劃平線，風(fēng)速記為0.0 m/s。每2.5 min記錄一次瞬時風(fēng)向。

1.1 系統(tǒng)設(shè)計基本思路

EL型電接風(fēng)自記紙跡線數(shù)據(jù)提取軟件系統(tǒng)設(shè)計的基本思路：利用基于Otsu算法的Canny邊緣檢測方法對風(fēng)自記紙圖像文件進行預(yù)處理、識別并自動跟蹤風(fēng)向風(fēng)速網(wǎng)格和跡線，提取出網(wǎng)格和跡線的坐標(biāo)點等參數(shù)，形成跡線坐標(biāo)數(shù)據(jù)文件，依據(jù)《風(fēng)自記紙記錄數(shù)字化 EL型》中單位時間風(fēng)速風(fēng)向計算算法，將跡線坐標(biāo)數(shù)據(jù)轉(zhuǎn)換生成標(biāo)準(zhǔn)的風(fēng)分鐘/小時數(shù)據(jù)文件。

1.2 系統(tǒng)處理流程

首先，完成包括輸入/出文件路徑配置、提取數(shù)據(jù)線型和顏色配置、提取節(jié)點稀疏閾值設(shè)定、跡線范圍等“系統(tǒng)設(shè)置”，通過“圖像管理”完成臺站信息管理；其次，在對文件名、文件可讀性以及連續(xù)性等進行檢查和圖像傾斜變形校正、圖像增強等預(yù)處理的同時，將掃描的風(fēng)自記紙圖像文件導(dǎo)入系統(tǒng)；接著(圖1)，通過批量自動處理或單張?zhí)幚矸绞剑砑踊虼蜷_圖像文件，設(shè)置圖像網(wǎng)格和自記跡線開始(結(jié)束)時間；識別并自動提取風(fēng)向風(fēng)速網(wǎng)格和跡線，通過人機交互方式回放檢查、修正完成風(fēng)自記紙記錄提取；最后，根據(jù)2.2章節(jié)中風(fēng)速風(fēng)向的計算算法，將提取到的跡線坐標(biāo)數(shù)據(jù)轉(zhuǎn)換生成分鐘、小時數(shù)據(jù)；對數(shù)據(jù)進行完整性、邏輯性檢查后轉(zhuǎn)換生成標(biāo)準(zhǔn)格式的分鐘/小時風(fēng)數(shù)據(jù)文件。

圖1 系統(tǒng)風(fēng)自記紙數(shù)字化處理流程

1.3 系統(tǒng)主要功能

系統(tǒng)包括系統(tǒng)管理、圖像管理、曲線提取和數(shù)據(jù)轉(zhuǎn)換4個模塊，實現(xiàn)的功能主要包括臺站信息、圖像文件和A文件等導(dǎo)入和系統(tǒng)參數(shù)配置，圖像文件預(yù)處理，單張或批量處理的風(fēng)向風(fēng)速網(wǎng)格和跡線的識別提取，跡線識別效果的回放檢查，風(fēng)分鐘/小時數(shù)據(jù)文件轉(zhuǎn)換生成和數(shù)據(jù)質(zhì)量檢查等。同時，通過調(diào)研全國風(fēng)自記紙的共性和特性，對存在的風(fēng)速跡線邁大(小)步(風(fēng)速筆尖正常跳動1次跡線上升或下降超過(或不足)1/3格)、同一張紙內(nèi)風(fēng)速跡線不連續(xù)或時間重疊的跡線中斷等情況，提供操作便捷的人工編輯功能。

2 核心算法

系統(tǒng)在對圖像進行處理的過程中主要采用Canny邊緣檢測算法對圖像進行邊緣檢測；利用Otsu算法計算圖像分割的最佳閾值，以實現(xiàn)圖像邊緣的自適應(yīng)提??；通過Radon變換技術(shù)求解圖像的旋轉(zhuǎn)角度并對圖像進行傾斜校正，這些圖像處理技術(shù)文中不再詳述。

2.1 跡線識別提取算法

2.1.1 圖像預(yù)處理

圖像預(yù)處理主要是對圖像進行尺寸、分辨率調(diào)整，以及糾偏、去污等預(yù)處理。首先，將讀入的彩色風(fēng)自記圖片轉(zhuǎn)成灰度圖；然后使用Canny邊緣檢測方法提取灰度圖中畫面的邊緣，得到網(wǎng)格邊緣的二值圖；最后，通過Radon變換對二值圖進行斜率計算，在(-3,3)的角度范圍內(nèi)旋轉(zhuǎn)，進行圖像校正并填充邊緣。

2.1.2 風(fēng)速風(fēng)向網(wǎng)格提取

以風(fēng)速網(wǎng)格提取為例：第一步，將風(fēng)自記紙圖片進行顏色空間的轉(zhuǎn)換，從BGR轉(zhuǎn)到CMYK。轉(zhuǎn)換后，Y(Yellow)是包含了網(wǎng)格的灰度圖，K(blacK)是包含了跡線的灰度圖，用Y減K(Y-K)去除跡線的干擾，得到網(wǎng)格的灰度圖，提取灰度圖下部的40%為風(fēng)速區(qū)域的灰度圖。第二步，用Otsu算法對風(fēng)速區(qū)域的灰度圖進行二值化，得到網(wǎng)格的二值圖，白點處為網(wǎng)格。第三步，設(shè)某行(列)白點數(shù)最大分別為maxR或maxL，統(tǒng)計每行(列)白點的數(shù)目。如果在某一行(列)，它的白點數(shù)目大于maxR*0.5或maxL*0.5，且比臨近的行(列)的白點數(shù)目多，那么此行(列)可能就是風(fēng)速網(wǎng)格中的某條橫線(縱線)所在的位置。第四步，根據(jù)風(fēng)自記紙的特征，對上面求出的橫線和縱線進行過濾，去除不合理的網(wǎng)格橫線和縱線。當(dāng)橫線和縱線組成的矩形，寬度大于整個自記紙圖片的寬度的60%，高度大于整個圖片的20%，則風(fēng)速網(wǎng)格邊框提取成功。當(dāng)橫線為21條，弧線為148條時，則風(fēng)速網(wǎng)格提取成功。當(dāng)風(fēng)速網(wǎng)格邊框與風(fēng)向網(wǎng)格邊框基本左右對齊，則自記紙網(wǎng)格邊框提取成功。

2.1.3 風(fēng)速風(fēng)向跡線提取

以風(fēng)速跡線提取為例：第一步，將風(fēng)自記紙圖片顏色空間從BGR轉(zhuǎn)換到CMYK，用K-(Y-K)去除網(wǎng)格的干擾，得到跡線的灰度圖，提取灰度圖下部的40%，得到風(fēng)速區(qū)域的灰度圖。第二步，用Otsu算法對風(fēng)速區(qū)域的灰度圖進行二值化，得到風(fēng)速跡線的二值圖，白點處為跡線。第三步，針對風(fēng)速跡線的二值圖，首先去除噪音和不合理的線段；其次進行兩次二值圖膨脹連接間距小的線段；然后連接距離小于

r

的線段(

r

為圖片的寬度/148)，保留寬度大于

s

的線段(

s

為圖片的寬度/24)，刪除其余的線段，得到最終的風(fēng)速跡線二值圖；最后采取張氏細化算法對風(fēng)速線二值圖進行細化，對細化線進行去毛刺處理。第四步，遍歷細化后的風(fēng)速線，記錄坐標(biāo)點列表，然后針對坐標(biāo)點列表依次進行稀疏化、單調(diào)性、峰值谷值等調(diào)整，得到最終的坐標(biāo)點列表。

2.2 風(fēng)速風(fēng)向計算方法

2.2.1 風(fēng) 速

因風(fēng)自記紙年代久遠，部分紙張風(fēng)速跡線模糊，試驗表明，與按人工整理記錄的“數(shù)跳”方法比較，利用“高度差”所計算的風(fēng)速數(shù)據(jù)準(zhǔn)確性更高。因此，系統(tǒng)采用計算高度差的方法計算風(fēng)速，具體如下：通過風(fēng)速跡線與風(fēng)速網(wǎng)格求交，獲取每一時刻的交點坐標(biāo)。以10 min為單位時間，將每10 min的風(fēng)速細化跡線的縱坐標(biāo)位置求差，并取其絕對值，然后除以網(wǎng)格分辨率得到當(dāng)前時刻

T

分鐘的10 min滑動平均風(fēng)速值。即假設(shè)在時間

t

時，交點坐標(biāo)為(

x

,

y

),在

t

前10 min的時刻

t

交點的坐標(biāo)為(

x

,

y

)。設(shè)

t

<

t

<

t

,且

t

時刻對應(yīng)的交點為(

x

,

y

)，那么高度差

h

=max(|

y

-

y

(

t

)|+|

y

-

y

(

t

)|)。設(shè)風(fēng)速網(wǎng)格的高度為

H

，則風(fēng)速網(wǎng)格的高度分辨率為

w

=

H

/20，于是風(fēng)速

v

=

h

/

w

。若小數(shù)部分≥1/6且<1/2，小數(shù)位記為3；若小數(shù)部分≥1/2且<5/6，小數(shù)位記為7；其他情況為0。

2.2.2 風(fēng) 向

與

T

時刻風(fēng)速相對應(yīng)的10 min內(nèi)，以風(fēng)自記紙上提取到的每一條風(fēng)向橫線為時間軸，以風(fēng)向豎線的時間間隔(2.5 min)為步長，對識別出的4條風(fēng)向橫線和多條風(fēng)向豎線，挑取出現(xiàn)次數(shù)最多的風(fēng)向為當(dāng)前時刻

T

分鐘的風(fēng)向，

T

分鐘的風(fēng)向為包括

T

分鐘在內(nèi)的前10 min的最多風(fēng)向，其挑選方法與文獻[20]中挑取各正點10 min最多風(fēng)向的步驟相同。

T

分鐘風(fēng)速為0.0 m/s時，風(fēng)向為靜風(fēng)。

2.2.3 小時和日最大風(fēng)速風(fēng)向

以每小時整點00分的風(fēng)速風(fēng)向數(shù)據(jù)作為該小時的風(fēng)速風(fēng)向。以20時為日界，從每日20:01到次日20:00的分鐘風(fēng)速記錄中比較挑選出最大值作為該日日最大風(fēng)速，并以此風(fēng)速相對應(yīng)的風(fēng)向及時間作為日最大風(fēng)速的風(fēng)向及出現(xiàn)時間。當(dāng)日最大風(fēng)速出現(xiàn)兩次或以上相同時，可任選其中一次的風(fēng)向和時間。日最大風(fēng)速可跨日、跨月、跨年挑取，但只能上跨。

3 應(yīng)用效果分析

隨機挑選榆林、呼中、茫崖、大同、無棣、林芝、萍鄉(xiāng)等7站約15 000張風(fēng)自記紙跡線提取數(shù)據(jù)，將原風(fēng)記錄視為“真值”，計算新風(fēng)與原風(fēng)記錄的小時和日最大風(fēng)速的平均偏差、平均絕對偏差、均方根誤差及風(fēng)向相符率，對比分析兩者差異，評估提取數(shù)據(jù)的精度和質(zhì)量，對評估結(jié)果中差異較大的數(shù)據(jù)通過分類比較、個例分析等方法，查找總結(jié)造成差異的可能原因。文中“偏差”和“絕對偏差”，均為新風(fēng)速減去原風(fēng)速值，并按16方位定量評價風(fēng)向相符率，當(dāng)兩風(fēng)向方位差在±一個方位以內(nèi)時，則認為兩者相符。即：風(fēng)向相符率=風(fēng)向相符次數(shù)/對比總次數(shù)*100%。

3.1 風(fēng)記錄對比分析

通過新風(fēng)與原風(fēng)日最大風(fēng)與小時風(fēng)記錄對比(表1)分析發(fā)現(xiàn)，平均來看，新日最大風(fēng)速大于原日最大風(fēng)速，平均偏差在0.01～0.33 m/s之間，平均絕對偏差在0.28～0.57 m/s之間，均方根誤差在0.43～0.96 m/s之間，風(fēng)向相符率為79.7%～89.6%。小時風(fēng)平均偏差在-0.08～0.01 m/s之間，平均絕對偏差0.13～0.22 m/s之間，均方根誤差在0.23～0.84 m/s之間，風(fēng)向相符率90.1%～99.9%。比較而言，新風(fēng)與原風(fēng)兩者之間差異小時風(fēng)明顯小于日最大風(fēng)，風(fēng)速偏差更小，風(fēng)向相符率更高。

表1 日最大/小時風(fēng)速風(fēng)向?qū)Ρ?/p>以榆林站為例，設(shè)原風(fēng)速為

v

，分

v

<3.0 m/s、3.0≤

v

<7.0 m/s和

v

≥7.0 m/s三個等級對比(表略)發(fā)現(xiàn)，無論日最大風(fēng)還是小時風(fēng)均表現(xiàn)出隨著風(fēng)速本身的增大，新風(fēng)與原風(fēng)兩者間風(fēng)速偏差變大，而風(fēng)向相符率提高的現(xiàn)象。

3.2 差異原因分析

除原風(fēng)風(fēng)速由人工“數(shù)跳”、新風(fēng)利用“高度差”計算所得有所區(qū)別外，對于風(fēng)自記紙記錄，無論是早期人工整理還是現(xiàn)今使用跡線數(shù)據(jù)提取軟件系統(tǒng)自動提取，均容易受到自記跡線模糊、粗細不均、中斷、墨跡污染等紙張質(zhì)量因素的影響。A文件中人工整理所得風(fēng)自記紙記錄因受觀測員個人視力、習(xí)慣、采樣時間不一致等影響，容易造成隨機誤差，而跡線數(shù)據(jù)提取軟件系統(tǒng)采用統(tǒng)一的技術(shù)標(biāo)準(zhǔn)，計算精度高，系統(tǒng)誤差穩(wěn)定可控，數(shù)據(jù)的準(zhǔn)確性和可靠性更高，長時間保存造成的紙張泛黃變色、潮濕或掃描變形對現(xiàn)今的跡線識別提取影響更大，這與文獻[23]中“讀取風(fēng)向錯誤主要與圖像中跡線的清晰程度有較大的關(guān)系，圖像跡線越清晰軟件讀取正確率越高，且正確率優(yōu)于人工讀取”的結(jié)論一致。

1971～2004年榆林站新風(fēng)與原風(fēng)小時風(fēng)速的平均偏差和絕對偏差有隨時間逐漸變小的趨勢(見圖2)，年平均偏差和絕對偏差分別在-0.19～0.0 m/s和0.12～0.23 m/s之間波動,20世紀70年代明顯大于1980年以后，直到2004年風(fēng)自記紙記錄僅作為臺站平行觀測記錄處理，偏差略有變大，也側(cè)面反映了紙張質(zhì)量、人工整理等因素對數(shù)據(jù)質(zhì)量的影響。

圖2 榆林站1971～2004年小時風(fēng)速的平均偏差、 絕對偏差和風(fēng)向相符率

4 結(jié)束語

該文簡要介紹了EL型電接風(fēng)自記紙跡線數(shù)據(jù)提取軟件系統(tǒng)設(shè)計的基本思路、處理流程和主要功能，以及系統(tǒng)所采用的應(yīng)用Otsu算法改進的Canny圖像邊緣檢測方法實現(xiàn)的跡線數(shù)據(jù)提取算法、風(fēng)向風(fēng)速數(shù)據(jù)計算算法，并通過比較A文件中人工整理的風(fēng)速風(fēng)向與利用軟件系統(tǒng)提取的風(fēng)數(shù)據(jù)之間的差異和分析差異原因，評估軟件提取數(shù)據(jù)的質(zhì)量和精度。應(yīng)用效果表明：系統(tǒng)對標(biāo)準(zhǔn)的EL型電接風(fēng)自記紙圖像識別率可達到100%，能夠滿足風(fēng)自記紙數(shù)字化對數(shù)據(jù)質(zhì)量和精度的要求。因?qū)Ρ确治鰯?shù)據(jù)樣本較少，尚不能全面代表風(fēng)自記紙跡線提取數(shù)據(jù)的質(zhì)量，有待在風(fēng)分鐘/小時標(biāo)準(zhǔn)數(shù)據(jù)集的制作過程中進一步進行質(zhì)量評估。