2025年2月13日木曜日

Non-Landslide Sampling

後輩君から非発生のサンプリングについて相談がありました。
個人的には崩壊斜面と非崩壊斜面の差を求めるのが最も厳しい条件なので、低地を除くことを勧めました(相対的に結果は悪化しますが)。前に話したバッファの話題が出なかったので忘れてしまったのかもしれませんが、非発生のサンプリングについて考えて迷ったのは良いことだ思います。

あらためて非発生のサンプリングにかかわる文献を見てみますと、ひどいサンプリングが多い。低地のみから非発生データを取得して見かけの判別精度を上げています。この方法は簡単に精度を上げることができるのですが、使い物になりません。「9割超えの結果はまず疑え」と、どこかで読んだことがありますが、素人さんからみると良い結果に見えるのでしょう。
文献の著者から見ると、ひどいと思う日本人の感覚の方がおかしいのかもしれません。たまたまその斜面が崩壊しなかっただけで、危険性が高いのは変わらない。そこから非崩壊サンプルを取得するのは合理性に欠ける。そもそも、地形量から推定できる崩壊・非崩壊の予測精度はそれほど大きくないですから、著者のように斜面は危ないと割り切ってしまう潔さも必要なのかもしれません。
力学の視点で LSM を作ると精度がよさそうに思えますが、これも含水率や土砂の厚さなど様々な条件を正しく設定できない限り当たらないはずです。この設定は相当困難。そう思うと、地形量のみから作成できる LSM のレベルがそれほど高くないことは容易に想像できるでしょう。達成できる精度を念頭に置きつつ、こだわる場所を選択すべきなのでしょうね。

現状でのサンプリング方法は以下でしょうか。
・ランダムサンプリング < バッファ < 斜面単位
・重心、コア < 頭部 < 滑落崖

非発生のサンプリング方法は重要です。引き続き考えて参りましょう。
以下、文献のメモ(感想付き)です。

感想〇
斜面単位でのサンプリングは、バッファ500m等よりも優位
予測結果をK-meansクラスタリングによりゾーニング。これが最も優れている。

感想X(保守的すぎる予測、機械学習不要)
過小評価されがちな重要な要素は、非発生地域のサンプリング(Rabby et al., 2023)。
LSMの信頼性と精度は、データ選択の卓越性とその形成に採用された方法論と本質的に関連している。
一般的に使用される非発生のサンプリング方法は次の通り。
1)非発生エリアからサンプリング
2)地形特性と専門家の判断によって決定される傾斜閾値よりも低いエリアの使用。
(Adnan et al., 2020)は2度、(Ali et al., 2021)は3度の閾値を使用。
3)地すべりから一定の距離を超えるエリアから選択。
(Taalab et al., 2018)は、200mバッファを使用。

当研究では以下の2種を比較、シナリオ1の方が高いF1ROC-AUCを示した。
1.傾斜角(シナリオ 1): 指定された傾斜角の閾値を下回るエリアに基づくサンプリング方法。
この研究では、10°の傾き閾値を採用。調査範囲内の地すべりのうち <10° の傾斜で発生したのは 2% に過ぎず、このクラスの頻度比はわずか 0.06 である。
この低頻度比は、地すべり発生との相関が極めて弱いことを示している。
2.バッファー(シナリオ 2): 本研究では、250mバッファを利用。

感想△(危険度抵からのみ取得するのは保守的)
LSA手順には、主に地すべり(ポジティブ)および非発生(ネガティブ)データセットのサンプリング、影響要因の決定、モデリングとマッピング、および結果の精度分析が含まれる(Barik et al., 2017)。
主な非発生のサンプリング手法は4つ
1)ランダムにサンプリング(例:Okalp and Akgün、2016;Bueechi et al., 2019;Azarafza et al., 2021)・・・最も一般的
2)バッファ(Xi et al., 2022)
3)自己組織化ニューラルネットワーク(Huang et al., 2017)、類似性に基づくアプローチ(Zhu et al., 2019)
4)傾斜角の小さい地形領域または平野地域からのサンプリング (Kavzoglu et al., 2014;Lucchese et al., 2021;Okalp and Akgün, 2022)
欠点
1) と 2) によって生成されたサンプルは、地すべりの発生しやすい急な斜面に配置されている可能性あり。
4)は地すべりの発生しやすい河川近くに配置される可能性あり。

当研究では以下を比較
SVM からの Very Low (VL) ゾーン
C5.0-DT からの VL ゾーン
LRからのVLゾーン
地すべりゼロエリア
バッファ距離 <200 m
バッファ距離 200–400 m
バッファ距離 >400 m

地すべりのコアのみを正のサンプルとして使用し、地すべりの危険性モデリング中に地すべりの境界をカバーするピクセルを破棄する方法が良い
バッファー距離の増加に伴って AUC の精度が向上
C5.0モデルのVLゾーンの精度が高い。

感想△(危険度抵・高から1:1で取得するのは、エリアにより保守的になる)
非発生サンプルを選択するための統一基準はない[22,23]。
①地すべりのない地域でのランダムサンプリング、地形的特徴に基づく低傾斜サンプリング、②地すべり境界のバッファー外のサンプリング等がある。
これらの方法は主観的な判断や特定の地理的要因に依存することが多く、潜在的な地すべり地域をネガティブサンプルとして誤って選択するリスクを完全に回避することはできない。
入力データが地すべり発生地域と非発生地域の両方の地理的特性を十分に表していることを確認することが必要。
現在の研究では、一般的に発生・非発生を1:1の比率で使用。ただし、このバランスの取れたサンプリング戦略は、実際のフィールド条件を完全には反映しいない[14]。
一部の研究者は、不均衡データを処理するためオーバーサンプリングやアンダーサンプリングなどの戦略を提案。それでも情報の損失やオーバーフィットのリスクに直面している。

当業務では③SHAP値を用いたサンプリングを実施。初期モデルとしてRFを構築。
シナリオ2:予測結果より、エリアを2分(危険度艇・高)、面積比1:4.23でサンプリング
シナリオ3:最も寄与の大きなNDVIのSHAP値分布により、エリアを2分(危険度艇・高)、そこから同数ずつ非発生を選択
過剰適合の可能性あり[14,35]。

感想〇
インベントリに15度未満の発生事例ナシ。
余裕をもって10度以上を対象とした。

2025年2月12日水曜日

Raspberry Pi - PC 直結

ラズパイと PC を LAN ケーブルで直結し操作する方法は、検索にてたくさん引っかかります。そのうちの一つを今日試しました(屋外でラズパイ単体を操作する必要に迫られました)。
Raspberry Pi 4をPCと有線LAN接続して外でも操作できるようにしてみた | DevelopersIO

RasPi
・LANケーブルを接続せずに起動

PC
・iPhone に Wi-Fi 接続
・USBアダプター[イーサーネット4] と LAN ケーブルで RasPi と直結
・iPhone の Wi‐Fiプロパティでネットワークを共有、ホームネットワーク接続[イーサーネット4]
・コマンド arp -a で マックアドレス b8 から始まるインターフェース(RasPi)の IP を確認
・VNC で RasPi に接続

RasPi
・USBメモリを差し込む
・wdisk 11 /media/pi/USB_Name/out


2025年2月6日木曜日

Windows11 への RDP接続

Windows11(Microsoft Entra 参加済み) への RDP接続で、うまくいかないWin PCが出てきました。

ストア版はつながる、OS付属の標準版?デスクトップ版?はダメな状況です。ストア版はマルチモニタに対応していない、再起動するまでキーマップがおかしくなるなど不満あり。だが仕方がない状況なので使っていたのですが、本日こちらもつながり難い状況に陥りました。仕方なく重い腰を上げて調べることにしました。

で、得られた答えはこちら。簡単に出てきました。もっと早く調べたらヨカッタ。

https://learn.microsoft.com/ja-jp/azure/virtual-desktop/rdp-properties

enablecredsspsupport
構文: enablecredsspsupport:i:<value>
説明: クライアントが認証に資格情報セキュリティ サポート プロバイダー (CredSSP) を使用するかどうかを指定します (使用可能な場合)。
サポートされている値:
0: オペレーティング システムが CredSSP をサポートしている場合でも、RDP では CredSSP を使用しません。
1: オペレーティング システムで CredSSP がサポートされている場合、RDP では CredSSP を使用します。
[既定値]: 1
適用対象:
Azure Virtual Desktop
リモート デスクトップ サービス
リモート PC 接続

クライアント側の rdp ファイルに enablecredsspsupport:i:0 を追記すればOK。
接続できました。


2025年2月2日日曜日

H/V の深度表現

H/VからS波構造、あるいはインピーダンス境界深度を求める文献をいくつか読んでいました。

単点観測結果から地下構造を推定する方法は、国内では流行っていませんが、海外では比較的利用されているようでした。一昨年のイタリア出張時でも類似事例を拝見しました。国内で利用されていない理由はわかりませんが、比較的古くからある手法のようです。

大別すると、多層線形地盤を対象にした1次元応答解析ベースの手法と、周波数や速度から疑似深度への変換関数を定義する簡易手法の2種です。両方試しましたが、私の実力では後者でとどめても良いいかという印象です。

前者は地すべり地で取得した30点以上のデータをソースに、OpenHVSRを適用しました。これはガイド付きモンテカルロ法(MC)に基づき最適モデルからランダムに摂動されたモデルを作成し、よりよい構造を探索していくアプローチを採用しています。複数の観測箇所のデータでも一度に処理できるところ、それらを補間して2次元、3次元表示してくれるところが長所です。
残念ながら、私の手元のデータではフィッティングがうまくいきませんでした。もっと現場の数を増やさないと、どこに問題があるのかわからないので、引き続きデータを取得して試してみましょう。

後者の簡易手法はEXCELでも計算できます。複数の文献では、速度と深度の関係を簡単に定式化し、1/4波長則を関連付けて定積分した式を利用していました。これは前者に比べて物理量を表現していない曖昧な絵になります。が、それはデータにも起因しているのだろうなと想像しています。

国内で使われない理由があるのでしょうから、今後はそれも含めて適用性を調べてみましょう。

2025年1月25日土曜日

win format

winformat
https://wwweic.eri.u-tokyo.ac.jp/WIN/man.ja/winformat.html

win file
├1秒目
│├─ブロックサイズ(4B)
│└─秒ブロック(可変長)
│    ├─秒ヘッダー(6B)
│    ├─第1チャネルブロック(1秒分)
│    │  ├─チャネルヘッダー(4B)
│    │  │  ├─チャネル番号 (2B)
│    │  │  ├─サンプルサイズ情報(0.5B=4bit)
│    │  │  └─サンプリングレート(1.5B=12bit)
│    │  └─データブロック 1チャネル・1秒分(可変長)
│    │       ├─先頭サンプル(4B)
│    │       └─差分(または生値) x (サンプリングレート-1)個
│    ├─・・・・
│    └─最終チャネルブロック(1秒分)
├2秒目
├・・・・
└最終秒


●ディスク/テープ上の win フォーマット
[ ブロックサイズ(4B) | 秒ブロック(可変長) ]
 ブロックサイズ(4B) | 秒ブロック(可変長) ]
     ... (秒数分だけ繰り返し) ...

●ブロックサイズ(4Byte)
ex.
1バイト目: 0x00
2バイト目: 0x00
3バイト目: 0x0B
4バイト目: 0xD4

big-endian として無符号整数に変換
val = int.from_bytes(b'\x00\x00\x0b\xd4', byteorder='big')
0x00000BD4 (16進数) = 3028 (10進数) 

0BD4 = 11 × 16^2 + 13 × 16^1 + 4 × 16^0
= 11 × 256 + 13 × 16 + 4 × 1
= 2816 + 208 + 4
= 3028
(3028-4)Byteが秒ブロック長

●秒ブロック
[ 秒ヘッダー(6B)| 第1チャネルブロック(1秒分) | ... | 最終チャネルブロック(1秒分) ]

●秒ヘッダー(6Byte)のBCD(Binary Coded Decimal)形式
ex.
項目 バイト数 BCD格納
1B 0x20 → 「20」 (2020年)
1B 0x10 → 「10」
1B 0x31 → 「31」
1B 0x13 → 「13」
1B 0x59 → 「59」
1B 0x08 → 「08」

●チャネルブロック(1秒分)
[ チャネルヘッダー(4B) | データブロック 1チャネル・1秒分(可変長) ]

●チャンネルヘッダー
[ チャネル番号 (2B) | サンプルサイズ情報(0.5B=4bit) | サンプリングレート(Hz)(1.5B=12bit)]

サンプルサイズ情報 サンプルサイズ(B) 差分or生値
0 0.5 差分値
1 1 差分値
2 2 差分値
3 3 差分値
4 4 差分値
5 4 生値

sample_size_info = (last2 >> 12) & 0x0F  # 上位4bit
sampling_rate    =  last2        & 0x0FFF # 下位12bit

(last2 >> 12) & 0x0F
まず last2 を右に12bitシフトして、最上位4bitを右端へ持ってくる。
& 0x0F で、その4ビットだけを取り出す (残りは 0)。
16進数:0x0F、10進数: 15、2進数: 0000 1111

last2 & 0x0FFF
0x0FFF は下位12ビットが1、それより上位4bitが0のパターン (16進定数)。
これで、last2 の下位12bitだけを取り出す。

●データブロック
[ 先頭サンプル(4B) | 差分(または生値) x (サンプリングレート-1)個 ]
すべて2の補数表示。
int.from_bytes(b, 'big', signed=True) とすることで、バイト列 b は「ビッグエンディアンの符号付き2の補数値」として解釈され、Python の整数 (int) に変換される。


2025年1月24日金曜日

win ファイル読み込み Python

 winファイルの振動データをPythonに取り込む必要に迫られました。

もともとは 既存の GUI プログラムでテキストに変換していたのですが、これだと CSV で保存してから Python で読みに行くといった無駄な時間が発生します。
C の共有メモリに Python でアクセスする方法が最も速いのですが、これは私には荷が重い。

ということで、刻々書き出される win ファイルを Python で読むことにしました。

といっても、今まで win format を読んだことがありません。確認すると、これがフクザツ。入れ子状態で説明が書かれています。

ということで、GPT-o1 の出番。
win format の解説と wdisk のソースをすべて読ませてから、読み込みスクリプトを作ってもらいました。

フムフム、ナルホド、コウイウコトダッタノカ。

import pandas as pd
from datetime import datetime, timedelta

def bcd_byte_to_int(b: int) -> int:
    """1 バイトBCDを整数に変換 (0x23 -> 23など)"""
    hi = (b & 0xF0) >> 4
    lo = (b & 0x0F)
    return hi * 10 + lo

def decode_bcd_time(bcd6: bytes) -> datetime:
    """
    秒ヘッダー6バイト (BCD: YY, MM, DD, HH, mm, SS) を datetime に変換。
    ここでは 0~49 => 2000~2049, 50~99 => 1950~1999 の例。
    """
    if len(bcd6) < 6:
        return None

    year  = bcd_byte_to_int(bcd6[0])  # 0 ~ 99
    month = bcd_byte_to_int(bcd6[1])  # 1 ~ 12
    day   = bcd_byte_to_int(bcd6[2])  # 1 ~ 31
    hour  = bcd_byte_to_int(bcd6[3])  # 0 ~ 23
    minu  = bcd_byte_to_int(bcd6[4])  # 0 ~ 59
    sec   = bcd_byte_to_int(bcd6[5])  # 0 ~ 59

    # 年の補正 (例: 年<50 => 2000年代, それ以外 => 1900年代)
    if year < 50:
        year_full = 2000 + year
    else:
        year_full = 1900 + year

    try:
        return datetime(year_full, month, day, hour, minu, sec)
    except ValueError:
        # 不正な日時など
        return None

def read_4bit_2s_comp_no_cross(data: bytes, bit_offset: int) -> int:
    """
    data バッファの bit_offset ビット目から 4 ビットを取り出し、
    2 の補数(4bit)として返す (バイト境界は跨がない前提)。
     - 0x0~0x7 => 0~7,  0x8~0xF => -8 ~ -1
    """
    byte_index = bit_offset // 8
    bit_in_byte = bit_offset % 8  # 0~7

    val = data[byte_index] if byte_index < len(data) else 0
    shift_amount = (8 - 4) - bit_in_byte
    nibble = (val >> shift_amount) & 0x0F

    if nibble >= 8:
        nibble -= 16  # => -8 ~ -1
    return nibble


def parse_channel_block(data: bytes, offset: int):
    """
    1つのチャネルブロックをパースして情報を返す。
    戻り値: (消費バイト数, dict(channel=<ch_num>,
                                 sample_size_info=<info>,
                                 sampling_rate=<rate>,
                                 samples=[<int値>, ...]))
    """
    if offset + 4 > len(data):
        return (0, None)
    chan_header = data[offset : offset+4]
    
    # チャネル番号
    ch_num_raw = int.from_bytes(chan_header[0:2], 'big', signed=False)
    ch_num = f"{ch_num_raw:04X}"  # 90EA のような16進大文字表記
    # サンプリングレート等
    last2 = int.from_bytes(chan_header[2:4], 'big', signed=False)
    sample_size_info = (last2 >> 12) & 0x0F
    sampling_rate    = last2 & 0x0FFF
    
    pos = offset + 4
    
    # 先頭サンプル 4バイト
    if pos + 4 > len(data):
        return (4, {
            'channel': ch_num,
            'sample_size_info': sample_size_info,
            'sampling_rate': sampling_rate,
            'samples': [],
        })
    first_sample_bytes = data[pos : pos+4]
    pos += 4
    
    current_value = int.from_bytes(first_sample_bytes, 'big', signed=True)
    samples = [current_value]
    
    diff_count = sampling_rate - 1
    if sample_size_info == 0:
        # 4bit(0.5バイト) 差分
        bit_offset = pos * 8
        for _ in range(diff_count):
            #nibble = read_4bit_2s_comp(data, bit_offset)
            nibble = read_4bit_2s_comp_no_cross(data, bit_offset)
            bit_offset += 4
            current_value += nibble
            samples.append(current_value)
        pos = (bit_offset + 7)//8
        
    elif sample_size_info in [1, 2, 3, 4]:
        # 差分
        sample_len = sample_size_info
        for _ in range(diff_count):
            if pos + sample_len > len(data):
                break
            diff_val = int.from_bytes(data[pos : pos+sample_len], 'big', signed=True)
            pos += sample_len
            current_value += diff_val
            samples.append(current_value)
            
    elif sample_size_info == 5:
        # 4バイト絶対値
        sample_len = 4
        for _ in range(diff_count):
            if pos + sample_len > len(data):
                break
            abs_val = int.from_bytes(data[pos : pos+sample_len], 'big', signed=True)
            pos += sample_len
            current_value = abs_val
            samples.append(current_value)
    
    ch_info = {
        'channel': ch_num,
        'sample_size_info': sample_size_info,
        'sampling_rate': sampling_rate,
        'samples': samples
    }
    consumed = pos - offset
    return (consumed, ch_info)

def read_win_file_as_dataframe(filename: str) -> pd.DataFrame:
    """
    ディスクファイル上の WINフォーマットを読み込み、
    ・行: 「各サンプル」の時刻 (秒ヘッダ + サンプルオフセット)
    ・列: チャネル番号
    ・値: 各サンプルの値
    として1つの DataFrame を返す
    """
    # "各サンプル"を行にするため、時刻単位(秒 + fraction)で紐づける。
    # 複数チャネルを同一時刻にマージするため、以下のような手順:
    #   - data_dict: { pd.Timestamp: { ch_name: sampleVal, ...}, ... }
    #   - 最後に data_dict -> DataFrame へ変換
    
    data_dict = {}
    
    with open(filename, 'rb') as f:
        while True:
            # 1) ブロックサイズ
            size_buf = f.read(4)
            if len(size_buf) < 4:
                break
            block_size = int.from_bytes(size_buf, 'big', signed=False)
            if block_size < 4:
                break
            
            # 2) 秒ブロック本体読み込み
            sec_block = f.read(block_size - 4)
            if len(sec_block) < (block_size - 4):
                break
            
            # a) BCD時刻 (先頭6バイト)
            if len(sec_block) < 6:
                continue
            bcd_part = sec_block[0:6]
            base_dt = decode_bcd_time(bcd_part)  # datetime
            if base_dt is None:
                continue
            
            # b) チャネルブロックを読み込む
            offset = 6
            channel_results = []
            while offset < len(sec_block):
                consumed, ch_info = parse_channel_block(sec_block, offset)
                if not ch_info or consumed <= 0:
                    break
                offset += consumed
                channel_results.append(ch_info)
            
            # c) チャネル結果 ⇒ サンプルの時刻と値を data_dict に格納
            for ch_info in channel_results:
                ch_num = ch_info['channel']
                sr = ch_info['sampling_rate']
                samples = ch_info['samples']
                # 1秒ブロック内で sr 個のサンプルがあると仮定(不足分は途中でbreakしてる場合もあり)
                # 各サンプル i 番は base_dt + i*(1/sr) 秒後とみなす
                #   → Timedelta等で補正
                # カラム名は16進表記でも10進表記でもOK
                ch_name = ch_num
                
                for i, val in enumerate(samples):
                    # サンプル時刻
                    if sr > 1:
                        frac_sec = i / sr
                    else:
                        frac_sec = 0.0
                    
                    sample_dt = base_dt + timedelta(seconds=frac_sec)
                    
                    # data_dict[sample_dt] に ch_name: val を入れる
                    if sample_dt not in data_dict:
                        data_dict[sample_dt] = {}
                    data_dict[sample_dt][ch_name] = val
    
    # 辞書 -> DataFrame
    df = pd.DataFrame.from_dict(data_dict, orient='index')
    # 時刻順にソート
    df.sort_index(inplace=True)
    return df,sr

fname = "input/25012123.36"
df, sr = read_win_file_as_dataframe(fname)


これでA/D値を取り込めます。テストしたファイルは sample_size_info = 2 のみでしたが、これまでのプログラムと出力が一致していることを確認しました。他は追々。

data frame に格納したので、物理量への変換はもちろん指定時刻から前〇〇秒とか、欠測秒をNanで埋めるとか後で容易に加工できます。

コードを読むことでファイル構造を理解するという逆の手順ですが、GPTの性能向上によって可能となりました。ラクナジダイニナッタ。

2025年1月22日水曜日

LSMに使用する地形

手を動かして気づいたのがLSMに使用すべき地形(DEM)です。

LSMで「地すべり地形分布図」のポリゴンを正解として使用する場合、地すべり発生前の地形データは存在しないので、地すべり発生後の地形(DEM)に対し特徴量を抽出することになります。そうすると出来上がるLSMは、過去に地すべりが発生した可能性が高い場所を地形以外の特徴量を加味して危険度として表現したものになります。

一方、LSMで表現したいのは既に発生した場所ではなく、今後地すべりや崩壊が発生する危険度です。この場合、発生した場所での特徴量を発生前の地形等から求め、他の未発生の地域に展開し危険度を表すという作業になります。こちらの方が需要は高いでしょう。

LSM作成の目的によって使用するデータの選択が必要です。
正解データ(すべり・崩壊ポリゴン)と地形(DEM)には、十分に留意が必要です。